Population structure of avian malaria parasites

Abstract The geographic distribution of genetic diversity in malaria parasite populations (Apicomplexa: Haemosporida) presumably influences local patterns of virulence and the evolution of host‐resistance, but little is known about population genetic structure in these parasites. We assess the distribution of genetic diversity in the partial Domain I of apical membrane antigen 1 (AMA1) in three mtDNA‐defined lineages of avian Plasmodium to determine spatial population structure and host–parasite genetic relationships. We find that one parasite lineage is genetically differentiated in association with a single host genus and among some locations, but not with respect to other hosts. Two other parasite lineages are undifferentiated with respect to host species but exhibit geographic differentiation that is inconsistent with shared geographic barriers or with isolation‐by‐distance. Additional differentiation within two other lineages is unassociated with host species or location; in one case, we tentatively interpret this differentiation as the result of mitochondrial introgression from one of the lineages into a second lineage. More sampling of nuclear genetic diversity within populations of avian Plasmodium is needed to rule out coinfection as a possible confounding factor. If coinfections are not responsible for these findings, further assessment is needed to determine the frequency of mitonuclear discordance and its implications for defining parasite lineages based on mitochondrial genetic variation. OPEN RESEARCH BADGES This article has earned an Open Data Badge for making publicly available the digitally‐shareable data necessary to reproduce the reported results. The data is available at Genbank https://www.ncbi.nlm.nih.gov/genbank/, accession numbers MK965548‐MK965653 and MK929797‐MK930264.


| INTRODUC TI ON
The ubiquity of avian malaria in the global avifauna makes it an accessible pathogen to investigate how the geographic distribution of parasite genetic diversity relates to patterns of virulence, the evolution of host-resistance, and the population structure and demography of hosts. Although more than 40 morphospecies of avian Plasmodium have been described (Valkiūnas, 2004) and hundreds of genetic lineages have been identified utilizing variation at the partial cytochrome b gene (Bensch, Hellgren, & Pérez-Tris, 2009), little is known of the genetic diversity within avian malaria parasite populations, or even concerning the boundaries of these populations. Whether individuals of a particular parasite lineage vary genetically relative to their hosts or to geographic location is poorly understood, but assessments of lineage distributions through time reveal dynamic patterns that differ among both location and host taxa Fallon, Ricklefs, Swanson, & Bermingham, 2003). Previous work in the West Indies has revealed disjunctions in the geographic distributions of parasite lineages (Ricklefs, Soares, Ellis, & Latta, 2016) and dynamic parasite community assembly processes, including differentiation of geographically separated parasite assemblages in as little as 2,500 years (Soares, Latta, & Ricklefs, 2017). However, information about the genetic diversity and structure of parasite populations, including effective population sizes and intra-lineage variation related to host species or geographic location, has not been reported. This lack of information reflects, in part, the difficulty of obtaining suitable genetic markers.
Genetic variation and population structures of two Plasmodium parasites infecting humans, P. vivax (Pv) and P. falciparum (Pf), have been characterized in detail, particularly at immunogenic loci. Among these loci is the apical membrane antigen 1 (AMA1), which plays a critical role in forming the junction between Plasmodium merozoites and host red blood cells (reviewed in Bai et al., 2005). AMA1 is also thought to interact directly with host immune systems; analyses of PfAMA1 have located one T-cell epitope and a second B/T-cell epitope in Domain I (Escalante et al., 2001), as well as several erythrocyte binding sites (Zakeri, Sadeghi, Abouie Mehrizi, & Dinparast Djadid, 2013). Analyses of variation in the three-dimensional structure of the protein suggest that the acquisition of several highly variable loops in Domain I is related to evasion of the host immune response (Collins, Withers-Martinez, Hackett, & Blackman, 2009). These direct interactions support Domain I of AMA1 as a marker related to host-specific immune pressures. Accordingly, the rapid accumulation of mutations as a result of diversifying selection provides information about fine-scale population structure. Assessments of PvAMA1 and PfAMA1 reveal distinct demographic patterns: PvAMA1 typically exhibits a differentiated structure consistent with an endemic pathogen (Neafsey et al., 2012;Taylor et al., 2013), while PfAMA1 most often exhibits an epidemic population structure with reduced diversity and little geographic differentiation (Arnott et al., 2014;Mueller, Kaiok, Reeder, & Cortés, 2002;Ord, Tami, & Sutherland, 2008), typical of frequent clonal outbreaks. Plasmodium vivax and P. falciparum exhibit these demographic differences despite infecting the same hosts in many of the same locations and being transmitted by the same vectors. Therefore, populations of avian malaria parasites might be expected to exhibit substantial variation in the degree of population structure within lineages as a result of the more complex relationships among hosts, vectors, and locations.
To develop a better understanding of the influence of host species and geographic location on parasite genetic diversity, we assess the distribution of genetic diversity in the partial Domain 1 of AMA1 of three mitochondrial lineages of avian Plasmodium commonly infecting birds in the West Indies and eastern North America. We assess genetic diversity and phylogenetic relationships within mitochondrial lineages and test whether parasite populations exhibit genetic structure related to hosts, location, or geographic distance.

