Phylogenetic congruence of Plasmodium spp. and wild ungulate hosts in the Peruvian Amazon

Malaria parasites are known to infect a variety of vertebrate hosts, including ungulates. However, ungulates of Amazonia have not been investigated. We report for the first time, the presence of parasite lineages closely related to Plasmodium odocoilei clade 1 and clade 2 in free-ranging South American red-brocket deer (Mazama americana; 44.4%, 4/9) and gray-brocket deer (Mazama nemorivaga; 50.0%, 1/2). We performed PCR-based analysis of blood samples from 47 ungulates of five different species collected during subsistence hunting by an indigenous community in the Peruvian Amazon. We detected Plasmodium malariae/brasilianum lineage in a sample from red-brocket deer. However, no parasite DNA was detected in collared peccary (Pecari tajacu; 0.0%, 0/10), white-lipped peccary (Tayassu pecari; 0.0%, 0/15), and tapir (Tapirus terrestris; 0.0%, 0/11). Concordant phylogenetic analyses suggested a possible co-evolutionary relationship between the Plasmodium lineages found in American deer and their hosts.


Introduction
The broad genus Plasmodium infects diverse animal taxa, from birds, reptiles and mammals, with human and non-human primates being the most studied hosts (Perkins and Schaer, 2016).Limited efforts worldwide have reported that ungulates are associated with parasites of the genus Plasmodium (Asada et al., 2018;Boundenga et al., 2016;Kaewthamasorn et al., 2018;Martinsen et al., 2016;Templeton et al., 2016a).The only reported molecular sequences of parasites in ungulates have been tentatively named based on the few morphological descriptions reported for many decades (Templeton et al., 2016b), including: a) the multiple lineages of Plasmodium odocoilei are closely related and are grouped into two clades (clade 1 and clade 2), found in white-tailed deer (Odocoileus virginianus) in different locations in the United States (Guggisberg et al., 2018;Martinsen et al., 2016); b) Plasmodium bubalis lineages, termed type I and II, from water buffalo (Bubalis bubalis) in Thailand and Vietnam (Templeton et al., 2016a); and c) Plasmodium caprae, a lineage from a Zambian goat (Kaewthamasorn et al., 2018).
The only study of malaria-associated ungulate parasites from South American reported Plasmodium lineages related to P. odocoilei clade 2 found in wild populations of the pampas deer (Ozotoceros bezoarticus) in the Brazilian Pantanal (Asada et al., 2018).It is likely that the expansion of these ungulate parasite lineages coincided with the radiation of ungulates in the Southern hemisphere (Duarte et al., 2008).Population level molecular studies in South American ungulates could help clarify the evolutionary history malarial parasites in the Americas.This is particularly true in the Amazon, a complex ecosystem where multiple hosts of malaria-associated parasites coexist in wildlife and live in close proximity to human populations in malaria-endemic areas (de Abreu et al., 2019;Amir et al., 2018;Perkins and Schaer, 2016;Voinson et al., 2022).
We detected the Plasmodium spp.infection in wild populations of ungulates in a nearly pristine, isolated and malaria-endemic geographical area in the Peruvian Amazon, using molecular approaches to detect the mitochondrial DNA (mtDNA) of Plasmodium species.
It is important to highlight that the study of malaria parasites in wild populations is beset by logistical and financial restrictions, resulting in incomplete information regarding the diversity of their vertebrate hosts and environmental influences on their prevalence and distribution.Considering that subsistence hunting is a cultural activity and an important source of dietary protein for indigenous and non-indigenous populations throughout the Amazon (Jacob et al., 2023), we took advantage of blood samples collected from ungulates under sustainable hunting conditions by local indigenous communities.This methodology is appropriate to the geographical and sociocultural context of the study area and collaborations with the local population were highly advantageous to overcoming common limitations of field work in the region (Aysanoa et al., 2017;Cáceres et al., 2017;González-Olvera et al., 2022;Mayor et al., 2017).

Study area
Samples were taken from an area covering 107,000 ha of contiguous, predominantly upland forest along the Yavari-Mirin River in the northeastern Peruvian Amazon (Fig. 1).There is only one human community within the study area, the indigenous Yagua community of Nueva Esperanza (04°19′53"S; 71°57′33"W; UT-5:00), with a population of 343 inhabitants in 2020 (167 males and 176 females).The Yavari-Mirin River basin is 302 km from the nearest city, Iquitos, and has the highest wildlife diversity observed in the Peruvian Amazon, especially mammals, including 14 species of non-human primates.The Peruvian Ministry of Health currently considers this region a low but stable malaria endemic area (Bernárdez-Rodriguez et al., 2021).The villagers depend on subsistence activities (hunting, fishing, other forestry activities, and small-scale agriculture), but opportunistically trade wood, fish, wild meat, and agricultural products.

