Introduction
Despite the progress toward malaria elimination from Latin America over the past 2 decades (Ferreira and Castro, Reference Ferreira and Castro2019), there were 548 000 clinical cases recorded across the region in 2023, with an estimated 139 million people currently at the risk of infection (World Health Organization, 2024). Three-fourths of malaria cases in the Americas are due to Plasmodium vivax and the remainder to Plasmodium falciparum; only 0.1% of the cases are caused by Plasmodium malariae. The Amazon Basin accounts for nearly 90% of the regional malaria burden, and Brazil contributed 29.8% of the laboratory-confirmed cases recorded on this continent in 2024 (World Health Organization, 2024).
Modern humans were already settled in the Americas 21 000 years ago (Pigati et al. Reference Pigati, Springer, Honke, Wahl, Champagne, Zimmerman, Gray, Santucci, Odess, Bustos and Bennett2023), but human malaria parasites are thought to have arrived much later (Bruce-Chwatt, Reference Bruce-Chwatt1965; Carter, Reference Carter2003). Current evidence suggests that human malaria was introduced in the Americas by settlers from southern Europe – mainly Portugal and Spain, where the disease was endemic at the time of the European colonization of the American continent (Bruce-Chwatt and de Zulueta, Reference Bruce-Chwatt and de Zulueta1980) – concurrently with the forced displacement of millions of enslaved people from West and Central Africa into South, Central and North America between the mid-1500s and mid-1800s (Bruce-Chwatt, Reference Bruce-Chwatt1965; Michel et al. Reference Michel, Skourtanioti, Pierini, Guevara, Mötsch, Kocher, Barquera, Bianco, Carlhoff, Coppola Bove, Freilich, Giffin, Hermes, Hiß, Knolle, Nelson, Neumann, Papac, Penske, Rohrlach, Salem, Semerau, Villalba-Mouco, Abadie, Aldenderfer, Beckett, Brown, Campus, Chenghwa, Cruz Berrocal, Damašek, Duffett Carlson, Durand, Ernée, Fântăneanu, Frenzel, García Atiénzar, Guillén, Hsieh, Karwowski, Kelvin, Kelvin, Khokhlov, Kinaston, Korolev, Krettek, Küßner, Lai, Look, Majander, Mandl, Mazzarello, McCormick, de Miguel Ibáñez, Murphy, Németh, Nordqvist, Novotny, Obenaus, Olmo-Enciso, Onkamo, Orschiedt, Patrushev, Peltola, Romero, Rubino, Sajantila, Salazar-García, Serrano, Shaydullaev, Sias, Šlaus, Stančo, Swanston, Teschler-Nicola, Valentin, van de Vijver, Varney, Vigil-Escalera Guirado, Waters, Weiss-Krejci, Winter, Lamnidis, Prüfer, Nägele, Spyrou, Schiffels, Stockhammer, Haak, Posth, Warinner, Bos, Herbig and Krause2024). The arrival of human malaria parasites in the Americas offers a remarkable example of relatively recent host–parasite–vector co-evolution (Rougeron et al. Reference Rougeron, Daron, Fontaine and Prugnolle2022). Parasites encountered new anopheline vector species in the New World – e.g. Anopheles darlingi across most of South America, but mainly in the Amazonian lowlands, and Anopheles albimanus along the Pacific Coast of South America and in Central America, both members of the subgenus Nyssorhynchus. These mosquito species are evolutionarily very distant from the dominant vectors in Africa (Anopheles gambiae complex, subgenus Celia) and southern Europe (Anopheles atroparvus, subgenus Anopheles) to which they had been previously exposed (Molina-Cruz et al. Reference Molina-Cruz, Zilversmit, Neafsey, Hartl and Barillas-Mury2016).
Once in the Americas, P. vivax also adapted to anophelines of the Kerteszia subgenus (Anopheles cruzii, Anopheles bellator and Anopheles homunculus) (de Azevedo et al. Reference de Azevedo, Lorenz, Chiaravalloti-Neto and Sallum2020) and to new vertebrate hosts along the Atlantic coast of Brazil – namely, local platyrrhine primates such as howler monkeys, woolly spider monkeys, capuchin monkeys, uakaris and titis (Duarte et al. Reference Duarte, Fernandes, Silva, Sicchi, Mucci, Curado, Fernandes, Medeiros-Sousa, Ceretti-Junior, Marrelli, Evangelista, Teixeira, Summa, Nardi, Garnica, Loss, Buery, Jr, Pacheco, Escalante, Sallum and Laporta2021; de Oliveira et al. Reference de Oliveira, Rodrigues, Duarte, Rona and Ferreira2021a). It has been proposed that 1 or more spillover or host-shift events gave rise to the P. vivax sister species Plasmodium simium, comprising P. vivax-related parasites that infect monkeys and are also associated with clusters of human malaria cases in the outskirts of major cities in Southeast Brazil, including the metropolitan areas of Rio de Janeiro and São Paulo (e.g., Brasil et al. Reference Brasil, Zalis, de Pina-costa, Siqueira, Júnior, Silva, Areas, Pelajo-Machado, de Alvarenga, da Silva Santelli, Albuquerque, Cravo, Santos de Abreu, Peterka, Zanini, Suárez Mutis, Pissinatti, Lourenço-de-Oliveira, de Brito, Ferreira-da-Cruz, Culleton and Daniel-Ribeiro2017; Duarte et al. Reference Duarte, Fernandes, Silva, Sicchi, Mucci, Curado, Fernandes, Medeiros-Sousa, Ceretti-Junior, Marrelli, Evangelista, Teixeira, Summa, Nardi, Garnica, Loss, Buery, Jr, Pacheco, Escalante, Sallum and Laporta2021). Limited genome sequence data show that P. simium remains nearly identical to P. vivax (Mourier et al. Reference Mourier, de Alvarenga, Kaushik, de Pina-costa, Douvropoulou, Guan, Guzmán-Vega, Forrester, de Abreu, Júnior, de Souza Junior, Moreira, Hirano, Pissinatti, Ferreira-da-Cruz, de Oliveira, Arold, Jeffares, Brasil, de Brito, Culleton, Daniel-Ribeiro and Pain2021; de Oliveira et al. Reference de Oliveira, Rodrigues, Early, Duarte, Buery, Bueno, Catão-Dias, Cerutti, Rona, Neafsey and Ferreira2021b), consistent with a sympatric speciation event at its incipience, likely facilitated by a host switch (de Oliveira et al. Reference de Oliveira, Rodrigues, Duarte, Rona and Ferreira2021a). The process(es) that enable this switch to a new vertebrate host remain unknown. Importantly, while the high genetic similarity between P. vivax and P. simium may lead to misidentification of the causative agent of malaria infections, the rate of such misclassifications is unknown.
Here, an expanded dataset of WGSs of P. simium is explored, including isolates from human malaria cases from Southeast Brazil originally attributed to P. vivax, and from global P. vivax isolates, to further investigate the origins of P. simium and genetic signatures of its adaptation to new hosts.
Materials and methods
Plasmodium simium sequence data
WGS data from 31 P. simium isolates (Mourier et al. Reference Mourier, de Alvarenga, Kaushik, de Pina-costa, Douvropoulou, Guan, Guzmán-Vega, Forrester, de Abreu, Júnior, de Souza Junior, Moreira, Hirano, Pissinatti, Ferreira-da-Cruz, de Oliveira, Arold, Jeffares, Brasil, de Brito, Culleton, Daniel-Ribeiro and Pain2021; de Oliveira et al. Reference de Oliveira, Rodrigues, Early, Duarte, Buery, Bueno, Catão-Dias, Cerutti, Rona, Neafsey and Ferreira2021b; Ibrahim et al. Reference Ibrahim, Manko, Dombrowski, Campos, Benavente, Nolder, Sutherland, Nosten, Fernandez, Vélez-Tobón, Castaño, Aguiar, Pereira, da Silva Santos, Suarez-Mutis, Di Santi, Baptista, Machado, Marinho, Clark and Campino2023) were downloaded from the Short Read Archive (SRA) of the National Center for Biotechnology Information of USA (Supplementary Table 1, Supplementary Materials). The data set comprises 6 isolates from nonhuman primates – 4 from the brown howler monkey Allouata clamitans (3 from the Cantareira Park, São Paulo State and 1 from Guapimirim, Rio de Janeiro State), 1 from a black-fronted titi Callicebus nigrifrons from the Ecological Park of Tietê, São Paulo State, and 1 from a black uakari Cacajao melanocephalus from Guapimirim, Rio de Janeiro State. There were 25 samples of human origin, from Santa Maria de Jequitibá, Espírito Santo State (n = 7), various locations in Rio de Janeiro State (n = 6) and various locations in São Paulo State (n = 12). All sampling sites map to areas that are currently, or were formerly, covered by the Atlantic Forest in Southeastern Brazil (Supplementary Figure 1). The 12 human-derived isolates from São Paulo State were originally described as part of a ‘highly clonal, potential P. simium cluster’ (Ibrahim et al. Reference Ibrahim, Manko, Dombrowski, Campos, Benavente, Nolder, Sutherland, Nosten, Fernandez, Vélez-Tobón, Castaño, Aguiar, Pereira, da Silva Santos, Suarez-Mutis, Di Santi, Baptista, Machado, Marinho, Clark and Campino2023); all harboured 1 or both single-nucleotide polymorphisms (SNPs) at the positions 3535 (T > C) and 3869 (A > G) of the mitochondrial genome that allows for the differentiation of P. simium from P. vivax (de Alvarenga et al. Reference de Alvarenga, Culleton, de Pina-costa, Rodrigues, Bianco, Silva, Nunes, de Souza, Hirano, Moreira, Pissinatti, de Abreu, Lisboa Areas, Lourenço-de-Oliveira, Zalis, Ferreira-da-Cruz, Brasil, Daniel-Ribeiro and de Brito2018).
Publicly available P. vivax sequence data
We analysed publicly available high-quality WGS data from 1050 P. vivax isolates originating from 4 continents. To this end, we first downloaded sequences from 779 P. vivax isolates from the MalariaGEN Community Pv4 dataset (Adam et al. Reference Adam, Alam, Alemu, Amaratunga, Amato, Andrianaranjaka, Anstey, Aseffa, Ashley, Assefa, Auburn, Barber, Barry, Batista Pereira, Cao, Chau, Chotivanich, Chu, Dondorp, Drury, Echeverry, Erko, Espino, Fairhurst, Faiz, Fernanda Villegas, Gao, Golassa, Goncalves, Grigg, Hamedi, Hien, Htut, Johnson, Karunaweera, Khan, Krudsood, Kwiatkowski, Lacerda, Ley, Lim, Liu, Llanos-Cuentas, Lon, Lopera-Mesa, Marfurt, Michon, Miotto, Mohammed, Mueller, Namaik-Larp, Newton, Nguyen, Nosten, Noviyanti, Pava, Pearson, Petros, Phyo, Price, Pukrittayakamee, Rahim, Randrianarivelojosia, Rayner, Rumaseb, Siegel, Simpson, Thriemer, Tobon-Castano, Trimarsanto, Urbano Ferreira, Vélez, Wangchuk, Wellems, White, William, Yasnot and Yilma2022 and references therein; Supplementary Table 2). Next, sequence data from 271 P. vivax isolates from the Americas were added to the dataset: Brazil (Hupalo et al. Reference Hupalo, Luo, Melnikov, Sutton, Rogov, Escalante, Vallejo, Herrera, Arévalo-Herrera, Fan, Wang, Cui, Lucas, Durand, Sanchez, Baldeviano, Lescano, Laman, Barnadas, Barry, Mueller, Kazura, Eapen, Kanagaraj, Valecha, Ferreira, Roobsoong, Nguitragool, Sattabonkot, Gamboa, Kosek, Vinetz, González-Cerón, Birren, Neafsey and Carlton2016; de Oliveira et al. Reference de Oliveira, Rodrigues, Menezes, Gonçalves-Lopes, Bastos, Lima, Barbosa, Gerber, Loss de Morais, Berná, Phelan, Robello, de Vasconcelos, Alves and Ferreira2017, Reference de Azevedo, Lorenz, Chiaravalloti-Neto and Sallum2020; Benavente et al. Reference Benavente, Manko, Phelan, Campos, Nolder, Fernandez, Velez-Tobon, Castaño, Dombrowski, Marinho, Aguiar, Pereira, Sriprawat, Nosten, Moon, Sutherland, Campino and Clark2021; Mourier et al. Reference Mourier, de Alvarenga, Kaushik, de Pina-costa, Douvropoulou, Guan, Guzmán-Vega, Forrester, de Abreu, Júnior, de Souza Junior, Moreira, Hirano, Pissinatti, Ferreira-da-Cruz, de Oliveira, Arold, Jeffares, Brasil, de Brito, Culleton, Daniel-Ribeiro and Pain2021; De Meulenaere et al. Reference De Meulenaere, Cuypers, Gamboa, Laukens and Rosanas-Urgell2023; Ibrahim et al. Reference Ibrahim, Manko, Dombrowski, Campos, Benavente, Nolder, Sutherland, Nosten, Fernandez, Vélez-Tobón, Castaño, Aguiar, Pereira, da Silva Santos, Suarez-Mutis, Di Santi, Baptista, Machado, Marinho, Clark and Campino2023; Kattenberg et al. Reference Kattenberg, Monsieurs, De Meyer, De Meulenaere, Sauve, de Oliveira, Ferreira, Gamboa and Rosanas-Urgell2024), Panamá (Buyon et al. Reference Buyon, Santamaria, Early, Quijada, Barahona, Lasso, Avila, Volkman, Marti, Neafsey and Obaldia2020) and Peru (Kattenberg et al. Reference Kattenberg, Nguyen, Nguyen, Sauve, Nguyen, Chopo-Pizarro, Trimarsanto, Monsieurs, Guetens, Nguyen, Esbroeck, Auburn, Nguyen and Rosanas-Urgell2022, Reference Kattenberg, Monsieurs, De Meyer, De Meulenaere, Sauve, de Oliveira, Ferreira, Gamboa and Rosanas-Urgell2024; De Meulenaere et al. Reference De Meulenaere, Cuypers, Gamboa, Laukens and Rosanas-Urgell2023) (Supplementary Table 3, Supplementary Materials). WGS data from 3 European P. vivax strains, the isolate Ebro1944 from Spain (van Dorp et al. Reference van Dorp, Gelabert, Rieux, de Manuel, De-dios, Gopalakrishnan, Carøe, Sandoval-Velasco, Fregel, Olalde, Escosa, Aranda, Huijben, Mueller, Marquès-Bonet, Balloux, Gilbert and Lalueza-Fox2020) and the isolates STR105 and STR185 from medieval/early modern Belgium (Michel et al. Reference Michel, Skourtanioti, Pierini, Guevara, Mötsch, Kocher, Barquera, Bianco, Carlhoff, Coppola Bove, Freilich, Giffin, Hermes, Hiß, Knolle, Nelson, Neumann, Papac, Penske, Rohrlach, Salem, Semerau, Villalba-Mouco, Abadie, Aldenderfer, Beckett, Brown, Campus, Chenghwa, Cruz Berrocal, Damašek, Duffett Carlson, Durand, Ernée, Fântăneanu, Frenzel, García Atiénzar, Guillén, Hsieh, Karwowski, Kelvin, Kelvin, Khokhlov, Kinaston, Korolev, Krettek, Küßner, Lai, Look, Majander, Mandl, Mazzarello, McCormick, de Miguel Ibáñez, Murphy, Németh, Nordqvist, Novotny, Obenaus, Olmo-Enciso, Onkamo, Orschiedt, Patrushev, Peltola, Romero, Rubino, Sajantila, Salazar-García, Serrano, Shaydullaev, Sias, Šlaus, Stančo, Swanston, Teschler-Nicola, Valentin, van de Vijver, Varney, Vigil-Escalera Guirado, Waters, Weiss-Krejci, Winter, Lamnidis, Prüfer, Nägele, Spyrou, Schiffels, Stockhammer, Haak, Posth, Warinner, Bos, Herbig and Krause2024), were also used for specific analyses.
Plasmodium vivax WGS data generation
We additionally generated WGS data from 36 new clinical isolates of P. vivax (Supplementary Table 3). One isolate was obtained in 2004 from a patient with imported P. vivax malaria diagnosed in the USA after a travel to India (Rodrigues et al. Reference Rodrigues, Alves, Santamaria, Calzada, Xayavong, Parise, da Silva and Ferreira2014). DNA was isolated from 200 μL of unprocessed whole blood, using the QIAamp blood DNA kit (Qiagen, Hilden, Germany). Strand displacement amplification technology, with the primer set pvset1 (Cowell et al. Reference Cowell, Loy, Sundararaman, Valdivia, Fisch, Lescano, Baldeviano, Durand, Gerbasi, Sutherland, Nolder, Vinetz, Hahn and Winzeler2017) and phi29 DNA polymerase (New England Biolabs, Ipswich, MA, USA), was used to enrich the sample for target parasite DNA prior to sequencing. The remaining 35 isolates were derived from patients with microscopy-confirmed P. vivax infection presenting between 2018 and 2019 in malaria clinics in the town of Mâncio Lima (Rodrigues et al. Reference Rodrigues, Johansen, Ladeia, Esquivel, Corder, Tonini, Calil, Fernandes, Fontoura, Cavasini, Vinetz, Castro and Ferreira2024) and the periurban settlement of Vila Assis Brasil (Fontoura et al. Reference Fontoura, Macedo, Calil, Corder, Rodrigues, Tonini, Esquivel, Ladeia, Fernandes, Johansen, Silva, Fernandes, Ladeia-Andrade, Castro and Ferreira2024), approximately 15 km apart, in the Juruá Valley region of Acre State. This region, next to the border with Peru, was the main malaria transmission hotspot of Brazil in the early 2000s (Ferreira and Castro, Reference Ferreira and Castro2016). QIAamp DNA investigator kits (Qiagen, Hilden, Germany) were used to isolate template DNA from 50 mL of venous blood that had previously been leukocyte-depleted, as described (de Oliveira et al. Reference de Oliveira, Rodrigues, Menezes, Gonçalves-Lopes, Bastos, Lima, Barbosa, Gerber, Loss de Morais, Berná, Phelan, Robello, de Vasconcelos, Alves and Ferreira2017). The presence of a single-species infection was confirmed by a species-specific TaqMan assay as described (Rodrigues et al. Reference Rodrigues, Johansen, Ladeia, Esquivel, Corder, Tonini, Calil, Fernandes, Fontoura, Cavasini, Vinetz, Castro and Ferreira2024). No genome amplification step prior to sequencing was applied to these samples. Illumina UDI libraries (Illumina, San Diego, CA) were prepared to generate paired-end 150 base pair-long sequence reads, on an Illumina NovaSeq 6000 platform, at the Institute for Genome Sciences, University of Maryland School of Medicine, Baltimore, USA. Raw fastq files were filtered for quality and those with mean quality scores ≤ 30 (expected base call accuracy ≤ 99.9%) were excluded. Reads of new WGS data described in this article have been deposited into the NCBI Sequence Read Archive (accession numbers in Supplementary Table 3).
WGS data processing and variant calling
Newly generated sequence data (from 36 P. vivax isolates) and those downloaded from SRA (from 1051 P. vivax and 31 P. simium isolates) were processed similarly. Quality-filtered fastq files were mapped onto the 20.8-megabase (Mb) core PvP01 genome (Auburn et al. Reference Auburn, Böhme, Steinbiss, Trimarsanto, Hostetler, Sanders, Gao, Nosten, Newbold, Berriman, Price and Otto2016), defined as in Table S2 of Daron et al. (Reference Daron, Boissière, Boundenga, Ngoubangoye, Houze, Arnathau, Sidobre, Trape, Durand, Renaud, Fontaine, Prugnolle and Rougeron2021), with the Burrows–Wheeler aligner (Li, Reference Li2013). Use of the PvP01 reference genome allows for comparisons with a large body of published genome data (Adam et al. Reference Adam, Alam, Alemu, Amaratunga, Amato, Andrianaranjaka, Anstey, Aseffa, Ashley, Assefa, Auburn, Barber, Barry, Batista Pereira, Cao, Chau, Chotivanich, Chu, Dondorp, Drury, Echeverry, Erko, Espino, Fairhurst, Faiz, Fernanda Villegas, Gao, Golassa, Goncalves, Grigg, Hamedi, Hien, Htut, Johnson, Karunaweera, Khan, Krudsood, Kwiatkowski, Lacerda, Ley, Lim, Liu, Llanos-Cuentas, Lon, Lopera-Mesa, Marfurt, Michon, Miotto, Mohammed, Mueller, Namaik-Larp, Newton, Nguyen, Nosten, Noviyanti, Pava, Pearson, Petros, Phyo, Price, Pukrittayakamee, Rahim, Randrianarivelojosia, Rayner, Rumaseb, Siegel, Simpson, Thriemer, Tobon-Castano, Trimarsanto, Urbano Ferreira, Vélez, Wangchuk, Wellems, White, William, Yasnot and Yilma2022). Alignments with an average genome sequence depth < 5 × were excluded. The remaining alignments were merged into Binary Alignment/Map (BAM) files and processed following the Genome Analysis Toolkit (GATK version 4.4.0) best practices (https://software.broadinstitute.org/gatk/best-practices/; McKenna et al. Reference McKenna, Hanna, Banks, Sivachenko, Cibulskis, Kernytsky, Garimella, Altshuler, Gabriel, Daly and DePristo2010). We used the GATK tool Mark Duplicates to identify and remove duplicated reads and base recalibrator and apply BQSR tools to detect and correct errors in base quality scores. Variant calling was carried out using the HaplotypeCaller module of GATK in the GVCF mode for joint genotyping of multiple samples. Variants were removed according to the following GATK variant filtration criteria: read depth (DP) < 5, variant confidence/quality by depth (QD) < 2.0, strand bias (FS) > 60.0, root mean square of the mapping quality (MQ) < 40.0, mapping quality rank sum (MQRankSum) < −12.5, read position rank sum (ReadPosRankSum) < −8.0, quality (QUAL) < 30.0. Non-biallelic SNPs and those with a minor allele frequency <0.001 (singletons) were also removed.
Complexity of infection and genetic diversity
We used the within-host diversity statistic F WS, calculated with the R package moimix (https://github.com/bahlolab/moimix), to distinguish between single- and multiple-clone infections. Multiple-clone infections, defined by F WS < 0.95 (Manske et al. Reference Manske, Miotto, Campino, Auburn, Almagro-Garcia, Maslen, O’Brien, Djimde, Doumbo, Zongo, Ouedraogo, Michon, Mueller, Siba, Nzila, Borrmann, Kiara, Marsh, Jiang, Su, Amaratunga, Fairhurst, Socheat, Nosten, Imwong, White, Sanders, Anastasi, Alcock, Drury, Oyola, Quail, Turner, Ruano-Rubio, Jyothi, Amenga-Etego, Hubbart, Jeffreys, Rowlands, Sutherland, Roper, Mangano, Modiano, Tan, Ferdig, Amambua-Ngwa, Conway, Takala-Harrison, Plowe, Rayner, Rockett, Clark, Newbold, Berriman, MacInnis and Kwiatkowski2012), were deconvoluted with DEploid (https://github.com/DEploid-dev/DEploid). Haploid genome data from single-clone infections and phased data from dominant clones within multiple-clone infections were used in all downstream analyses. VCFtools (Danecek et al. Reference Danecek, Auton, Abecasis, Albers, Banks, DePristo, Handsaker, Lunter, Marth, Sherry, McVean and Durbin2011) was used to estimate π, defined as the average number of pairwise differences per site between pairs of DNA sequences, as a measure of genome-wide nucleotide diversity in P. simium and P. vivax populations. The ratio of π values in P. simium to P. vivax was calculated for 3615 orthologous gene pairs to identify genes with highest or lowest diversity in P. simium relative to P. vivax.
Clustering and differentiation among P. simium and P. vivax genomes
Principal component analysis (PCA), as implemented in PLINK version 1.90 (Purcell et al. Reference Purcell, Neale, Todd-Brown, Thomas, Ferreira, Bender, Maller, Sklar, de Bakker, Daly and Sham2007), was used to determine clustering patterns among P. simium and P. vivax isolates based on genome-wide genotypes. Missing genotypes were imputed using Beagle version 5.4 (Browning et al. Reference Browning, Zhou and Browning2018). Variant sites in linkage disequilibrium (LD) (r 2 value > 0.2) within 50-kb windows (step size of 10 base pairs, bp) were pruned to obtain a set of unlinked SNPs. Additional analyses run with the smartPCA software (Patterson et al. Reference Patterson, Price and Reich2006) included low-coverage genomes from European P. vivax isolates. We used eigenvectors computed for high-quality WGS data and projected low-coverage data from European isolates onto the axes of variation, as described elsewhere (Michel et al. Reference Michel, Skourtanioti, Pierini, Guevara, Mötsch, Kocher, Barquera, Bianco, Carlhoff, Coppola Bove, Freilich, Giffin, Hermes, Hiß, Knolle, Nelson, Neumann, Papac, Penske, Rohrlach, Salem, Semerau, Villalba-Mouco, Abadie, Aldenderfer, Beckett, Brown, Campus, Chenghwa, Cruz Berrocal, Damašek, Duffett Carlson, Durand, Ernée, Fântăneanu, Frenzel, García Atiénzar, Guillén, Hsieh, Karwowski, Kelvin, Kelvin, Khokhlov, Kinaston, Korolev, Krettek, Küßner, Lai, Look, Majander, Mandl, Mazzarello, McCormick, de Miguel Ibáñez, Murphy, Németh, Nordqvist, Novotny, Obenaus, Olmo-Enciso, Onkamo, Orschiedt, Patrushev, Peltola, Romero, Rubino, Sajantila, Salazar-García, Serrano, Shaydullaev, Sias, Šlaus, Stančo, Swanston, Teschler-Nicola, Valentin, van de Vijver, Varney, Vigil-Escalera Guirado, Waters, Weiss-Krejci, Winter, Lamnidis, Prüfer, Nägele, Spyrou, Schiffels, Stockhammer, Haak, Posth, Warinner, Bos, Herbig and Krause2024).
ADMIXTURE version 1.3.0 (Alexander and Lange, Reference Alexander and Lange2011) was used to investigate the genetic ancestry of P. simium and Latin American P. vivax isolates. An unsupervised ADMIXTURE analysis was run to assign isolates to K putative ancestral populations according to SNP frequencies. The LD-pruned data set with imputed missing genotypes was used in the analysis. We calculated cross-validation error rates to determine the most likely number of ancestral populations, with K between 1 and 15 (Alexander et al. Reference Alexander, Novembre and Lange2009). Wright’s Fixation Index, F ST (Weir and Cockerham implementation in VCFtools), was used to determine differentiation between P. simium and Latin American populations of P. vivax. To further explore the genetic relationships between P. simium and present-day P. vivax populations from Latin America, we used the R package admixr (Petr, Reference Petr2024) to compute the f4 statistic, f4 (W, X; Y, Z) (Patterson et al. Reference Patterson, Moorjani, Luo, Mallick, Rohland, Zhan, Genschoreck, Webster and Reich2012). We compared the proportion of derived alleles shared between P. simium (Y) and pairs of P. vivax populations (W and X), using P. cynomolgi strain B (RefSeq assembly: GCF_000321355.1; Tachibana et al. Reference Tachibana, Sullivan, Kawai, Nakamura, Kim, Goto, Arisue, Palacpac, Honma, Yagi, Tougan, Katakai, Kaneko, Mita, Kita, Yasutomi, Sutton, Shakhbatyan, Horii, Yasunaga, Barnwell, Escalante, Carlton and Tanabe2012) as the outgroup (Z). The P. vivax population from Mexico (X), the genetically most distant population from P. simium samples in PCA and F ST analysis, was used in all pairwise comparisons.
Genetic relatedness and positive selection among P. simium isolates
We used identity-by-descent (IBD) analysis as implemented in hmmIBD (Schaffner et al. Reference Schaffner, Taylor, Wong, Wirth and Neafsey2018) to measure the genetic relatedness between pairs of P. simium samples from monkeys and from humans, and between P. simium isolates from different states (São Paulo, Rio de Janeiro and Espírito Santo). hmmIBD implements a hidden Markov model-based approach that accounts for recombination to systematically search for genomic segments that are inferred to have descended from a common ancestor without intervening recombination and to estimate the proportion of shared ancestry between genomes. A recombination rate of 13.5 kilobases (kb) per centiMorgan was assumed (Miles et al. Reference Miles, Iqbal, Vauterin, Pearson, Campino, Theron, Gould, Mead, Drury, O’Brien, Ruano Rubio, MacInnis, Mwangi, Samarakoon, Ranford-Cartwright, Ferdig, Hayton, Su, Wellems, Rayner, McVean and Kwiatkowski2016) and the -n option was set to 25, to call IBD segments from common ancestors within the last 25 generations (approximately 12.5 years, assuming 2 generations per year). Genetic relatedness networks were drawn to connect sample pairs with at least 25% (equivalent to half-siblings), 50% (equivalent to meiotic siblings) or 90% of genomes IBD. We estimated the distribution of IBD coverage (i.e., the number of IBD segments overlapping each position) and identified IBD peaks (i.e., chromosome regions with IBD sharing > 2 standard deviations above the mean) across the P. simium genome using the ibdutils command line tool (Guo et al. Reference Guo, Borda, Laboulaye, Spring, Wojnarski, Vesely, Silva, Waters, O’Connor and Takala-Harrison2024; https://github.com/bguo068/ishare).
Demographic history of P. simium and P. vivax populations from Brazil
We used pairwise sequentially Markovian coalescent (PSMC) analysis (Li and Durbin, Reference Li and Durbin2011; Mather et al. Reference Mather, Traves and Ho2019) to infer the historical effective population size (N e) of P. simium and P. vivax populations from Latin America. We ran PSMC analysis (https://github.com/lh3/psmc) assuming a mutation rate (μ) of 1 × 10−9 and a generation time (g) of 0.18 (Daron et al. Reference Daron, Boissière, Boundenga, Ngoubangoye, Houze, Arnathau, Sidobre, Trape, Durand, Renaud, Fontaine, Prugnolle and Rougeron2021), with the 5 samples with the highest sequence coverage for each species. We performed 100 bootstrap replicates by randomly resampling the variants used in the analysis.
Results
WGS data for P. simium and P. vivax
We analysed WGS data from 31 P. simium isolates, with an average of 1.1 × 106 sequence reads per sample, for an average depth of coverage of 49.7 × (range between 3.4 × and 233.0 × among samples). We identified 55 682 high-confidence SNPs (Supplementary Table 4, Supplementary Materials). Only 3 (9.7%) samples contained 2 or more clones; those sequence data sets were deconvoluted, with the dominant genotype retained for further analysis. Although P. simium sequences had been generated from unprocessed blood (i.e. leukocytes had not been removed prior to DNA extraction), most (50.7%) reads mapped to the reference P. vivax genome. On average, 73.9% (range: 29.4–93.6%) of the PvP01 core genome was covered with a read depth ≥ 5 × .
A total of 35 new P. vivax WGS data sets from Brazil were generated in this study. We obtained an average of 10.6 × 106 reads per sample, resulting in an average read depth of 32.4 × (range: 5.5 × to 55.6 ×), and a total of 60 568 high-confidence SNPs were identified after filtering (Supplementary Table 5, Supplementary Materials). Eight (22.8%) clinical samples contained 2 or more clones (F ws ≤ 0.95) and the sequence data were deconvoluted. On average, 81.4% of the reads obtained from these leukocyte-depleted samples mapped to the reference PvP01 core genome, with an average breadth of coverage of 95.9% (range, 80.0% to 98.2%) of the reference genome mapped with a read depth ≥5 × .
Plsmodium simium originated in South America and is closest to the P. vivax population from Brazil
As in previous studies (Adam et al. Reference Adam, Alam, Alemu, Amaratunga, Amato, Andrianaranjaka, Anstey, Aseffa, Ashley, Assefa, Auburn, Barber, Barry, Batista Pereira, Cao, Chau, Chotivanich, Chu, Dondorp, Drury, Echeverry, Erko, Espino, Fairhurst, Faiz, Fernanda Villegas, Gao, Golassa, Goncalves, Grigg, Hamedi, Hien, Htut, Johnson, Karunaweera, Khan, Krudsood, Kwiatkowski, Lacerda, Ley, Lim, Liu, Llanos-Cuentas, Lon, Lopera-Mesa, Marfurt, Michon, Miotto, Mohammed, Mueller, Namaik-Larp, Newton, Nguyen, Nosten, Noviyanti, Pava, Pearson, Petros, Phyo, Price, Pukrittayakamee, Rahim, Randrianarivelojosia, Rayner, Rumaseb, Siegel, Simpson, Thriemer, Tobon-Castano, Trimarsanto, Urbano Ferreira, Vélez, Wangchuk, Wellems, White, William, Yasnot and Yilma2022; Kattenberg et al. Reference Kattenberg, Monsieurs, De Meyer, De Meulenaere, Sauve, de Oliveira, Ferreira, Gamboa and Rosanas-Urgell2024; Michel et al. Reference Michel, Skourtanioti, Pierini, Guevara, Mötsch, Kocher, Barquera, Bianco, Carlhoff, Coppola Bove, Freilich, Giffin, Hermes, Hiß, Knolle, Nelson, Neumann, Papac, Penske, Rohrlach, Salem, Semerau, Villalba-Mouco, Abadie, Aldenderfer, Beckett, Brown, Campus, Chenghwa, Cruz Berrocal, Damašek, Duffett Carlson, Durand, Ernée, Fântăneanu, Frenzel, García Atiénzar, Guillén, Hsieh, Karwowski, Kelvin, Kelvin, Khokhlov, Kinaston, Korolev, Krettek, Küßner, Lai, Look, Majander, Mandl, Mazzarello, McCormick, de Miguel Ibáñez, Murphy, Németh, Nordqvist, Novotny, Obenaus, Olmo-Enciso, Onkamo, Orschiedt, Patrushev, Peltola, Romero, Rubino, Sajantila, Salazar-García, Serrano, Shaydullaev, Sias, Šlaus, Stančo, Swanston, Teschler-Nicola, Valentin, van de Vijver, Varney, Vigil-Escalera Guirado, Waters, Weiss-Krejci, Winter, Lamnidis, Prüfer, Nägele, Spyrou, Schiffels, Stockhammer, Haak, Posth, Warinner, Bos, Herbig and Krause2024), PCA revealed the distinct geographic structure of P. vivax populations from around the globe. The first 2 principal components (PCs) captured one-third of the overall genetic variation and defined 3 main clusters: 1 with samples from Latin America, clearly separated from samples from Africa, South and West Asia and from a third cluster with samples from East and Southeast Asia and Oceania (Figure 1). Importantly, all P. simium samples clustered close to Latin American P. vivax populations in the PCA space. Historical P. vivax samples from Europe that were suitable for population genetic analysis – Ebro1944 from Spain (van Dorp et al. Reference van Dorp, Gelabert, Rieux, de Manuel, De-dios, Gopalakrishnan, Carøe, Sandoval-Velasco, Fregel, Olalde, Escosa, Aranda, Huijben, Mueller, Marquès-Bonet, Balloux, Gilbert and Lalueza-Fox2020) and STR105 and STR185 from Belgium (Michel et al. Reference Michel, Skourtanioti, Pierini, Guevara, Mötsch, Kocher, Barquera, Bianco, Carlhoff, Coppola Bove, Freilich, Giffin, Hermes, Hiß, Knolle, Nelson, Neumann, Papac, Penske, Rohrlach, Salem, Semerau, Villalba-Mouco, Abadie, Aldenderfer, Beckett, Brown, Campus, Chenghwa, Cruz Berrocal, Damašek, Duffett Carlson, Durand, Ernée, Fântăneanu, Frenzel, García Atiénzar, Guillén, Hsieh, Karwowski, Kelvin, Kelvin, Khokhlov, Kinaston, Korolev, Krettek, Küßner, Lai, Look, Majander, Mandl, Mazzarello, McCormick, de Miguel Ibáñez, Murphy, Németh, Nordqvist, Novotny, Obenaus, Olmo-Enciso, Onkamo, Orschiedt, Patrushev, Peltola, Romero, Rubino, Sajantila, Salazar-García, Serrano, Shaydullaev, Sias, Šlaus, Stančo, Swanston, Teschler-Nicola, Valentin, van de Vijver, Varney, Vigil-Escalera Guirado, Waters, Weiss-Krejci, Winter, Lamnidis, Prüfer, Nägele, Spyrou, Schiffels, Stockhammer, Haak, Posth, Warinner, Bos, Herbig and Krause2024) – clustered with Latin American P. vivax, to which they are more similar than to P. simium (Supplementary Figure 2). These findings are consistent with the origin of P. simium in Latin America.
Figure 1. Global P. vivax and P. simium population structure revealed by standard PCA. Data analysed corresponds to linkage disequilibrium-pruned biallelic SNPs. We show the first 2 PCs, which together account for 35.9% of the overall variance. Each symbol – circles for P. simium and triangles for P. vivax – represents a single isolate and was coloured according to the country of origin of the sample.
A regional PCA was done to further investigate the origin of P. simium. All P. simium isolates clustered together in the regional PCA, regardless of their state of origin in Brazil (Figure 2). The P. simium cluster includes 12 human-derived samples from São Paulo that were originally labelled as P. vivax (Ibrahim et al. Reference Ibrahim, Manko, Dombrowski, Campos, Benavente, Nolder, Sutherland, Nosten, Fernandez, Vélez-Tobón, Castaño, Aguiar, Pereira, da Silva Santos, Suarez-Mutis, Di Santi, Baptista, Machado, Marinho, Clark and Campino2023). PCA also revealed a closer affinity of P. simium to lineages of P. vivax from the Amazon Basin of Brazil (states of Acre, Amapá, Amazonas, Pará and Rondônia), compared to P. vivax populations from Peru and Colombia and those from Central America (Panama, Nicaragua and El Salvador) and Mexico. Moreover, P. vivax isolates from Europe clustered together with the P. vivax populations from Brazil and all other Latin American countries, but not with P. simium (Supplementary Figure 3), reflecting the extent of divergence between present-day P. simium and the founding European lineages of P. vivax introduced in Brazil.
Figure 2. Plasmodium simium and P. vivax population structure in Latin America revealed by PCA). Analysis included a total of 495 isolates (P. simium: n = 31; P. vivax: n = 464). We display the first 3 PCs, which together account for 45.9% of the overall variance. Each symbol – circles for P. simium and squares for P. vivax – represents a single isolate and was coloured according to the country or state (within Brazil) of origin of the sample. Locations of each state in Brazil are shown in Supplementary Figure 1.
Unsupervised ADMIXTURE analysis was used to examine shared ancestry patterns of P. simium and P. vivax at the regional and country level. Cross-validation error decreased with increasing K until K = 15 (Supplementary Figure 4). We chose to display results of the analyses with the smallest K value that captures most of the geographic structure in the data (K = 4), the K-value that maximizes clustering on a country level (K = 10) and the K-value associated with the lowest cross-validation error (K = 15) (Figure 3). At K = 4, 3 regional P. vivax clusters can be seen: Brazil (apple green); Peru, Panama and Nicaragua (light sea green); and Colombia, El Salvador and Mexico (purple). At K = 10, most P. vivax isolates were assigned to country-specific populations but some samples from Brazil, Panama and Peru appear to be admixed. Indeed, samples from Brazil appear to have ancestry in 1 of 2 different ancestral populations (blue and green), or both (admixed samples), while some samples from Peru appear to share ancestry with the green population from Brazil. At K = 15, the P. simium population appeared to comprise 2 distinct subpopulations, 1 from São Paulo and the other comprising isolates from Rio de Janeiro and Espírito Santo. Importantly, at all K values, P. simium isolates were assigned to separate population(s), with negligible admixture with regional or country-specific P. vivax populations (Figure 3). Moreover, the elevated estimates of population differentiation in pairwise comparisons of P. simium with P. vivax populations from Mexico (F ST = 0.29), Panama (F ST = 0.22), Colombia (F ST = 0.20), Peru (F ST = 0.20) and Brazil (F ST = 0.16) indicate limited historical gene flow between species (Supplementary Figure 5, Supplementary Materials).
Figure 3. Unsupervised ADMIXTURE analysis of P. simium and P. vivax from Latin America. Three layers correspond to k = 4, k = 10 and k = 15 populations. In ADMIXTURE bar plots, each isolate is represented by a bar that is coloured to indicate the proportion of the genome (from 0 to 1), with ancestry from each of k putative ancestral source populations. Admixed samples (those with ancestry from more than one source population) are represented by bar segments of different colours.
Previous analyses had suggested that P. simium was most closely related to present-day P. vivax samples from Mexico, compared to other locations in Latin America (Mourier et al. Reference Mourier, de Alvarenga, Kaushik, de Pina-costa, Douvropoulou, Guan, Guzmán-Vega, Forrester, de Abreu, Júnior, de Souza Junior, Moreira, Hirano, Pissinatti, Ferreira-da-Cruz, de Oliveira, Arold, Jeffares, Brasil, de Brito, Culleton, Daniel-Ribeiro and Pain2021; de Oliveira et al. Reference de Oliveira, Rodrigues, Early, Duarte, Buery, Bueno, Catão-Dias, Cerutti, Rona, Neafsey and Ferreira2021b). However, f4 statistics with the expanded data set comprising 31 P. simium samples from 3 states in southwest Brazil indicated that this species shares significantly more derived alleles with P. vivax samples from the Amazon Basin of Brazil, followed by Peru, Colombia and Panama, compared to P. vivax samples from Mexico (Supplementary Figure 6A). These findings, consistent with our regional PCA (Figure 2) and F ST results (Supplementary Figure 5), point to Brazil as the most likely birthplace of the novel parasite lineage that resulted from parasite jumps from humans to platyrrhine monkeys (de Oliveira et al. Reference de Oliveira, Rodrigues, Early, Duarte, Buery, Bueno, Catão-Dias, Cerutti, Rona, Neafsey and Ferreira2021b; Mourier et al. Reference Mourier, de Alvarenga, Kaushik, de Pina-costa, Douvropoulou, Guan, Guzmán-Vega, Forrester, de Abreu, Júnior, de Souza Junior, Moreira, Hirano, Pissinatti, Ferreira-da-Cruz, de Oliveira, Arold, Jeffares, Brasil, de Brito, Culleton, Daniel-Ribeiro and Pain2021).
Plasmodium simium has a decreasing effective population size, smaller than P. vivax in Latin America
PSMC analysis revealed distinct historical trends for N e estimates of P. simium and P. vivax populations (Figure 4). First, looking backwards in time, P. vivax has a much deeper genealogy than P. simium, consistent with the hypothesis that P. simium has a much more recent origin. Second, P. vivax shows a substantial decrease over time, followed by a recent sharp increase in N e, suggestive of a sharp population expansion following a bottleneck. In contrast, the P. simium population shows a stable decrease in N e over the entire period (Figure 4; see Supplementary Figure 7, for bootstrap replicates), consistent with a fairly small effective population size after the split from its common ancestor with P. vivax.
Figure 4. History of effective population size, N e, for P. simium and P. vivax. Pairwise sequentially Markovian coalescent (PSMC) model estimation of historical population size changes in P. simium and P. vivax populations from Brazil. We used a mutation rate (μ) of 1 × 10−9 and generation time (g) of 0.18. Results are shown for the 5 samples with the highest sequence coverage for each species.
We next compared the overall nucleotide diversity of present-day P. simium (n = 31) with that of Latin American populations of P. vivax (n = 465). In line with previous studies (Mourier et al. Reference Mourier, de Alvarenga, Kaushik, de Pina-costa, Douvropoulou, Guan, Guzmán-Vega, Forrester, de Abreu, Júnior, de Souza Junior, Moreira, Hirano, Pissinatti, Ferreira-da-Cruz, de Oliveira, Arold, Jeffares, Brasil, de Brito, Culleton, Daniel-Ribeiro and Pain2021; de Oliveira et al. Reference de Oliveira, Rodrigues, Early, Duarte, Buery, Bueno, Catão-Dias, Cerutti, Rona, Neafsey and Ferreira2021b), we found a low diversity in the core genome of P. simium (mean π = 2.5 × 10−4), approximately 2 times lower than that of P. vivax populations from Latin America (mean π = 5.0 × 10−4, across loci) and 1.5–2.4 times lower than the diversity calculated for P. vivax populations from Brazil (mean π = 5.0 × 10−4), Colombia (mean π = 5.9 × 10−4), Mexico (mean π = 5.8 0 × 10−4) and Peru (mean π = 3.8 × 10−4). The average π across loci among isolates from Panama (π = 2.4 × 10−4) was similar to that for P. simium. Interestingly, nucleotide diversity estimates were nearly identical for P. simium samples from humans (mean π = 2.94 × 104; n = 25 samples) and platyrrhine monkeys (mean π = 2.90 × 104; n = 6 samples), as would be expected from repeat sampling from the same panmictic population.
Compared to a null distribution centred on 0, expected from identical nucleotide diversity per locus for both species, the density distribution of log ratios of π estimates is shifted to the left (Figure 5). This indicates less diversity in P. simium genes compared to their orthologs in P. vivax, with an average log ratio of −0.24 (standard deviation = 0.5). Interestingly, the right-hand tail of the distribution in Figure 5, which comprises genes with higher diversity in P. simium, includes loci that encode proteins involved in host–parasite interactions in P. vivax (Supplementary Table 6, Supplementary Materials). Among them are those encoding: the vacuolar protein sorting-associated protein 46 (PVP01_0704600), which may be linked to the production of extracellular vesicles in parasitized red blood cells, with a role in intercellular communication (Toda et al. Reference Toda, Diaz-Varela, Segui-Barber, Roobsoong, Baro, Garcia-Silva, Galiano, Gualdrón-López, Almeida, Brito, de Melo, Aparici-Herraiz, Castro-Cavadía, Monteiro, Borràs, Sabidó, Almeida, Chojnacki, Martinez-Picado, Calvo, Armengol, Carmona-Fonseca, Yasnot, Lauzurica, Marcilla, Peinado, Galinski, Lacerda, Sattabongkot, Fernandez-Becerra and Del Portillo2020; Avalos-Padilla et al. Reference Avalos-Padilla, Georgiev, Lantero, Pujals, Verhoef, Borgheti-Cardoso, Albertazzi, Dimova and Fernàndez-Busquets2021); the ookinete maturation gene OMG1 (PVP01_0609300), crucial for the invasive stage in the mosquito midgut in P. berghei (Nishi et al. Reference Nishi, Kaneko, Iwanaga and Yuda2022); the oocyst capsule protein Cap380, essential for oocyst development, sporozoite differentiation and malaria transmission in P. berghei (Srinivasan et al. Reference Srinivasan, Fujioka and Jacobs-Lorena2008; Nakayama et al. Reference Nakayama, Kimura, Kitahara, Soga, Haraguchi, Hakozaki, Sugiyama, Kusakisako, Fukumoto and Ikadai2021); the liver merozoite formation protein (PVP01_1146600), vital for the maturation of liver merozoites in P. berghei (Haussig et al. Reference Haussig, Matuschewski and Kooij2011), which in P. falciparum appears to be essential for sporozoite formation within the oocyst (Franke-Fayard et al. Reference Franke-Fayard, Marin-Mogollon, Geurten, Chevalley-Maurel, Ramesar, Kroeze, Baalbergen, Wessels, Baron, Soulard, Martinson, Aleshnick, Huijs, Subudhi, Miyazaki, Othman, Kolli, Lamers, Roques, Stanway, Murphy, Foquet, Moita, Mendes, Prudêncio, Dechering, Heussler, Pain, Wilder, Roestenberg and Janse2022); and the AP2 domain transcription factor AP2-G5 (PVP01_0940100), crucial for gametocyte maturation in P. falciparum (Shang et al. Reference Shang, Shen, Tang, He, Zhao, Wang, He, Guo, Liu, Wang, Zhu, Yang, Jiang, Zhang, Yu, Han, Culleton, Jiang, Cao, Gu and Zhang2021).
Figure 5. Nucleotide diversity in P. simium and P. vivax. Empirical density distribution of the log ratios of nucleotide diversity (π) estimates for P. simium and P. vivax ortholog genes (bars). Red line shows the null distribution centred on 0 as expected under identical nucleotide diversity between species. Data correspond to a total of 3615 orthologous gene pairs between P. vivax and P. simium.
Contemporary P. simium infects both monkeys and humans
Plasmodium simium samples from humans from São Paulo, which were originally labelled as P. vivax (Ibrahim et al. Reference Ibrahim, Manko, Dombrowski, Campos, Benavente, Nolder, Sutherland, Nosten, Fernandez, Vélez-Tobón, Castaño, Aguiar, Pereira, da Silva Santos, Suarez-Mutis, Di Santi, Baptista, Machado, Marinho, Clark and Campino2023), were found to share more derived alleles with P. simium from humans and from monkeys from Espírito Santo and Rio de Janeiro than with any of the P. vivax populations from Latin America (Supplementary Figure 6B). These findings further confirm that the samples of P. vivax-related parasites infecting humans from São Paulo belong to the P. simium clade.
Parasite relatedness networks revealed some recent gene flow between P. simium populations across states. Although no sample pair from different states displayed ≥50% of the genome IBD (Supplementary Figure 8), examples of ≥25% IBD between pairs of isolates were found in a single cluster comprising 6 samples from São Paulo and 3 from Rio de Janeiro (Figure 6). We found only a very low proportion of IBD sharing between pairs of P. simium isolates originating from sites >450 km apart (Supplementary Figure 9). However, we found a pair of very closely related isolates, with ≥90% of the genome IBD, circulating within the state of São Paulo (Supplementary Figure 10), while there were no clonal lineages of P. simium in which isolates shared ≥99% of the genome IBD.
Figure 6. Relatedness network of Plasmodium simium samples. Samples were collected from humans and platyrrhine monkeys from the states of São Paulo, Rio de Janeiro and Espírito Santo, Southeastern Brazil (n = 31). Nodes represent individual samples that are coloured according to the state of origin; edges connect samples with mean pairwise ancestry sharing ≥ 0.25 (equivalent to half-siblings). Unconnected nodes indicate isolates that do not share at least 25% of their genome-wide ancestry with other isolates from the same or different states.
Notably, we found 2 clusters of parasites sharing ≥25% of the genome IBD that were derived from different mammalian hosts (Supplementary Figure 11). The larger cluster comprised 5 samples from monkeys (4 from São Paulo and 1 from Rio de Janeiro) and 4 samples from humans (2 from São Paulo and 2 from Rio de Janeiro), while the smaller cluster comprised 1 sample from a monkey and 4 from humans, all from Rio de Janeiro. These findings are consistent with recent gene flow between parasites from human and nonhuman hosts.
Positive selection in P. simium reveals signals of adaptation to new vertebrate host and mosquito vectors
We next searched for signatures of positive selection across the P. simium genome. Positive selection is expected to increase IBD sharing at the target locus and neighbouring sites, generating peaks of within-population IBD sharing (Guo et al. Reference Guo, Borda, Laboulaye, Spring, Wojnarski, Vesely, Silva, Waters, O’Connor and Takala-Harrison2024). A genome-wide scan identified 7 validated IBD peaks, possibly associated with selective sweeps during the parasite′s adaptation to new hosts (Figure 7). IBD peaks mapped to chromosomes 3, 9, 12 and 14 and comprise several annotated genes of potential interest (Supplementary Table 7, Supplementary Materials), including some transcription factors encoding an AP2 domain (DNA-binding) that may be involved in red blood cell invasion, gametocytogenesis, oocyst formation and sporozoite formation (reviewed by Singhal et al. Reference Singhal, Prata, Bonnell and Llinás2024). For example, the chromosome 3 peak comprises genes that encode the AP2 domain transcription factor AP2-SP2 (PVP01_0303400), which is crucial for oocyst maturation in the vector (Modrzynska et al. Reference Modrzynska, Pfander, Chappell, Yu, Suarez, Dundas, Gomes, Goulding, Rayner, Choudhary and Billker2017), and the 6-cysteine protein P36 (PVP01_0303700), required by sporozoites for invasion and establishment of the parasitophorous vacuole within hepatocytes (Arredondo et al. Reference Arredondo, Swearingen, Martinson, Steel, Dankwa, Harupa, Camargo, Betz, Vigdorovich, Oliver, Kangwanrangsan, Ishino, Sather, Mikolajczak, Vaughan, Torii, Moritz and Kappe2018), while one of the chromosome 9 peaks contains the gene for the red blood cell ligand apical membrane antigen 1 (PVP01_0934200), a major vaccine candidate antigen (Drew et al. Reference Drew, Wilson, Weiss, Yeoh, Henshall, Crabb, Dutta, Gilson and Beeson2023).
Figure 7. Genomic regions in P. simium under positive selection. Domains putatively under strong positive selection revealed by identity-by-descent (IBD) analysis of the P. simium genome. We display peaks of IBD coverage and proportion of shared ancestry along the 14 chromosomes. Red shading indicates validated peaks likely associated with selective sweeps (Guo et al. Reference Guo, Borda, Laboulaye, Spring, Wojnarski, Vesely, Silva, Waters, O’Connor and Takala-Harrison2024).
Discussion
Plasmodium simium offers a compelling example of zoonotic malaria parasite originating from recent host-shift speciation events. Two main scenarios have been suggested for the origin of P. simium: (1) human-to-monkey shifts took place after P. vivax lineages were introduced in Brazil, mostly by European settlers but perhaps also by enslaved populations displaced from Africa (de Oliveira et al. Reference de Oliveira, Rodrigues, Duarte, Rona and Ferreira2021a), and (2) P. vivax stocks from different source populations from Europe, Africa and/or Asia were introduced in the New World. While ‘generalist’ lineages were able to infect local platyrrhine monkeys in addition to humans, giving rise to P. simium, more ‘specialist’ lineages are expected to remain limited to humans (Rougeron et al. Reference Rougeron, Daron, Fontaine and Prugnolle2022). Our findings support the first scenario. First, it was found low within-species nucleotide diversity among present-day lineages of P. simium circulating in both humans and monkeys, consistent with its recent origin from the same source population. Separate host-shift events may have occurred elsewhere in South America – e.g., in Colombia, where P. vivax-related parasites appear to infect wild nonhuman primates (Rondón et al. Reference Rondón, León, Link and González2019), but monkey-derived P. vivax samples outside Southeast Brazil have not been characterized at the genome level. Second, our global population structure analysis shows that all P. simium and Latin American P. vivax isolates cluster together (Figure 1) and both display genetic similarity to now-extinct European lines of P. vivax (Supplementary Figure 2, Supplementary Materials). The seemly common geographic origin of P. simium and Latin American P. vivax isolates argues against the hypothesis that separate introductions of parasites, from different regions or continents, originated monkey-adapted vs. exclusively human parasites found nowadays in Brazil. Plasmodium simium challenges current malaria control and elimination efforts. About 0.05% of the malaria cases recorded each year in Brazil are acquired outside the Amazon Basin (Ferreira and Castro, Reference Ferreira and Castro2016), especially in forest fringes, where nonhuman primates serve as parasite reservoirs for spillback events (Abreu et al. Reference Abreu, Santos, Mello, Gomes, Alvarenga, Gomes, Vargas, Bianco-Júnior, Pina-Costa, Teixeira, Romano, Manso, Pelajo-Machado, Brasil, Daniel-Ribeiro, Brito, Ferreira-da-Cruz and Lourenço-de-Oliveira2019; Duarte et al. Reference Duarte, Fernandes, Silva, Sicchi, Mucci, Curado, Fernandes, Medeiros-Sousa, Ceretti-Junior, Marrelli, Evangelista, Teixeira, Summa, Nardi, Garnica, Loss, Buery, Jr, Pacheco, Escalante, Sallum and Laporta2021) and human–vector contact is favoured by environmental change (Medeiros-Sousa et al. Reference Medeiros-Sousa, Laporta, Coutinho, Mucci and Marrelli2021). The vast majority of autochthonous human malaria cases in the Brazilian Atlantic Forest biome are reportedly caused by P. vivax (Garcia et al. Reference Garcia, Abrahão, Oliveira, Henriques, de Pina-costa, Siqueira and Ramalho2022), but available molecular evidence suggests that most, if not all, cases attributed to P. vivax in these settings are indeed caused by the zoonotic parasite P. simium (Brasil et al. Reference Brasil, Zalis, de Pina-costa, Siqueira, Júnior, Silva, Areas, Pelajo-Machado, de Alvarenga, da Silva Santelli, Albuquerque, Cravo, Santos de Abreu, Peterka, Zanini, Suárez Mutis, Pissinatti, Lourenço-de-Oliveira, de Brito, Ferreira-da-Cruz, Culleton and Daniel-Ribeiro2017; Buery et al. Reference Buery, Rodrigues, Natal, Salla, Loss, Vicente, Rezende, Duarte, Fux, Malafronte, Falqueto and Cerutti2017; de Oliveira et al. Reference de Oliveira, Rodrigues, Early, Duarte, Buery, Bueno, Catão-Dias, Cerutti, Rona, Neafsey and Ferreira2021b). Here, further support to this hypothesis is provided by showing that the genome sequences of 12 parasites originally labelled as P. vivax circulating among humans in São Paulo (Ibrahim et al. Reference Ibrahim, Manko, Dombrowski, Campos, Benavente, Nolder, Sutherland, Nosten, Fernandez, Vélez-Tobón, Castaño, Aguiar, Pereira, da Silva Santos, Suarez-Mutis, Di Santi, Baptista, Machado, Marinho, Clark and Campino2023) cluster together with P. simium lineages of human and nonhuman origin and can be clearly differentiated from P. vivax populations from the Amazon (Figure 2 and Supplementary Figure 6B). Importantly, zoonotic malaria transmission to humans occurs in the vicinity of Rio de Janeiro, São Paulo and other major cities in Southeast Brazil and can undermine the country′s elimination efforts (Fornace et al. Reference Fornace, Drakeley, Lindblade, Jelip and Ahmed2023a; Fornace et al. Reference Fornace, Laporta, Vythilingham, Chua, Ahmed, Jeyaprakasam, Duarte, Amir, Phang, Drakeley, Sallum and Lau2023b).
Previous analyses of substantially smaller P. simium sequence datasets had suggested that this parasite was most closely related to present-day P. vivax samples from Mexico, a relatively inbred population that has experienced a steady decline in recent years, compared to those from Brazil and other locations in Latin America (Mourier et al. Reference Mourier, de Alvarenga, Kaushik, de Pina-costa, Douvropoulou, Guan, Guzmán-Vega, Forrester, de Abreu, Júnior, de Souza Junior, Moreira, Hirano, Pissinatti, Ferreira-da-Cruz, de Oliveira, Arold, Jeffares, Brasil, de Brito, Culleton, Daniel-Ribeiro and Pain2021; de Oliveira et al. Reference de Oliveira, Rodrigues, Early, Duarte, Buery, Bueno, Catão-Dias, Cerutti, Rona, Neafsey and Ferreira2021b). However, these findings must not be over-interpreted, since P. vivax sequence data from Mexico originated from relatively isolated foci of residual malaria transmission in the southern part of the country (Hupalo et al. Reference Hupalo, Luo, Melnikov, Sutton, Rogov, Escalante, Vallejo, Herrera, Arévalo-Herrera, Fan, Wang, Cui, Lucas, Durand, Sanchez, Baldeviano, Lescano, Laman, Barnadas, Barry, Mueller, Kazura, Eapen, Kanagaraj, Valecha, Ferreira, Roobsoong, Nguitragool, Sattabonkot, Gamboa, Kosek, Vinetz, González-Cerón, Birren, Neafsey and Carlton2016). The present analysis combines additional P. simium sequences and a wide range of P. vivax sequence data from Latin America to show that P. simium is clearly more closely related to P. vivax from Brazil (Supplementary Figure 6A). However, evidence from PCA (Figure 2), ADMIXTURE (Figure 3) and F ST analyses (Supplementary Figure 5) shows a clear divergence between P. simium and P. vivax from Brazil and other Latin American countries. At the species level, population differentiation follows an isolation-by-distance model, with little, if any, IBD sharing between P. simium samples from sites >450 km apart (Supplementary Figure 8) and across states (Figure 6). These results are consistent with focal parasite transmission among nonhuman primates, and occasionally humans, within discontinuous forest fragments intermingled with malaria-free areas in Southeast Brazil (Ferreira and Castro, Reference Ferreira and Castro2016; Duarte et al. Reference Duarte, Fernandes, Silva, Sicchi, Mucci, Curado, Fernandes, Medeiros-Sousa, Ceretti-Junior, Marrelli, Evangelista, Teixeira, Summa, Nardi, Garnica, Loss, Buery, Jr, Pacheco, Escalante, Sallum and Laporta2021). However, the substantial IBD sharing (≥25%) found between parasites from humans and platyrrhine monkeys (Supplementary Figure 11, Supplementary Materials) suggests occasional cross-species transmission events (Su and Wu, Reference Su and Wu2021).
Mutations, insertions and deletions in key erythrocyte invasion ligands of malaria parasites are commonly seen in host-shift speciation events (de Oliveira et al. Reference de Oliveira, Rodrigues, Duarte, Rona and Ferreira2021a). Accordingly, all P. simium isolates analysed so far display a deletion of >40% of the coding sequence of the locus encoding the reticulocyte binding protein 2a, rbp2a (Mourier et al. Reference Mourier, de Alvarenga, Kaushik, de Pina-costa, Douvropoulou, Guan, Guzmán-Vega, Forrester, de Abreu, Júnior, de Souza Junior, Moreira, Hirano, Pissinatti, Ferreira-da-Cruz, de Oliveira, Arold, Jeffares, Brasil, de Brito, Culleton, Daniel-Ribeiro and Pain2021; de Oliveira et al. Reference de Oliveira, Rodrigues, Early, Duarte, Buery, Bueno, Catão-Dias, Cerutti, Rona, Neafsey and Ferreira2021b), which encodes a reticulocyte-specific parasite ligand (Malleret et al. Reference Malleret, El Sahili, Tay, Carissimo, Ong, Novera, Lin, Suwanarusk, Kosaisavee, Chu, Sinha, Howland, Fan, Gruszczyk, Tham, Colin, Maurer-Stroh, Snounou, Ng, Chan, Chacko, Lescar, Chandramohanadas, Nosten, Russell and Rénia2021) and may be involved in the adaptation to binding to and/or entrance into red blood cells of platyrrhine monkeys. Our search also revealed additional putative genomic signatures of parasite–host-vector adaptation in P. simium: (1) an elevated nucleotide polymorphism in genes putatively involved in gametocytogenesis and ookinete and oocyst development (Supplementary Table 6) and (2) an IBD peak suggestive of a selective sweep in a genomic domain containing the gene encoding an AP2 domain transcription factor family member, AP2-SP2, which appears to modulate oocyst maturation in the vector (Figure 7). In addition, it has been suggested that adaptive changes in the gamete surface protein (P47) – the P. vivax ortholog of Pfs47, a surface protein that allows P. falciparum to evade the mosquito immune system (Molina-Cruz et al. Reference Molina-Cruz, Zilversmit, Neafsey, Hartl and Barillas-Mury2016) – could have enhanced the compatibility between P. vivax and P. simium and New World vectors. However, no IBD peak was found around the p47 gene and relatively few additional sequences are currently available for testing this hypothesis (de Oliveira et al. Reference de Oliveira, Rodrigues, Early, Duarte, Buery, Bueno, Catão-Dias, Cerutti, Rona, Neafsey and Ferreira2021b). Somewhat surprisingly, both P. simium and P. vivax from the Amazon appear to infect efficiently anophelines from across the globe, in addition to mosquitoes found in the Amazon and along the Atlantic Coast of South America (Collins et al. Reference Collins, Sullivan, Galland, Williams, Nace, Williams and Barnwell2005; Shaw-Saliba et al. Reference Shaw-Saliba, Clarke, Santos, Menezes, Lim, Mascarenhas, Chery, Gomes, March, Bhatia, Rathod, Ferreira, Catteruccia and Duraisingh2016).
The present study has some limitations. First, it was analysed genome sequence data from only 31 P. simium isolates from 3 states in Southeast Brazil. No WGS data were available from the southernmost range of the current distribution of P. simium (states of Santa Catarina and Rio Grande do Sul). Second, no WGS data was available from parasites labelled as P. vivax that were found to infect wild (Rondón et al. Reference Rondón, León, Link and González2019) or captive platyrrhine monkeys (Silva et al. Reference Silva, Barros, Bahia, Sampaio Junior, Santos, Inoue, Gonçalves, Chiesorin Neto, Faria, Tochetto, Viana, Monteiro, Góes-Cavalcante and Scofield2019) outside the currently known geographic range of P. simium. Whether these parasites share recent ancestry with P. simium lineages from Southeast Brazil remains to be investigated. Third, PCA and ADMIXTURE analyses of the genetic relatedness between P. simium and global or regional populations of P. vivax are largely exploratory and may be affected by geographic biases in parasite sampling. For example, no P. vivax sequence data from West and Central Africa, where P. vivax is rare but not absent (Baird, Reference Baird2022), were available for analysis. Importantly, P. vivax infections acquired in Angola – the origin of nearly two-thirds of the more than 4 million enslaved Africans displaced to Brazil over 3 centuries – have been occasionally described in migrants and travellers (e.g., Haiyambo et al. Reference Haiyambo, Uusiku, Mumbengegwi, Pernica, Bock, Malleret, Rénia, Greco and Quaye2019; Martins et al. Reference Martins, Marques, Nieto-Andrade, Kelley, Patel, Nace, Herman, Barratt, Ponce de León, Talundzic, Rogier, Halsey and Plucinski2020). Consequently, whether African lineages of P. vivax have contributed to the ancestry of P. vivax and P. simium lineages circulating nowadays in the New World, remains undetermined. In addition, isolates from the Western Amazonian state of Acre were overrepresented in our dataset of P. vivax sequences from Brazil (143 of 203, or 70.4% of the sequences analysed), which included data from only 23 isolates from the Eastern Amazonian states of Amapá and Pará.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182025100310.
Data availability
Reads of new WGS data described in this article have been deposited into the NCBI Sequence Read Archive and are publicly available under BioProject accession code PRJNA1242290.
Acknowledgements
We are grateful to all patients who donated blood samples for genome sequencing; to Ajucilene (Joice) G. Mota, Francisco Melo and their team at the Health Secretary of Mâncio Lima for their overall logistic support during fieldwork in Brazil; and to Juliana Tonini for her support in the laboratory and the field.
Author contributions
N.R.M.A., J.C.S. and M.U.F. conceived and designed the study. W.A.L., P.T.R. and T.C.O. gathered genome sequence data. N.R.M.A., R.J.S., A.D., T.C.S. and J.C.S. performed bioinformatic analysis. N.R.M.A., J.C.S. and M.U.F. wrote the first draft of the article, which was reviewed and approved by all authors.
Financial support
This research was supported by grants from the National Institute of Allergy and Infectious Diseases, National Institutes of Health, United States of America (grant U19 AI089681, subcontract to M.U.F.; grant R01 AI141900 to J.C.S.) and the Fundação de Amparo à Pesquisa do Estado de São Paulo, Brazil (FAPESP; grant 2016/18740-9 to M.U.F.). We also acknowledge post-doctoral FAPESP fellowships from FAPESP to N.R.M.dA. (2022/10056-2 and 2023/12394-5), P.T.R. (2018/03902-9) and T.C.d.O. (2021/01017-0); a doctoral FAPESP scholarship to W.A.L. (2023/15369-1); a senior researcher scholarship from the Conselho Nacional de Desenvolvimento Científico e Tecnológico of Brazil (CNPq; 301011/2019-2) to M.U.F.; and institutional support from the Fundação para a Ciência e Tecnologia of Portugal (FCT), through the projects UID/04413/2020 to the Global Health and Tropical Medicine Research Center and LA-REAL LA/P/0117/2020. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Competing interests
The authors declare there are no conflicts of interest.
Ethical standards
Clinical samples for genome sequencing were collected under a research protocol approved by the Institutional Review Board of the Institute of Biomedical Sciences, University of São Paulo and the National Committee on Ethics in Research of the Ministry of Health of Brazil (CAAE number, 6467416.6.0000.5467). Written informed consent was obtained from all study participants or their parents or legal guardians.
Open access