| DNA extraction and sequencing
Samples for this study were obtained over several years from diverse localities in the Americas (see Figure 1). Birds were captured with mist-nets and ca. 10 μl of blood was collected by sub-brachial venipuncture (field techniques described in Latta & Ricklefs, 2010). DNA was extracted from blood using the isopropanol precipitation technique described in Svensson and Ricklefs (2009), and all samples were screened for avian malaria using primers 343F and 496R . Samples that screened positive were genotyped to identify the cytochrome b lineage of the infection (Bensch et al., 2009) using a variety of primers and PCR conditions described in Perkins and Schall (2002), Ricklefs et al. (2005), and Waldenström, Bensch, Hasselquist, and Östman (2004) of the partial Domain 1 of AMA1 were amplified using nested primers Pg_AMA1F1/Pg_AMA1R1 and Pg_AMA1F2/Pg_AMA1R2 and PCR protocols described in Lauron et al. (2014). Negative controls F I G U R E 1 Overall detections of lineages OZ01 (blue triangles), OZ04 (red diamonds), and OZ14 (green circles). OZ01 is broadly distributed in the Neotropics and Nearctic; OZ04 is primarily restricted to the West Indies with limited detection in North and South America; OZ14 is distributed in the eastern United States and locally in the Greater Antilles. All three lineages are known to cooccur on Jamaica, Hispaniola, and in the Ozarks region of Missouri were included in each PCR reaction, and products were verified by visualization on 1% TBE agarose gels with ethidium bromide. PCR product was cleaned using the ExoSAP-IT protocol (Bell, 2008) and sequenced by Eurofins Genomics (Louisville, KY). Contigs were aligned and edited in Mega6 (Tamura, Stecher, Peterson, Filipski, & Kumar, 2013), and chromatograms were checked by eye to confirm polymorphisms. Heterozygous positions were denoted with IUPAC ambiguous base codes, and haplotypes were reconstructing using PHASE in DNAsp v. 5 (Librado & Rozas, 2009).

| Statistical analyses of variation
Summary statistics describing genetic variation within parasite populations, based either on geography or host species, were calculated using DNAsp v.5 (Librado & Rozas, 2009). These statistics include number of haplotypes (h), haplotype diversity (Hd), and nucleotide diversity (π). Hd is the probability that two haplotypes chosen at random differ (Nei, 1987), while π is the average number of nucleotide differences per site for two sequences chosen at random (Nei & Li, 1979). We calculated π for the entire sequence, and to examine patterns of nucleotide substitution within the AMA1 gene, we additionally calculated π using a sliding window of 30 bp and step size of 9 bp (Lauron et al., 2014). Synonymous and nonsynonymous substitutions, and minimum number of recombination events (Rm, using the four-gamete test after Hudson & Kaplan, 1985), were also calculated in DNAsp v.5. Pairwise F ST values were calculated in Arlequin v. 3.5 (Excoffier & Lischer, 2010) for a reduced dataset containing only sequences from locations or hosts with 10 or more samples.