Samples
Local hunters were trained to collect blood samples from all species hunted for subsistence purposes using discarded material from legal hunting.From 2019 to 2020, blood samples were collected from 47 ungulates belonging to five species: red-brocket deer (Mazama americana, n = 9), gray-brocket deer (Mazama nemorivaga, n = 2), white-lipped peccary (Tayassu pecari, n = 15), tapir (Tapirus terrestris, n = 11), and collared peccary (Pecari tajacu, n = 10).Hunters impregnated FTA ® cards (Scheilcher & Schuell, Germany) with blood from the cranial or caudal vena cavaduring postmortem processing of each animal.Local hunters were also trained to record the date, species, and sex of all animals of all taxa hunted as part of their normal subsistence activities.There was no incentive or payment to encourage additional hunting of wildlife.This sampling strategy was integrated into a community participatory program to improve the conservation and sustainability of natural resources and indigenous livelihoods (González-Olvera et al., 2022).This approach allows for the determination of wildlife population health status without promoting hunting beyond subsistence levels.
Collected biological material was sealed in individual plastic bags with desiccant and stored in the collection area at room temperature for a minimum of two weeks and a maximum of six months before being transferred to −70 °C for preservation.DNA was extracted from blood spots on filter paper using the AllPrep DNA/RNA Mini Kit (Qiagen, Germany) according to the manufacturer's instructions at the Malaria Laboratory of the Institute of Global Health, Barcelona.The DNA was used for malaria parasite detection by nested PCR amplification, and the products were sequenced for a portion of the Plasmodium cytb gene.

Ethical review
The research protocols for the sampling of wild animals were approved by the Peruvian Forest and Wildlife Service (N° 258-2019-MINAGRI-SERFOR-DGGSPFFS) and the Institutional Animal Use Ethics Committee of the Universidad Peruana Cayetano Heredia (ref. 102,142).Dried wildlife blood samples on filter paper were exported from Peru to Spain with the approval of the Ministry of Agriculture and Irrigation (MINAGRI) through the Peruvian Forestry and Wildlife Service -SERFOR (N° 003258/SP, N° 003260/SP, N° 003568-SERFOR, N° 003579-SERFOR) according to the Nagoya Protocol.

