INTRODUCTION
Ecological studies of avian haemosporidian parasites based largely on PCR amplification and direct sequencing of short segments of parasite mitochondrial genes have grown exponentially in number over the past decade (Bensch et al. Reference Bensch, Hellgren and Pérez-Tris2009). Single nucleotide substitutions have been used to identify unique lineages that are often used as discrete taxa for studies of geographical distribution, host specificity and parasite–vector associations. These studies have revealed an unprecedented amount of diversity among these organisms that is often difficult to reconcile with more traditional species concepts based on morphological and biological characteristics. As a rule of thumb, haemosporidians that differ by more than 5% in mitochondrial sequence appear to be distinguishable by morphological features of erythrocytic stages of the parasites (Hellgren et al. Reference Hellgren, Krizanauskiene, Valkiūnas and Bensch2007; Valkiūnas et al. Reference Valkiūnas, Iezhova, Loiseau, Smith and Sehgal2009), with most diversity captured by sequencing a 479 bp region of the cytochrome b gene (Hellgren et al. Reference Hellgren, Krizanauskiene, Valkiūnas and Bensch2007). In reality though, we know little about how to define these parasites as biological species because so little is known about their basic life cycles, vectors and epizootiology. This problem is illustrated by the remarkable diversity of mtDNA lineages that have been associated with Plasmodium relictum, one of the most widely distributed species of avian Plasmodium both in terms of host range and geographical distribution (Beadell et al. Reference Beadell, Ishtiaq, Covas, Melo, Warren, Atkinson, Bensch, Graves, Jhala, Peirce, Rahmani, Fonseca and Fleischer2006).
In an attempt to place isolates of P. relictum from the Hawaiian Islands into a global phylogenetic context, Beadell et al. (Reference Beadell, Ishtiaq, Covas, Melo, Warren, Atkinson, Bensch, Graves, Jhala, Peirce, Rahmani, Fonseca and Fleischer2006) identified 31 distinct lineages from around the world in two well-supported clades which included lineages identified as P. relictum by established morphological criteria. A single lineage of Plasmodium (lineage 15, previously identified as lineage GRW4, Bensch et al. Reference Bensch, Stjernman, Hasselquist, Ostman, Hansson, Westerdahl and Torres-Pinheiro2000) was identified in 75 infected introduced and endemic forest birds from Hawaii, Maui, Oahu and Kauai. The authors concluded that a single lineage of the parasite had been introduced to the islands and that the recent expansion of native forest bird populations in some low-elevation habitats (Woodworth et al. Reference Woodworth, Atkinson, LaPointe, Hart, Spiegel, Tweed, Henneman, LeBrun, Denette, DeMots, Kozar, Triglia, Lease, Gregor, Smith and Duffy2005) had likely not been facilitated by cryptic introduction of different parasite lineages of lower virulence. However, their report relied on results of traditional direct sequencing, which fails to detect related lineages making up less than 5–25% of genetically mixed infections (Kapoor et al. Reference Kapoor, Jones, Shafer, Rhee, Kazanjian and Delwart2004; Hunt et al. Reference Hunt, Fawcett, Carter and Walliker2005).
It has been suggested that traditional cloning and sequencing techniques be used when attempting to determine the presence of cryptic lineages in avian malarial infections (Perez-Tris and Bensch Reference Perez-Tris and Bensch2005). Such techniques are generally more efficient for the detection of minority variants in mixed infections than direct sequencing of PCR products (Lindner and Banik, Reference Lindner and Banik2009). Based on either cloning and sequencing or SNP analysis of partial nuclear trap (thrombospondin-related anonymous protein) genes of P. relictum, we have recently documented genetically diverse, mixed infections of P. relictum in the Hawaiian Islands (Jarvi et al. Reference Jarvi, Farias and Atkinson2008; Farias et al. Reference Farias, Atkinson, LaPointe and Jarvi2012). Based on these results, we hypothesized that mtDNA diversity in P. relictum in Hawaii may be higher than that detected by direct sequencing. This is of particular interest because co-evolutionary interactions among variants may affect both the virulence and transmission of the parasites across steep environmental gradients in the islands and within the diverse communities of native and non-native forest birds that occur along these gradients (Farias et al. Reference Farias, Atkinson, LaPointe and Jarvi2012). In recent years, next-generation amplicon sequencing has developed into an approach that offers several advantages over both direct sequencing or cloning and sequencing. We therefore re-examined diversity of mitochondrial lineages of P. relictum in Hawaii by using the Roche 454 platform for deep sequencing of a 346 bp fragment of the mitochondrial genome that encompasses a portion of the 479 bp region used by Hellgren et al. (Reference Hellgren, Krizanauskiene, Valkiūnas and Bensch2007) and other investigators. We illustrate here that while the predominant lineage appears to be GRW4, carefully confirmed rare lineages exist in Hawaii that are likely not due to multiple independent introductions.
MATERIALS AND METHODS
Sample collection and DNA extraction
Samples for this study were collected as part of a larger investigation of the ecology of avian malaria on the windward slopes of Mauna Loa and Kilauea Volcanoes (NSF DEB 0083944). Native and non-native forest birds were captured with mist nets on a monthly sampling schedule (minimum 2–4 days/month/site) at nine 1 km2 study sites (Fig. 1) in 2002, 2003 and 2004 (Woodworth et al. Reference Woodworth, Atkinson, LaPointe, Hart, Spiegel, Tweed, Henneman, LeBrun, Denette, DeMots, Kozar, Triglia, Lease, Gregor, Smith and Duffy2005). Birds were measured, sexed and aged, and a blood sample was obtained by jugular venipuncture. Birds were examined for the presence or absence of Avipoxvirus-like lesions on the feet, legs, beak and around the eyes at the time of capture and infection status was recorded. Samples were screened for infection with P. relictum using a combination of serology and microscopy as previously described (Woodworth et al. Reference Woodworth, Atkinson, LaPointe, Hart, Spiegel, Tweed, Henneman, LeBrun, Denette, DeMots, Kozar, Triglia, Lease, Gregor, Smith and Duffy2005). Genomic DNA was extracted from packed red blood cells stored in lysis buffer at −80 °C using a Qiagen DNeasy Tissue Extraction Kit following manufacturer's protocols. DNA concentrations were determined using a NanoDrop ND-1000 Spectrophotometer.
Map of eastern Hawaii Island showing locations of 1 km2 study sites. Solid lines indicate 250 m elevation contours. The high-elevation site (>1500 m above sea level) is CJ Ralph (CJR). The mid-elevation sites (1000 to 1300 m) are Ainahou Ranch (AIN), Crater Rim (CRA), Cooper's (COO), Puu Unit (PUU) and Waiakea (WAI). The low-elevation sites (<300 m) are Bryson's (BRY), Nanawale (NAN) and Malama Ki (MAL).
PCR amplification and 454 sequencing
We designed and optimized 454 fusion primers (Roche, 2009) for preparation of amplicon sequencing of the cytochrome b gene of P. relictum. A total of 14 forward and 14 reverse primers each incorporated one of 14 unique 10 nt-long multiplex identifier (MID) sequences. The cytochrome b forward primer basic sequence is: 5′-CGTATCGCCTCCCTCGCGCCATCAG(MID)AAAGCACATTTAATTCATTATCC-3′, and the reverse primer: 5′-CTATGCGCCTTGCCAGCCCGCTCAG(MID)TAGCACCCCAGAAACTCAT-3′. The sequencing primer (A or B) is upstream of the MID, and the cytochrome b specific sequence is downstream of the MID. A unique combination of the forward and reverse primers was used for amplification from each individual to allow for sequence assignment. The target-specific portion of the reverse fusion primers fall within the region of cytochrome b analysed by Hellgren et al. (Reference Hellgren, Krizanauskiene, Valkiūnas and Bensch2007), so that the 3′ ends of the sequences reported here overlap the 5′ ends of sequences reported in the MalAvi database by 161 nucleotides. To produce tagged amplicons, 1–7 μL gDNA were used as template in a 25 μL reaction containing 1X reaction buffer, 1 mm MgCl2, 200 μ m each dNTP, 0·8 μ m each forward and reverse fusion primers, and 1·25 units FastStart High-Fidelity Enzyme Blend (Roche). Successful amplicons were isolated on a 2% low-melt agarose gel (Fisher), excised, and purified using the NucleoSpin Extract II Kit (Machery-Nagel) according to manufacturer's protocols. Purified amplicons were quantified either by BioSpec-Nano (Shimadzu) or PicoGreen (Life Technologies) analysis with signal detection on a StepOne Plus Real Time PCR system (Life Technologies). Equimolar amounts of each amplicon were mixed for sequencing with a maximum possible of 196 uniquely tagged individuals per 1/4 sequencing run. Sequences were run at ASGPB, University of Hawaii at Manoa on a Roche Next Generation Sequencer GS FLX system.
Data management and analysis
Initial parsing of sequence data was completed using a new, user-friendly and publicly available portal called Integroomer. The Integroomer portal was developed specifically to handle the requirements of the current project but was designed and implemented with the purpose of flexibly demultiplexing and deconvoluting any DNA construct. The pipeline was designed with performance in mind and can easily adapt to run in a grid or cloud environment or be integrated into the Galaxy platform (Blankenberg et al. Reference Blankenberg, Von Kuster, Coraor, Ananda, Lazarus, Mangan, Nekrutenko and Taylor2010). Integroomer was developed in Python using the Django framework (www.djangoproject.com). A free instance of the Integroomer pipeline can be used at http://courge.ics.hawaii.edu/integroomer. For this study, the Roche primers A and B were trimmed using the sffInfo script in the SFF Tools suite (Roche). The reads were then demultiplexed and parsed using Integroomer.
In Integroomer, a construct is defined as a series of hyphen-delimited labels representing the adapters and a specific region of interest, termed GENSEQ, to be extracted. The labels comprising a construct are either user- or keyword-defined. User-defined labels represent physical adapters whose sequences are fixed such as MIDs and target specific priming sequences, while keyword-defined labels represent subsequences of unknown length and/or value. The keyword GENSEQ denotes the region of interest, while the REGEX and UNKNOWN keywords are used to describe a class of subsequences or to simply refer to sequences of unknown length or value. As an example, all cytochrome b amplicon sequences originating from a single individual and amplified using forward MID01 and reverse MID01 would be parsed with the following construct:
KEY-MID_01_FOR-CYTB_A_FOR-GENSEQ-CYTB_B_REV-MID_01_REV-UNKNOWN where CYTB_A_FOR and CYTB_B_REV refer to the target-specific portion of the forward and reverse fusion primers given above, and the label KEY refers to the Roche sequencing key (TCAG). The postfix substring is encoded using the reserved label UNKNOWN, allowing for reads which continue beyond the reverse primer for a varying number of nucleotides. The output of running this construct in Integroomer is an .sff file containing all reads matching the specified adapters and criteria.
To allow for the highest level of flexibility in the parsing process, adapters are assigned attributes to describe how they can be matched to the sequencing reads. These attributes are: (1) the direction, used to determine whether the adapter should be matched in the forward or reverse/complement orientation; (2) a Boolean value describing whether the occurrence of the adapter is required; (3) minimum overlap length; (4) the maximum number of allowed errors; and (5) whether the adapter is required to align to the end of the sequence (blunt alignment), or whether the sequence is allowed to continue beyond the adapter. Given that users specify the stringency and minimum quality required for every hit, Integroomer automatically filters the reads whose lengths are inferior to the total sum of the minimal required adapter overlap lengths.
Integroomer output files containing parsed sequences were converted to fastq files using the Python script SFF_extract (http://bioinf.comav.upv.es/) and quality filtered using PRINSEQ (Schmieder and Edwards, Reference Schmieder and Edwards2011) and Sequencher 5.0 (GeneCodes). Sequences containing errors in the priming sites, ambiguous bases, or indels were discarded, as were sequences with a mean Q score of less than 30 following trimming of primers and sequences in which more than 10% of nucleotides had a Q score of 20 or less. To determine the probability that the observed polymorphisms were due to sequencing artefacts, each nucleotide, i, was assigned a phred quality score, Q, which was used to infer its probability P(i) of being erroneously called. The probability P(i) was computed as follows (Ewing and Green, Reference Ewing and Green1998): P(i) = 10 −Q/10. Considering that the phred quality scores of n bases belonging to the same column of a multiple sequence alignment are independent, we can infer P(variant), the probability of variant being erroneous as:
$$P(variant) = \prod\limits_{i = 1}^n {P(i)} $$Data were analysed using Sequencher 5.0 for alignment, and Minitab® 16.2.1 for statistical analysis. Phylogenetic and molecular evolutionary analyses such as construction of a neighbour-joining tree, average nucleotide divergence and nucleotide content were conducted using MEGA version 4 (Tamura et al. Reference Tamura, Dudley, Nei and Kumar2007) with the gamma parameter estimated in MEGA version 5 (Tamura et al. Reference Tamura, Peterson, Peterson, Stecher, Nei and Kumar2011).
Non-random distribution of variants
Since Avipoxviruses are known to be immunosuppressive (Smith and Kotwal, Reference Smith and Kotwal2002), we evaluated mean variant frequency (# variants/total number of hosts) as well as variant prevalence (# hosts with the variant/total # hosts) in Apapane and Amakihi with and without Avipoxvirus-like lesions and tested for mean differences with a two-sample t-test (Minitab® 16.2.1). If variants are real rather than artefactual, we hypothesized that infection with Avipoxvirus would lead to significant changes in variant frequency and prevalence. We classified lesions as Avipoxvirus-like based on morphological criteria used in previous studies to identify active lesions (van Riper et al. Reference van Riper, van Riper and Hansen2002).
RESULTS
Cytochrome b sequences from parasite populations from a total of 97 malaria-positive Hawaii Amakihi and 42 malaria-positive Apapane were evaluated in this study. Samples were collected across an altitudinal gradient of 3 low-elevation (<300 m), 3 mid-elevation (1200 m) and 1 high-elevation (1500 m) study sites (Fig. 1, Table 1). Because abundance of both native species varies with elevation, we were not able to obtain balanced sample sizes of each species at each elevational zone. Blood samples were collected from Hawaii Amakihi at 3 low-elevation and 3 mid-elevation study sites with most (84%, 81/97) of the samples originating from 3 low-elevation sites. Apapane samples were collected from 2 low-, 3 mid- and 1 high-elevation study sites with 90% (38/42) of the samples originating from mid-elevation sites. A total of 26 238 reads were obtained and after quality filtering a total of 16 636 sequences were further evaluated. Sequences that appeared only once (singlets) were removed, as were sequences from individuals that had <10 usable sequences/individual. Also removed from the final dataset were loci that did not meet the 95% criterion rule for defining polymorphic loci in populations described in Hartl and Clark (Reference Hartl and Clark1997), thus only variants that were detected in ⩾7 individual hosts (⩾5% of the host population of 139) were included in this study.
Total number of GRW04 and variant sequences detailed by species and geographic site of origin, and calculated probability that the variant sequences are due to sequencing error. Site name abbreviations are: BRY (Brysons), MAL (Malama Ki), NAN (Nanewale), PUU (Pu`u unit), AIN (Ainahou), CRA (Crater Rim), COO (Coopers), WAI (Waikea) and CJR (CJ Ralph)
a nv, Total number of individuals with variant.
Our final dataset consisted of 15 953 reads of which 15 742 were GRW4 (Table 1). A total of 211 variant sequences were further analysed. Sequences were imported into MEGA 4·0 and in a neighbour-joining tree (Kimura 2-parameter) 23 variant clusters each containing identical sequences originating from 8–15 individuals per cluster were delineated (not shown). The numbers of sequences found at each geographical site and in each species for the 23 variants and GRW4 are described in Table 1. All variants are single nucleotide transition substitutions.
For each variant, the quality scores of the bases belonging to its multiple sequence alignment were used in equation (1) (see methods) to compute the probability that the variant is the product of sequencing error (right column, Table 1). The extremely low probabilities obtained are clear evidence that the observed polymorphisms are highly likely not sequencing artefacts.
All of the 139 avian hosts possessed GRW4. A total of 49 hosts possessed only GRW4 (35%), GRW4 plus 1 additional variant were detected in 34 hosts (24%), 2 additional variants in 27 hosts (19%), 3 additional variants in 14 hosts (10%), 4 additional variants in nine hosts (6·5%), 5 additional variants in 5 hosts (3·6%), and 6 additional variants were detected in 1 host (0·7%). A total of 132 variant sequences were obtained from 121 Hawaii Amakihi and 79 variant sequences from 72 Apapane. The hypothesized amino acid sequences show that nucleotide substitutions in 9 of the 23 variants are non-synonymous resulting in amino acid changes, 11 are synonymous substitutions. Three substitutions resulted in stop codons and these sequences were eliminated from further analysis. A summary of nucleotide and amino acid substitutions in all variants as compared with GRW4 is presented in Table 2. Average nucleotide divergence was estimated at 0·006±0·001SE (Kimura 2-parameter, г-corrected with the gamma parameter estimated in MEGA 5, standard error estimate obtained by bootstrap with 5000 replicates) and was based on the number of base substitutions per site (n = 346) over all sequence pairs (n = 21, i.e. 20 variants and GRW4). Total G–C content was 25·3%. All sequences are available in GenBank accession nos KF428723–KF428746, or as an alignment file upon request.
Summary of nucleotide (nt) and amino acid (aa) substitutions in all variants as compared with GRW4
a N designates non-synonymous substitutions, S designates synonymous substitutions.
b Variants resulting in stop codons.
Variant frequency and variant prevalence comparisons were made among Apapane and Hawaii Amakihi with or without Avipoxvirus-like lesions. Variant frequency (#variants/#hosts) was higher (P = 0·006) in malaria-infected Amakihi with lesions (n = 15, mean = 2·07), than in Amakihi with no lesions (n = 57, mean = 0·88). However, no significant differences (P = 0·691) were detected between Apapane that were lesion positive (n = 23, mean 1·52) and those that were lesion negative (n = 16, mean 1·31). Similarly, variant prevalence (# of individual hosts harbouring a variant/the total number of hosts) in malaria-infected Amakihi with lesions (mean 0·1033) was significantly higher (P = 0·001) than in Amakihi with no lesions (mean 0·0456) (Fig. 2). However, no significant differences (P = 0·497) were detected between Apapane that were lesion positive (mean 0·0761) and those that were lesion negative (mean 0·0656) (data not shown).
Average variant prevalence in Amakihi with Avipoxvirus-like lesions appears greater than in Amakihi without lesions (P = 0·001). Prevalence was calculated as the number of individual hosts harbouring that variant/the total number of hosts.
DISCUSSION
We document previously undetected levels of variation of P. relictum in Hawaii based on next-generation sequencing of a 346 bp fragment of the cytochrome b gene from 97 Hawaii Amakihi and 42 Apapane. Until now, direct sequencing of Plasmodium mtDNA fragments (cytochrome b) from 75 Hawaiian forest birds identified only the most common lineage (GRW4) (Beadell et al. Reference Beadell, Ishtiaq, Covas, Melo, Warren, Atkinson, Bensch, Graves, Jhala, Peirce, Rahmani, Fonseca and Fleischer2006). The application of new sequencing methodologies that do not require cloning and allow for increased depth of coverage provide a viable alternative and, as expected, allow for more thorough evaluation.
We detected 23 additional lineages based on single nucleotide transitional changes at 23 distinct sites. Single nucleotide polymorphisms (SNPs) appearing in only one or two individuals might be attributable to PCR or sequencing error, but each variant reported here was detected in 7–15 different individuals. These SNPs were not concentrated in any one region, nor did they occur in regions of high G–C or A–T content. For additional verification, we calculated the probabilities that the SNPs defining each variant were due to sequencing artefacts. The extremely low probabilities (in Table 1) provide clear evidence that these polymorphisms are highly unlikely to be due to sequencing error. Small differences introduced by singlet polymorphisms in the cytochrome b gene have been shown to have significant biological consequences. For example, in Plasmodium falciparum, resistance to atovaquone has been associated with specific point mutations in the cytochrome b gene in the region spanning codons 271–284 (Korsinczky et al. Reference Korsinczky, Chen, Kotecka, Saul, Rieckmann and Cheng2000; Schwobel et al. Reference Schwobel, Alifrangis, Salanti and Jelinek2003; Ekala et al. Reference Ekala, Khim, Legrand, Randrianarivelojosia, Jambou, Fandeur, Menard, Assi, Henry, Rogier, Bouchier and Mercereau-Puijalon2007). Similarly minor changes have been documented as multiple lineages that differ from GRW4 by one nucleotide (PADOM10, MYIALE02, PRBQ22, PRBQ21), two nucleotides (LINOL101, LINOL102) or five nucleotides (PRBQ20, Bensch et al. Reference Bensch, Hellgren and Pérez-Tris2009). Within this framework, the sequences reported here may therefore be considered unique mitochondrial lineages.
We designed our primers to amplify a cytochrome b region that appeared more variable than the 479 bp fragment analysed by Hellgren et al. (Reference Hellgren, Krizanauskiene, Valkiūnas and Bensch2007). Based on previous cloning studies of a 1193 bp fragment from six P. relictum isolates from four native Hawaiian species on the islands of Kauai and Hawaii (54 sequences), a greater number of polymorphic sites was found outside of this region than within, especially upstream of the 5′ primer of this 479 bp fragment (Jarvi et al. unpublished results). The 3′ end (161 bp) of the 346 bp fragment we amplified overlaps with the 479 bp fragment. Thirteen of the 23 polymorphic sites were found within the overlap while we identified an additional 10 lineages, five of which resulted in non-synonymous substitutions that would have been missed by focusing on the 479 bp region alone.
The rare, distinct lineages reported here likely either evolved from the ancestral type (i.e. GRW4) via mutation events or are the result of multiple introductions of parasites to the Hawaiian Islands. Multiple lineages may have been introduced to Hawaii either through introduction of multiple infected hosts or through introduction of a few individual hosts infected with multiple parasite lineages. In the case of multiple introductions, GRW4 clearly appears highly competitive. Based on currently available evidence, P. relictum has been in Hawaii for no more than 100 years (van Riper et al. Reference van Riper, van Riper, Goff and Laird1986). Ricklefs and Outlaw (Reference Ricklefs and Outlaw2010) report a rate for cytochrome b evolution in avian malaria parasites that corresponds to ∼1·2% sequence divergence per million years. If one presumes a steady rate of evolution over time, one might then expect 0·00012% sequence divergence per 100 years. This is much lower than our estimate of average evolutionary divergence of 0·006 (0·001) or 0·6% among GRW4 and the 20 variants (excluding variants possessing stop codons, i.e. potential pseudogenes). If the variation observed in this study is presumed due to mutations of a single introduced lineage (i.e. GRW4), the observed level of divergence in this 346 bp fragment having occurred within this time period seems highly unlikely based on predicted divergence rates. It could be that the cytochrome b gene in P. relictum in this disease system is evolving more rapidly. This might be expected given the year-round transmission, the relatively high levels of parasitaemias, and the lower selective pressure placed on the parasites by susceptible native Hawaiian forest birds. Also, co-infection with Avipoxvirus can result in immunosuppression in avian hosts, thus potentially further reducing selective pressure on parasite populations.
We found higher mean variant frequency and prevalence in low-elevation Hawaii Amakihi co-infected with Avipoxvirus, which provides additional support that these variants are likely not artefacts and that they may have distinct biological characteristics. Poxviruses possess genes that especially modify the innate mechanisms of the host immune response, therefore sustaining viral persistence during early infections (Smith and Kotwal, Reference Smith and Kotwal2002). A release of immunological regulation of Plasmodium parasites in chronic malaria infections due to Avipoxvirus co-infection may result in recrudescence of malaria with potential impacts on malaria diversity, which could provide an explanation for our finding of a higher mean variant prevalence in Hawaii Amakihi with Avipoxvirus-like lesions than in those without lesions. These data mirror earlier preliminary (unpublished) studies of P. relictum diversity in low-elevation Hawaii Amakihi and mid-elevation Apapane based on SNP analysis of the nuclear trap gene of P. relictum. The trap gene is a single copy gene in Plasmodium spp. (Robson et al. Reference Robson, Hall, Davies, Crisanti, Hill and Wellems1990), thus the detection of multiple trap sequences within a single host reflects multiple variants. We found the rate of multiple trap infection was much higher in malaria-infected Hawaii Amakihi with lesions (70%) as compared with malaria-infected Hawaii Amakihi without lesions (46%), while the rates of multiple trap infection were very similar in malaria-infected Apapane with lesions (39%) versus those infected with malaria only (35%). We know from experimental studies that birds exhibit premunition to P. relictum, i.e. exposure to one malarial isolate provides immunological protection against a second isolate (Atkinson et al. Reference Atkinson, Dusek and Lease2001). It is unlikely then, that chronically infected hosts would acquire new variants. Upon co-infection with Avipoxvirus, suppression of the host immune response may lead to recrudescence of existing infections, and/or possibly allow for acquisition of new variants through mosquito-borne transmission. If co-infection with Avipoxvirus results in increased malaria diversity through immunosuppression mechanisms, we would expect the effects to be most evident in populations or species exhibiting stronger immunity to malaria, i.e. more resistant population such as low-elevation Amakihi. Our preliminary data strongly suggest that this is the case, at least in low-elevation Hawaii Amakihi. We do not see the same effect upon co-infection of mid-elevation Apapane. Malaria and pox transmission levels vary between elevations, being year-round and endemic at low elevations and seasonal at mid-elevations (Atkinson and LaPointe, Reference Atkinson, Lapointe, Pratt, Atkinson, Banko, Jacobi and Woodworth2009). Malaria diversity based on trap is greater at low elevations than mid-elevations (Farias et al. Reference Farias, Atkinson, LaPointe and Jarvi2012), which is consistent with expectations of greater parasite diversity in areas of higher transmission (Conway, Reference Conway2007). Species differences may play a role in the increase in malaria diversity observed upon co-infection of low-elevation Hawaii Amakihi but not mid-elevation Apapane. However, in this comparison, species effects might be confounded by potential elevation-based transmission effects and cannot be determined at this time.
A total of 35% of hosts tested possessed only GRW4 and 65% of hosts possessed GRW4 along with at least one other variant, but no hosts harboured variants in the absence of GRW4. If multiple, independent introductions were responsible for the observed variation, one might expect to find variants independent of GRW4. To date, the mtDNA genomes of P. yoelii, P. falciparum, P. vivax, P. gallinaceum and P. juxtanucleare appear to be ∼6 kb, and exhibit a linear, tandem, repeated structure (Vaidya and Arasu, Reference Vaidya and Arasu1987; Feagin, Reference Feagin1992; Sharma et al. Reference Sharma, Pasha and Sharma1998; Omori et al. Reference Omori, Sata, Isobe, Yukawa and Murata2007). Given the tandemly repeated structure of the Plasmodium mitochondrial genome and the finding of three variants with stop codons, it is possible that some of these GRW4 concatenated repeats might contain cytochrome b variants. Another possibility is that these variants exist with GRW4 in heteroplasmy. Lastly, it is also possible that these variants are nuclear mtDNA (NUMTs) copies, however this is least likely since Plasmodium have low numbers of NUMTs (Richly and Leister, Reference Richly and Leister2004) and only one mitochondrion per merozoite during early stages of the asexual erythrocytic cycle (Divo et al. Reference Divo, Geary, Jenson and Ginsburg1985; Hopkins et al. Reference Hopkins, Fowler, Krischna, Wison, Mitchell and Bannister1999). These scenarios all offer viable explanations for the lack of detection of variants independent of GRW4.
While the source of variation has yet to be documented, multiple, independent introductions seem unlikely. Instead, we hypothesize that these variants, whether they exist within the concatenated repeat structure of the Plasmodium mitochondrial genome or exist in heteroplasmy, were either introduced to Hawaii with GRW4 or arose due to an exceptionally rapid mutation rate potentially occurring during replication of the mtDNA genome. Whole genome sequencing of these variants might further our understanding of the nature of mitochondrial diversity. A comparison of next generation sequencing of parasites in ancient versus modern-day populations as well as examination of the geographical range and diversity of these lineages across the length of the archipelago with continental populations may further our understanding of the origin and evolution of this parasite in the Hawaiian Islands.
ACKNOWLEDGEMENTS
We would like to thank Staffan Bensch, Robert Ricklefs and Andrew DeWoody for early manuscript review. We thank Julie Lease, Matt Reiter and Becky Imdieke (Ainahou Ranch) and Pat Hart, Bethany Woodworth, Erik Tweed, Carlie Henneman, Jaymi LeBrun, Tami Denette and field interns of the Biocomplexity of Avian Diseases Project (all other sites) for collecting blood specimens.
FINANCIAL SUPPORT
Funding for this project was provided by Biocomplexity of Introduced Avian Diseases in Hawaii: Threats to Biodiversity of Native Forest Ecosystems (NSF DEB 0083944) (SIJ, CTA, DAP), NIH/NCRR IDeA Networks of Biomedical Research Excellence (INBRE) the National Institute of General Medical Sciences (8 P20 GM103466-11) (SIJ) and the National Center for Research Resources (5P20RR016467-11)(SIJ), and P20GM103516 (MB) from the National Institutes of Health, the National Park Service (CTA, DAP), Natural Resource Protection Program (NRPP) (CTA, DAP), and the US Geological Survey Wildlife and Terrestrial Resources Program (CTA, DAP). Any use of trade, product or firm names in this publication is for descriptive purposes only and does not imply endorsement by the US Government. The content is solely the responsibility of the authors and do not necessarily represent the official views of the National Institutes of Health.