| Phylogeographic analysis
A median-joining haplotype network was generated using Popart (Leigh & Bryant, 2015) to compare and visualize relationships among haplotypes at each sample location and within each host species or family. Pairwise F ST values were calculated with confidence intervals based on 1000 permutations (Arlequin v.3.5, Excoffier & Lischer, 2010) for all locations with more than 10 sequences and for the two clusters identified in the haplotype networks of lineages OZ01 and OZ14. A Mantel test comparing pairwise F ST and geographic distance was employed to determine whether the population structure detected among some locations in lineage OZ01 is consistent with isolation-by-distance. The Mantel test was implemented in the R package "ade4" version 1.7-11 with 9,999 permutations.
Sliding window analysis of π revealed congruence in the distribution of diversity among lineages, particularly within the range of nucleotides from positions 1-130 ( Figure 4). We detected the highest polymorphism at nucleotide positions (nt) 94-130. This region aligns with PfAMA1 amino acid residues 259-271 (Escalante et al., 2001), which constitute a hypervariable region encompassing a Tcell epitope. Lineage OZ14 exhibits three additional peaks in π; one is shared with OZ04 (nt 204-249 in both lineages) and corresponds to a B/T-cell epitope in PfAMA1 at amino acid residues 279-288, a third is at OZ14AMA1 nt 150-170, and the final peak is at OZ14AMA1 nt 276-303. The latter two are not known to have immunogenic functions. These lineages share an amino acid insertion detected in other avian malaria parasites at PfAMA1 amino acid residue 188 (Lauron et al., 2014), and lineage OZ01 and a subset of OZ04 exhibit a second

| Phylogeographic analysis
Haplotype networks ( Figure 5) reveal that lineages OZ04 and OZ14 each contain a lineage-specific common and widespread AMA1 haplotype, and that OZ01 contains two common haplotypes. However, in OZ01 these haplotypes are shared with individuals belonging to the mitochondrial OZ04 lineage (visualized in Figure 5A). One of the OZ01 haplotypes corresponds to a group of parasites infecting many hosts in many locations, and a second

| D ISCUSS I ON
The avian malaria parasite lineages assessed here have been defined based on their similarity at the partial cytochrome b gene, but whether they represent good phylogenetic species has not been de- Outlaw . In the absence of morphological data, species delimitation is thought to be best addressed by comparing multiple gene trees (Bensch, Pérez-Tris, Waldenström, & Hellgren, 2004;Galen et al., 2018;. The AMA1 gene tree produced here contradicts the mitochondrial lineage assignment in some cases. Specifically, all parasites from mitochondrial lineage OZ01 possess AMA1 alleles similar to one another, but 31 parasites from mitochondrial lineage OZ04 also possess AMA1 alleles identical or very similar to OZ01 alleles; the remaining parasites of lineage OZ04 possess AMA1 alleles that are similar to one another and more than 60 mutations divergent from the OZ01 cluster ( Figure 6). We interpret these relationships as supporting TA B L E 2 Nucleotide diversity (π), percentage of synonymous variation (S%), and recombination estimates (Rm) for avian and human Note: Lineage OZ04 estimates of π are for all sequences, including suspected cases of introgression, major OZ04 includes only the major cluster, and minor OZ04 includes only the apparently introgressed samples. Locations of samples for congeneric comparisons: a Africa (for all sequences above, for a reduced dataset excluding seven sequences of uncertain taxonomic identity below, from Lauron et al., 2014), b Venezuela (Grynberg et al., 2008), c Brazil (Figtree et al., 2000), d Sri Lanka (Gunasekera et al., 2007), e Iran (Abouie Mehrizi, Sepehri, Karimi, Djadid, & Zakeri, 2013), f 10 locations distributed globally (Escalante et al., 2001), g Venezuela (Ord et al., 2008).  Bensch et al., 2009), and a morphological assessment of lineages OZ01 and OZ04 following the protocols and keys described in Valkiūnas, Iezhova, Loiseau, and Sehgal (2009) and Valkiūnas and Iezhova (2018) revealed no discernible difference between lineages.