Amplification of Plasmodium sequences by polymerase Chain reaction (PCR)
Ungulates blood spots samples on filter paper were screened for Plasmodium cytb sequences (~776 bp) as previously described (Martinsen et al., 2006), with some modifications to the PCR program conditions, as these were samples with low-density, submicroscopic infections (Guggisberg et al., 2018;Martinsen et al., 2016).Nested PCR was performed using DW2 [5 μM] 5-TAA TGC CTA GAC GTA TTC CTG ATT ATC CAG-3 and DW4 [5 μM] 5-TGT TTG CTT GGG AGC TGT AAT CAT AAT GTG-3 in the first round followed by genus-specific primers: FP3 [5 μM] 5-TAT ATA ACT TAT TTT TTG ATA TG-3 and RP3 [5 μM] 5-GTR ATW GCA TTA TCT GGA TGT GA-3 in the second round of PCR (Martinsen et al., 2006).For the first round, 5 μL of genomic DNA and 10 μL of 2× Platinum II Hot-Start PCR Master Mix (Invitrogen, USA) were used in a 20 μL reaction volume.Cycling conditions included an initial denaturation step of 94 °C for 3 min, followed by 40 cycles of denaturation (20 s at 94 °C), annealing (30 s at 62 °C), and elongation 1 min at 72 °C, followed by a final elongation step of 10 min at 72 °C.For the second PCR round, 2 μL of the primary PCR product was diluted 1/20 in DNAase and RNAase free water and used as template for nested PCR with 10 μL of 2× Platinum II Hot-Start PCR Master Mix (Invitrogen, USA) in a 20 μL reaction volume.Cycling conditions included an initial denaturation step of 3 min at 94 °C, followed by 40 cycles of denaturation (94 °C, 20 s), annealing (52 °C, 30 s) and elongation (72 °C, 1 min), followed by a final elongation step of 10 min at 72 °C.Amplified products were gel-purified using the QIA quick Gel Extraction Kit (Qiagen, Germany) and sequenced by Macrogen (Korea) to confirm Plasmodium infection using FP3 and RP3 primers (Protocol published in: https:// www.protocols.io/private/6653CE1D8DD611ECA1850A58A9FEAC02).
Positive controls for Plasmodium spp.were performed on DNA extracted from P. falciparum strain 3D7, and molecular grade water was included as a contamination control.The second round of PCR (nPCR) was performed in a separate laboratory from the first round.The bench, instruments and materials were disinfected with DNAZap (Invitrogen) after each PCR analysis.No pipettes or racks were shared between the two rounds of PCRs.

Phylogenetic analysis
The chromatograms for each sequence were corrected for errors associated with saturation at the beginning of the sequences and low signal at the ends.For amplicons which were weak, duplicates were sequenced and the forward and reverse direction sequences were assembled using the Geneious Prime software (v.2023.0.4.Biomatters Ltd.) (Kearse et al., 2012).Multiple sequence alignments were carried out using the CLUSTAL-OMEGA algorithm implemented in Geneious Prime (Thompson et al., 1994).We used different models of sequence evolution, which had the best fit to our dataset as inferred from the Akaike Information Criterion (AIC) value calculated in MEGA (v.11.0.13).Phylogenetic relationships were inferred using maximum likelihood (ML) methods using FastTree 2.1.11(Price et al., 2009) and PhyML (Guindon et al., 2010).Nodal support was assessed by bootstrap using 2000 pseudo-replicates.Additionally, Bayesian inference (BI) using the program MrBayes 3.2.6 (Huelsenbeck and Ronquist, 2001) implemented in Geneious Prime was also used to reconstruct phylogenetic relationships.Parameters for the evolutionary model were the same as in ML and the search was carried out in two simultaneous runs of one million generations, sampled each 1000 generations, with a burn-in of 25%.
The Bayesian Markov Chain Monte Carlo (MCMC) method implemented in the BEAST v2.7.4 package was used for divergence time calculations, with a molecular clock based on the mutation rate used for Plasmodium mtDNA of 1.5% per million years (Pacheco et al., 2011;Silva et al., 2011).General Time Reversal (GTR) was used as the initial evolutionary model, and an uncorrelated relaxed clock model.To follow one of the assumptions of the birth-death model, only one haplotype per lineage was allowed for molecular clock analysis to avoid overestimation of divergence times.A homologous sequence from Theileria taurotragi (ref.NC053926.1)was used as an outgroup and calibrated against the estimated time of divergence of the human malaria parasite P. malariae (~23 MYA) (Pacheco et al., 2011).All divergence estimates recovered with support values greater than 0.92 in all major branches.Phylogenetic analyses were based on ~776 bp of cytb sequences.We additionally analyzed 14 complete or nearly complete cytb genes, obtained from the NCBI nucleotide database) of: P. odocoilei strains obtained from samples of Anopheles punctipennis (accession number, KU133755.1 and KU133748.1),P. odocoilei strains recovered from white-tailed deer from USA (accession numbers, MG709249.1,KU133750.1,MG709243.1,MG709250.1,KU133754.1 and MH177860.1)and a P. odocoilei strain BRDeer-89 from Brazil (accession numbers, LC326034.1).In addition, we added existing sequences from Plasmodium species that cause malaria in ungulates.These include: Plasmodium bubalis strains of water buffalo (B.bubalis) from Thailand (accession number, LC090213.1 and OL672209.1)and Plasmodium caprae strain of Zambian goat (accession number, LC090215.1).Two Plasmodium species of humans and non-human primates, P. malariae (accession number, AB489194.1)and P. brasilianum strain Bolivian I (accession number, CM043783.1).
Reference genome sequence and annotation files were downloaded from GenBank.
Phylogenetic analysis of the partial Plasmodium cytb sequences from the brocket deer in the present study and sequences from North American P. odocoilei showed that 83.3% (5/6) of the sequences in the Amazonian brocket deer are related to North American ungulate parasite species.Sequences in two red-brocket deer were monophyletic with P. odocoilei clade 1 sequences obtained from North American white-tailed deer and Anopheles mosquitoes (with ML/BPP of 98/0.99,respectively), and sequences of three brocket deer (two red-and one gray-brocket deer) formed a monophyletic clade with P. odocoilei clade 2, with low support (54/0.96,respectively).Three, four and five nucleotides out of 979 bp (ref.Divergence estimates suggest P. odocoilei clade 1 and 2 lineages of the North and South American deer diverged approximately 2.6-3.5 MYA ago.The estimated divergence of lineages related to P. odocoilei clade 1 found in the red-brocket and the white-tailed deer was more recent (0.8-1.0 MYA), as was the time of divergence between P. odocoilei clade 2 lineages in the same host species (0.6-1.1 MYA).When other deer species were included, the divergence time estimate between P. odocoilei clade 2 lineages of the white-tailed deer and those in the Amazonian-brocket deer was 1.1-1.4MYA, and in the Amazonian-brocket deer and the Brazilian pampas deer, 1.4-1.9MYA (Fig. 3).

Discussion
The current study indicates that malaria parasite infection is frequent widespread in wild Amazonian-brocket deer (55.6%).These findings highlight features of the ecology of Plasmodium lineages closely related to P. odocoilei species that naturally infect Amazonian deer.First, studies suggested that these parasites are not host specific, as they can infect different species (Asada et al., 2018;Guggisberg et al., 2018;Martinsen et al., 2016;Templeton et al., 2016a).However, their host range is apparently restricted to New World deer, as the lineages have not been detected in any other ungulates.The parasites may have an ecogeographic range with multiple potential hosts, and therefore these Plasmodium spp.lineages are likely of concern in the Americas (Gutiérrez et al., 2017).
We report for the first time the presence of lineages related to clade 1 of P. odocoilei in South America, found in wild Amazonian red-brocket deer.In addition, three haplotypes related to clade 2 of P. odocoilei were identified in wild Amazonian red and gray-brocket deer.Until the present study, only lineages related to clade 2 of P. odocoilei had been described in South America in wild pampas deer (O.bezoarticus) from Brazilian Pantanal (Asada et al., 2018).These results are consistent with previous studies that identified up to seven different cytb alleles of P. odocoilei clade 2 in the white-tailed deer (O.virginianus) in North America (Guggisberg et al., 2018;Martinsen et al., 2016) and consistent with a monophyletic grouping of all ungulate malaria parasites (Templeton et al., 2016a).
The rate of divergence suggest that the P. odocoilei clade 1 and 2 lineages diverged approximately 2.6-3.5 MYA, a result consistent with that reported by Asada and collaborators (Asada et al., 2018).Therefore, one possible evolutionary scenario that would explain the existence of lineages related to both parasite clades in the American deer is that the divergence of these lineages occurred after American deer divergence from a common ancestor with other deer lineages (~5 MYA) (Duarte et al., 2008).The Great American Biotic Interchange -GABI (~3 MYA) between North America and South America (Carrillo et al., 2020;Domingo et al., 2020) may have been the event that separated these lineages and their parasites.However, due to the lack of fossil records and lack of reports of Plasmodium lineages from ungulates in South America it is currently not possible to test this hypothesis (Asada et al., 2018).
The estimated date of divergence of the parasitic lineages coincides with the invasion of at least eight distinct ancestral forms of land deer that migrated from North America to South America and have been grouped into two major clades: the gray clade, which includes the gray brocket (Mazama gouazoubira and M. nemorivaga), the marsh deer (Blastocerus dichotomus), the huemul (Hippocamelus bisulcus), the taruka (Hippocamelus antisensis), and the pampas deer (O.bezoarticus); and the red clade, which includes the species of the red deer group (Mazama bororo, Mazama nana, Mazama temama and M. americana) and the genus Odocoileus (Carrillo et al., 2020;Duarte et al., 2008).We hypothesize that the differentiation of P. odocoilei clade 1 lineages may have occurred in ancestral host species that gave rise to the "red clade" of American deer, which could indicate that these parasites have adapted to the red clade.However, the absence of infection of clade 1 lineages in "gray clade" hosts may reflect insufficient sampling and not necessarily the absence of the specific lineage in this group of ungulates.
An alternative hypothesis is that, given the longer divergence time of lineages within clade 2 relative to clade 1 of P. odocoilei (~1. 9 MYA), common ancestral parasite lineages of both clades in South America originated from clade 2 lineages found in a North American deer ancestor (Asada et al., 2018), possibly in separate transitions of hosts representing both the "red clade" and the "gray clade" of South American ungulates during the GABI.The radiation and adaptation of ungulates in the South American area could have led to the speciation of lineages belonging to both clades in South America.Furthermore, sequences found of the Plasmodium species in the "red clade" (M. americana and O. virginianus) show divergence within each of the clades (1 and 2) with a possible recent divergence (~1 MYA).This could reflect the separation of South American red deer during the Pleistocene (Domingo et al., 2020;Duarte et al., 2008), resulting in the isolation of these ungulates after their colonization of South America.
A sequence identical to P. malariae/brasilianum, was identified in red-brocket deer suggesting that ungulates are probably incidental hosts of other Plasmodium lineages, possibly through habitat sharing with other vertebrate hosts.However, further studies are needed to determine parasite host range (Chaves et al., 2022;Grigg and Snounou, 2017;Joyner et al., 2015;Silva et al., 2019).Even if only incidental hosts, it would be important to determine how often ungulates are sink hosts for simian malarias, given that frequent infections could facilitate the evolution of novel parasites in a group where there have been frequent host-shifts (Wolfe et al., 2007).The Plasmodium spp.lineages could be transmitted among different mammalian wildlife taxa by generalist vectors than by species with greater parasite specificity (Boundenga et al., 2016;Galen et al., 2018;Lutz et al., 2016).Such vector behavior has been proposed as an adaptive trait in response to temporal fluctuations in host diversity and density in forest environments such as Amazonia (Boundenga et al., 2016).However, this trait could increase the possibility of interspecific transfer of parasite lineages that have a propensity to infect different host species, which is a concern for global health.

Conclusions
We conclude that Plasmodium species from both clades 1 and 2 of P. odocoilei are present in Amazonian red-brocket deer.The phylogenetic history of Plasmodium lineages infecting South American deer may represent phylogenetic associations with their hosts.The estimated time of lineage divergence was 2.6-3.5 MYA, which coincides with the migration of artiodactyls from North to South America, and this event may have contributed to the presumed initial speciation of clade 1-related lineages and multiple clade 2-related lineages of the P. odocoilei parasite in South American deer.The possibility that ungulates are incidental hosts of other Plasmodium species cannot be excluded.
MH177860) and 607 bp (ref.KU133755 and KU133748), separated the red-brocket deer Plasmodium sequences from the P. odocoilei clade 1 sequences, respectively.Up to three related P. odocoilei clade 2 strains co-circulated in the Amazonian-brocket deer population studied.One of the sequences in the red-brocket deer differed by three nucleotides out of 419 bp from the P. odocoilei allele D lineage (ref.MG709250).A three nucleotides difference from the P. odocoilei allele B lineage (ref.KU133750) out of 607 bp was also observed.The Plasmodium sp.sequence in the only positive gray-brocket deer differed by three nucleotides out of 607 bp amplified from the P. odocoilei allele E/F lineage (ref.KU133754).One red-brocket deer Plasmodium sp.sequence was identical to Plasmodium malariae/brasilianum (Fig.2).

Fig. 1 .
Fig. 1.Map of study area, Nueva Esperanza indigenous community located on the slopes of the Yavari-Mirin River.

Fig. 2 .
Fig. 2. Phylogeny of Amazonian ungulate Plasmodium lineages.(A) Maximum likelihood (ML) tree for all ungulate Plasmodium lineages studied, obtained using PhyML with 2000 replicates.(B).Bayesian inference (BI) tree obtained using MrBayes, Bayesian posterior probability (BPP) is indicated for each internal branch.The lineages found in this study are shown in blue, red, and pink.The blue stripe represents clade 2 with lineages related to P. odocoilei reported in ungulates host: O. virginianus, O. bezoarticus, M. nemorivaga y M. americana.The red stripe indicates clade 1, which consists of the P. odocoilei-related lineages present in the host species: O. virginianus and M. americana, and the pink lineage indicates the only P. malariae/brasilianum lineage found in M. americana.Host silhouettes in which the studied parasite lineages were found is shown.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3. Time-calibrated phylogeny of maximum cladistic credibility based on mtDNA of Plasmodium species lineages from ungulates.TMRCA confidence intervals are shown at nodes.The blue stripe represents clade 2 with lineages related to P. odocoilei and the red stripe indicates clade 1 with lineages related to P. odocoilei.The reference lineages are those with accession numbers.Next to each lineage is the name of the host in which the parasite lineage was found.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)