OZ01
As a result, coinfections could not be identified or ruled out at present by inspecting blood smears. Though we cannot definitively rule out coinfections in these cases, we proceed in the following with a tentative interpretation of mitochondrial introgression among lineages and encourage increased attention to this matter in future research. The introgressed OZ01 parasites detected here were most often recovered in hosts that are rarely parasitized by OZ01 but commonly parasitized by OZ04, namely Coereba flaveola and Tiaris bicolor.
These patterns suggest that mitochondrial introgression from OZ04 to OZ01 may have facilitated this host shift. Moreover, the higher relative abundance of introgressed OZ01 parasites on Jamaica, the co-occurrence of OZ01 and OZ04 on the same island, and the relationship of the Jamaican haplotypes to the other introgressed haplotypes detected (Table 5, Figures 7 and 9) support Jamaica as the likely site of introgression. One possible hypothesis for this occurrence is that an individual mosquito vector on this island was infected by both OZ01 and OZ04 gametes (from either multiple blood meals or a single co-infected host) which hybridized during sexual reproduction in the mosquito vector. The hybrid progeny then repeatedly backcrossed with OZ01 parasites to produce the mitonuclear discordance we detect in contemporary populations.
The introgressed parasites exhibit several geographic disjunctions that might be related to migratory movement of hosts.
Specifically, one or more migrants might have carried the parasite from Jamaica to the North American Midwest (e.g., Missouri Ozarks), where it was transmitted to other hosts, which then carried it to Puerto Rico and the Lesser Antilles in following migrations, skipping the island of Hispaniola. Further dispersal within the Lesser Antilles may have been facilitated by movements of common endemic species, including C. flaveola and T. bicolor. Other intra-lineage divisions, for example in lineage OZ14 and in P. lucens (Lauron et al. (2014), discussed below), may also represent mitonuclear discordance resulting from introgression, although further sampling will be required to substantiate this. If substantiated, these findings imply that defining lineages based on CYTB similarity might often be problematic.
The comparatively high estimate of nucleotide diversity (π) for samples of P. lucens prompted our further investigation into the relationships among AMA1 haplotypes within that parasite species.
We constructed a median-joining haplotype network for P. lucens ( Figure 10)

| Congeneric comparisons
Plasmodium vivax and P. falciparum represent contrasting demographic patterns and impacts on hosts, and these parasites may provide an informative context in which to consider the findings presented here.

| CON CLUS IONS
Analyses at more loci and with wider taxon sampling will be necessary to uncover the complexities of these relationships, but findings here provide a glimpse into the host distribution, spatial distribution, and diversity of avian malaria AMA1 in natural host communities.
The complex spatial patterns and differentiation in relation to host genus described here suggest several possible influences on population structure, including host immune pressure, host dispersal, and migratory movements, as well as vector dispersal and host feeding preferences. Moreover, the mitonuclear discordance detected here warrants further investigation to assess the role of coinfections and to determine the frequency of such introgression events and implications for defining parasite lineages based on mitochondrial genetic variation.

ACK N OWLED G M ENTS
The author acknowledges Whitney R. Harris World Ecology Center, National Geographic Society, and the Curators of the University of Missouri for funding; Eldredge Bermingham, John Faaborg, Irby Lovette, Maria Pil, and Leticia Soares for help in the field.

CO N FLI C T O F I NTE R E S T
All authors declare no conflict of interest.

AUTH O R CO NTR I B UTI O N S
RER and MBH designed the research; RER collected samples; MBH and MTS performed laboratory and statistical analyses; MBH wrote the manuscript with contributions and revisions by all authors. All authors approved the final manuscript.

DATA ACCE SS I B I LIT Y
All sequences deposited to Genbank, accession numbers MK965548-MK965653 and MK929797-MK930264.

OPEN RESEARCH BADGES:
This article has earned an Open Data Badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at Genbank https ://www.ncbi.