Contrasting the seasonal and elevational prevalence of generalist avian haemosporidia in co‐occurring host species

Abstract Understanding the ecology and evolution of parasites is contingent on identifying the selection pressures they face across their infection landscape. Such a task is made challenging by the fact that these pressures will likely vary across time and space, as a result of seasonal and geographical differences in host susceptibility or transmission opportunities. Avian haemosporidian blood parasites are capable of infecting multiple co‐occurring hosts within their ranges, yet whether their distribution across time and space varies similarly in their different host species remains unclear. Here, we applied a new PCR method to detect avian haemosporidia (genera Haemoproteus, Leucocytozoon, and Plasmodium) and to determine parasite prevalence in two closely related and co‐occurring host species, blue tits (Cyanistes caeruleus, N = 529) and great tits (Parus major, N = 443). Our samples were collected between autumn and spring, along an elevational gradient in the French Pyrenees and over a three‐year period. Most parasites were found to infect both host species, and while these generalist parasites displayed similar elevational patterns of prevalence in the two host species, this was not always the case for seasonal prevalence patterns. For example, Leucocytozoon group A parasites showed inverse seasonal prevalence when comparing between the two host species, being highest in winter and spring in blue tits but higher in autumn in great tits. While Plasmodium relictum prevalence was overall lower in spring relative to winter or autumn in both species, spring prevalence was also lower in blue tits than in great tits. Together, these results reveal how generalist parasites can exhibit host‐specific epidemiology, which is likely to complicate predictions of host–parasite co‐evolution.


| INTRODUC TI ON
Epidemiological dynamics and virulence evolution are dependent on the selective pressures that parasites face across their infection landscape (Alizon, Hurford, Mideo, & Van Baalen, 2009;Mackinnon & Marsh, 2010;Rigaud, Perrot-Minnot, & Brown, 2010). These pressures will vary dependant on the host species infected (as hosts differ in their physiology, immunology, behavior, and distribution [Innes, 1997;Li, Cowles, Cowles, Gaugler, & Cox-Foster, 2007;Palinauskas, Valkiunas, Križanauskiene, Bensch, & Bolshakov, 2009;Perkins & Swayne, 2001]), and on seasonal and spatial variation in abiotic conditions; some effects of which may be mediated through the host (Amundsen, Knudsen, Kuris, & Kristoffersen, 2003;Cosgrove, Wood, Day, & Sheldon, 2008;Tack, Thrall, Barrett, Burdon, & Laine, 2012). In parasite ecology, the existence of spatial and temporal structuring of parasite prevalence is well established (Altizer et al., 2006;Penczykowski, Laine, & Koskella, 2016). Climatic variability, such as over an elevation gradient or between seasons, is a particularly common driver of parasite prevalence because abiotic conditions (e.g., temperature) can influence parasite transmission and development (Harvell et al., 2002;LaPointe, Goff, & Atkinson, 2010;Sehgal, 2015). For example, sporogony (development of infective sporozoites) of the protozoan blood parasite Plasmodium relictum within its mosquito vector was shown to significantly decrease below 21°C (LaPointe et al., 2010). In addition, seasonal variation in infection risk coinciding with key phenological events, such as the onset of host reproduction and vector emergence, can give rise to temporal variation in prevalence (Altizer et al., 2006;Cosgrove et al., 2008). Rodent tick-borne encephalitis infection dynamics, for example, are driven by the synchronous postwinter activity of larval and nymphal ticks, which is in turn determined by the rapidity of autumnal ground temperature declines (Randolph, Green, Peacey, & Rogers, 2000). Collectively, these abiotic and biotic variables sculpt intricate infection dynamics which are then likely made more complex for generalist parasites that infect multiple hosts. However, empirical studies of the spatial and temporal variation in prevalence of generalist parasites across multiple host species remain rare (Penczykowski et al., 2016;Rigaud et al., 2010).
Avian haemosporidia (encompassing the agents of avian malaria) are a species-rich assemblage of blood parasites which often occur as complex communities within host populations (Harvey & Voelker, 2019;Oakgrove et al., 2014;Ricklefs et al., 2005;van Rooyen, Lalubin, Glaizot, & Christe, 2013b). The most prevalent genera are as follows: Haemoproteus, Leucocytozoon, and Plasmodium ). Our ability to detect these parasites using molecular methods has enabled a better understanding of haemosporidian community ecology and epidemiology (Bensch et al., 2000Hellgren, Waldenström, & Bensch, 2004;Waldenström, Bensch, Hasselquist, & Ostman, 2004), and recent studies have demonstrated the value of using avian haemosporidia to explore fundamental questions in parasitology (Drovetski et al., 2014;Jenkins, Delhaye, & Christe, 2015;Mata, da Silva, Lopes, & Drovetski, 2015). Haemosporidia adopt both host-generalist and host-specialist infection strategies, yet few studies have attempted to resolve how these host-specific strategies might impact parasite ecology Ventim et al., 2012). Climatic variation across an elevational gradient and throughout the year is predicted to impact the spatial and temporal distribution of haemosporidia, primarily via an influence on vector ecology and parasite development (Atkinson et al., 2014;Harrigan et al., 2014;LaPointe et al., 2010). For example, additional to the thermal requirements for P. relictum parasite development (i.e., shown by LaPointe et al., 2010), elevation has also been linked to host infection risk. In Hawaii, mosquito vectors decrease in numbers at high elevations such that populations of honeycreepers living above 1,000 m interact less with mosquitoes and have decreased exposure to Plasmodium (Atkinson & LaPointe, 2009). Plasmodium relictum prevalence in vectors has also been shown to vary seasonally, either in response to within-vector competition between parasites or due to the availability of infected hosts (Lalubin, Delédevant, Glaizot, & Christe, 2013). However, many studies which identify these spatial and temporal patterns in haemosporidian prevalence have focused on either single host or parasite species, or have explored generalizations across the host or parasite community (i.e., by pooling infection data for host or parasite species) (Cosgrove et al., 2008;Huang, Dong, Zhang, & Zhang, 2015;Oakgrove et al., 2014;Pulgarín-R, Gómez, Robinson, Ricklefs, & Cadena, 2018;Santiago-Alarcon, Bloch, Rolshausen, Schaefer, & Segelbacher, 2011).
To examine the extent to which the spatiotemporal infection dynamics of generalist parasites is equivalent among host species, we surveyed haemosporidian parasites in populations of great tits (Parus major) and blue tits (Cyanistes caeruleus) during a three-year period and across an elevational gradient in the French Pyrenees. These two host species are commonly infected by a diversity of haemosporidia, including representatives of the three major genera: Plasmodium, Haemoproteus, and Leucocytozoon, with community structure varying in both time and space (Glaizot et al., 2012;Jenkins & Owens, 2011;la Puente et al., 2010;van Rooyen, Lalubin, Glaizot, & Christe, 2013a;Schumm et al., 2019;Wood et al., 2007). Existing molecular detection methods, which amplify a region of the parasite's mitochondrial genome, are efficient at detecting Plasmodium and Haemoproteus (Hellgren et al., 2004). However, nontarget amplification has been reported to occur in blue tits where Leucocytozoon parasites also prevail (Cosgrove, Day, & Sheldon, 2006;Cosgrove et al., 2008). To address this issue, we designed a new amplification method aimed at specifically detecting Plasmodium and Haemoproteus infections to the exclusion of Leucocytozoon infections. In conjunction with a previous method designed for Leucocytozoon amplification (Hellgren et al., 2004), we surveyed the haemosporidian community with two main aims. First, we report the prevalence and diversity of haemosporidian parasites within our two host species, revealing genetically diverse parasites, which appear to adopt host-specific or host-generalist strategies. Next, while controlling for important infection predictors (i.e., host age), we identified whether the prevalence of generalist parasites varies similarly across season and elevation between the two host species.

| Study populations and sampling
We captured adult great tits (N = 443) and blue tits (N = 531) in the Ariège Pyrenees in France. Four study sites of nest box populations have been established within 14 km of one another close to the Station for Theoretical and Experimental Ecology in Moulis (42°57′29″N, 1°05′12″E) and cover an elevational range from 430 to 1,530 m. As expected, temperatures decrease with increasing elevation at our study sites (Bründl, 2018 N-42°55′07″N, 01°05′40″E-01°03′43″E). With each 1,000 m increase in elevation, we could expect an approximate decrease in surface air temperature of 5.5°-6.5°C (Rolland, 2003). This variation is seen most starkly at temperature minima, for example, in 2016-2017 a difference of 5°C was recorded between the low elevation site Moulis (565 m, min to max = −7.6 to 32.1°C) and the higher elevations of the site Castera (1,502 m, min to max = −12.6 to 30.8°C). The landscape of this study is predominantly mixed deciduous woodland, interspersed by small patches of conifers and open fields used for low-intensity pastoral farming.
Birds were caught and banded at these sites during the breeding (May and June) and nonbreeding (October-March) seasons (Bründl, 2018). During breeding, birds were captured in their nest boxes using spring traps, while mist nets were deployed near artificial feeders in the nonbreeding season. Individuals used in this study were captured between May 2015 and May 2017, and therefore, these data encompass three breeding seasons and two overwintering periods. On capture, morphological data were recorded (sex and age). Sex was determined based on plumage characteristics and the presence or absence of the brood patch. Because they are less sexually dimorphic, the sex of blue tits was not recorded for individuals captured outside of the breeding season. Individuals were categorized into either of two age classes: as yearlings or as postsecond year adults based on contrast between the greater and primary coverts (Svensson, 1992). On first capture, all individuals were fitted with a CRBPO metal identification ring and unique colored band combination. This allowed for the recording of recaptured individuals. A ~35 μl blood sample was obtained by brachial venipuncture.

| Molecular analysis
We extracted DNA from blood samples using the DNeasy Blood & Tissue extraction kit (QIAGEN ® ) following the manufacturer's protocol pertaining to nucleated blood. We standardized these extractions to working concentrations of 25 ng/μl. To detect haemosporidians, samples were subjected to two nested-polymerase chain reaction (PCR) methods which target a specific region of the blood parasites' mitochondrial cytochrome b (Cytb) gene. To amplify the DNA of Leucocytozoon spp., we followed the protocol developed by Hellgren et al. (2004), using primers and cycle conditions shown in Table S1.  (Bensch et al., 2000;Cosgrove et al., 2008;Glaizot et al., 2012;Jenkins & Owens, 2011). These primers were designed to be specific to Plasmodium and Haemoproteus sequences, while still encompassing the region of the apicomplexan's mitochondrial Cytb gene now widely used to classify these parasites (Bensch et al., 2000. As the primer binding regions selected here are conserved between Plasmodium and Haemoproteus species, these primers are suitable for amplifying a range of species within these genera (e.g., between three dissimilar species of haemosporidia: P. cathemerium (AY377128), P. elongatum (AF069611), and H. multipigmentatus (KY653756)) the primers only mismatch the target region by a maximum of 2bp). We adopted this approach as our previous application of Plasmodium-and Haemoproteus-specific primers (as found in Hellgren et al. (2004)) had led to unintentional amplification of Leucocytozoon DNA. This issue has been reported elsewhere (Cosgrove et al., 2006) and is likely most noticeable in our study population due to the very high prevalence rates of Leucocytozoon species (>90%). We concur with Cosgrove et al. (2006)  This mixture was incubated at 37°C for 40 min, heated to 80°C for 10 min, and then cooled to 4°C. Cleaned product was then diluted to nucleic acid concentrations of approximately ≤10 ng/μl, 5 μl of which was sequenced either bidirectionally (using both primers) or unidirectionally using the primers HaemFL (for Leucocytozoon amplicons) or HaemRP (for Plasmodium/Haemoproteus amplicons) (Eurofins sequencing service, Eurofins-MWG).
Sequences were assembled, aligned, and analyzed using Geneious (Geneious ® 9.1.5, Biomatters, Ltd., New Zealand). Parasite lineages were identified by carrying out a BLAST search on the MalAvi database . Sequencing results for Leucocytozoon lineages could at times not be completely resolved due to slightly shorter sequence reads or more commonly due to multiple infections preventing decisive identification. In most cases, however, it was possible to identify the clade or paraphyletic group of lineages to which the lineage or constituent lineages belonged. Multiple infections (infection by more than one lineage of haemosporidian) were identified, as they result in double-peaks on the sequence chromatograms. We identified double-peaks using the Heterozygote Plugin (Geneious ® 9.1.5), marking each contentious peak with the appropriate IUPAC ambiguity code. We then compared these sequences to a library of all previously detected haemosporidian lineage within our population, and through a process of elimination (exclusion of lineages due to the absence of their sequence signals on the target chromatogram), it was possible to determine the most parsimonious combination of lineages which would produce the peak pattern observed (Drovetski et al., 2014). Although it was often not possible to completely resolve these mixed-infections to just two distinct lineages, in most cases we could still identify the presence or absence of the lineage groups common in our populations (e.g., Leucocytozoon lineages clustered into four groups).

| Phylogenetic and statistical analyses
Two phylogenetic analyses were carried out to determine evolutionary relationships and to contextualize the diversity of parasite lineages detected within our populations. In the first phylogeny, all verified haemosporidian lineages from our sampling were included, and sequences trimmed and standardized to 390 bp. The first phy- was generated using a Bayesian reconstruction performed using the MrBayes plugin (Ronquist & Huelsenbeck, 2003) in Geneious and a GTR + I + G substitution model (recommended for our alignment by JModelTest [Darriba, Taboada, Doallo, & Posada, 2012]). Two runs were conducted, both of 10 million generations and with sampling set to every 200 generations. Following a "burn-in" of 25%, the remaining trees were used to calculate posterior probabilities. This approach is used elsewhere in studies of avian haemosporidia (Jenkins & Owens, 2011;Oakgrove et al., 2014). Our second phylogeny was constructed to contextualize two focal Leucocytozoon clades detected in our populations within the genus as a whole. For this phylogeny, P. relictum (this study) was used as outgroup and all sequence data available for Leucocytozoon lineages were downloaded from MalAvi (version 2.4.2, downloaded December 9, 2019). After sequence cropping and removal of partial sequences 447 lineages remained, in addition to two new lineages detected in this study.
Construction was performed using the same parameters as above.
We applied a two-proportion Z test to identify significant differences in parasite prevalence between the two host species. We then tested for significant predictors of infection using six logistic regression models (logit function) with infection status by  (Table 1 for Leucocytozoon groups). We opted not to analyze Leucocytozoon groups B and C, as they were both paraphyletic and some of their constituent lineages showed contrasting distributions within our host species. Of the 972 birds captures, 181 were captured more than once, although the vast majority were only recaptured once. Because few birds have repeat samples and when they do the number of repeats is modally 1, we elected to remove recaptures from the analyses and to perform more robust nonmixed modeling than attempting to control for a lack of statistical independence arising from repeat samples by fitting individual identity as a random intercept (Harrison et al., 2018). Host species, host age, capture season, capture elevation (continuous), capture year, and interaction terms between TA B L E 1 Prevalence of haemosporidian lineages detected in great tits and blue tits in our study population Note: Lineage names were identified using MalAvi . "Count" refers to the total number of infections by that lineage or clade. "As %" provides the prevalence of that lineage or clade in each host species. Two novel lineages are indicated with an asterisk and italicized. For classifications in bold we provide totals including infections that were not identified to lineage, but which could be classified as belonging to that group, species, or genera. Significant z values in bold reflect host-specific differences in parasites prevalence, shading indicates the more common host species; yellow for great tits, blue for blue tits.
host species and elevation, host species and season, and season and year were included as explanatory terms. Terms were dropped if doing so reduced the Akaike information criteria (AIC) estimator by at least 2.0, improving model fit (Zuur, Ieno, Walker, Saveliev, & Smith, 2009 lineage-specific host prevalence's, as such we treated these two clades as distinct parasites in our analyses. We contextualized these two groups within the genus and revealed both that their divergence from one another was substantial (Figure 2). Prevalence for 12 of the 27 lineages and all four of the Leucocytozoon groups were significantly higher in one host, with most lineages (8 of 12) found at higher prevalence in great tits ( lence was more pronounced in blue tits than great tits, as older blue tits were more than twice as likely to be infected than first year birds.

| Elevational and seasonal distribution of infections
For two of the four parasites, capture elevation was retained in the models, although the direction of this effect was parasite specific ( In three of the four parasites there was significant variation in infection prevalence between seasons (Table 2). For H. majoris, spring brought a substantial peak in infections, with prevalence recorded as: spring = 15%, autumn = 3%, and winter = 0.6%. Meanwhile, P. relictum prevalence was reduced in spring (14%) compared to autumn or winter (28%), although this was not determined to be significant ( Figure 4c,d).

| Host-dependent distributions
For most parasites, we found significant host-specific differences in infection distribution relating to either elevation or season. For P. relictum, as well as host-dependent prevalence, we found that the decrease in prevalence at higher elevations was driven primarily by infections in great tits (Figure 3). Host-specific differences in seasonal prevalence patterns were particularly striking; while the prevalence of Leucocytozoon group A in great tits decreased between autumn and spring, it increased in blue tits (Figure 4a). For group D, the opposite was true: prevalence in spring was greater in great tits relative to autumn (60% of group D infections in great tits were detected in the spring), while prevalence was unchanged in blue tits ( Figure 4b).   (Drovetski et al., 2014;Jenkins & Owens, 2011;Mata et al., 2015) and, as reported here, are found almost exclusively in blue tits (Cyanistes spp.). Mixed-infections comprised of multiple lineages from this clade were common (40%) and between-lineage sequence divergences were low (<1.1%, maximum of 4 bp difference between lineages) lending weight to the interpretation that group D may represent within-species variation attributable to a highly host-specific Leucocytozoon species. However, additional morphological and genetic analysis is required to confidently determine the biological relevance of these groupings (Nilsson et al., 2016;Sehgal et al., 2006;Valkiūnas, Sehgal, Iezhova, & Hull, 2010). Contrasting other studies of haemosporidia in these hosts, we found lower prevalence of both P. relictum and H. majoris, which may be driven by local vector abundance (Tomás et al., 2007). Studies conducted at lower altitude sites close to larger water bodies (ideal breeding grounds for mosquitoes; (Glaizot et al., 2012)) find Plasmodium prevalence to be at least twice that reported here for either blue tits (Knowles et al., 2011) or great tits (van Rooyen et al., 2013a).

TA B L E 2 Predictors of infection probability
There is therefore considerable variation in haemosporidia prevalence for these host species across their ranges and additional to this, as shown here and by Schumm et al. (2019), variation also exists between the host species within a shared habitat.
We detected primarily generalist lineages (i.e., lineages detected in the two host species). There was not, however, equivalence in prevalence between host species, as 11 lineages and five pooled classifications (P. relictum and the four Leucocytozoon group) had significantly higher prevalence in one of the hosts. Leucocytozoon parasites showed a host preference that was lineage group (A-D) dependent. As previously mentioned, group D showed high host specificity for blue tits, with more than 20-fold higher prevalence.
Meanwhile Leucocytozoon lineages PARUS16 and PARUS17 in group C were almost exclusively found in great tits, and group A prevalence was three times higher in great tits than in blue tits. These patterns could represent differential exposure to vectors, however to our knowledge, no study has identified host preferences in relevant vectors (e.g., species of Culicidae, Simuliidae, or Ceratopogonidae) extending to such similar host species, nor between host species with such similar ecologies. As such, it seems more likely that host-dependent prevalence of different lineage groups represents host specialization of the parasites per se. An interesting example is Leucocytozoon group A (predominately lineage PARUS22), which has a wide species range; eight host species across five genera ). This group is therefore a broad host generalist, yet we find here host-bias in prevalence, with great tits more typically parasitized than blue tits. This finding is reminiscent of recent studies which have revealed haemosporidian generalists to have host-species specific prevalence and infection intensities within their host ranges Huang, Ellis, Jönsson, & Bensch, 2018).
Perhaps predictably, the haemosporidia we detected displayed distinct elevational and seasonal distributions as both are expected to play a significant role in the probability of infection. Leucocytozoon group A prevalence increased at higher elevations (in both host species) while Leucocytozoon group D prevalence did not change across the elevational gradient (Figure 3b). It seems reasonable to hypothesize that these two Leucocytozoon groups may be vectored  . In contrast to elevation, parasite seasonal prevalence patterns were often specific to the host species infected.
For great tits, Leucocytozoon group A decreased in prevalence between autumn and spring, while in blue tits spring marked the peak in prevalence (Figure 4a). Similarly, the highly blue tit-specific parasite Leucocytozoon group D was detected primarily in the spring for great tits; when prevalence decreased in blue tits (Figure 4b).
Assuming these parasites are each transmitted by the same vector to both host species, these prevalence patterns are contrary to the expected (e.g., peak prevalence mirroring vector activity (Cosgrove et al., 2008)). They suggest instead that host-specific infection dynamics may be driving population level prevalence rates. These dynamics could result from host-specific environmental conditions; such as immune responses linked to seasonally variable conditions (Dowell, 2001), or from interactions occurring between coinfecting and host-specific parasite communities (Lello, Boag, Fenton, Stevenson, & Hudson, 2004;Read & Taylor, 2001). Whatever the driver of these host-specific seasonal prevalence patterns, they will likely be consequential for the trajectories of host-parasite co-evolution in this system.
To conclude, by contrasting the complete haemosporidia community of two geographically overlapping bird species, we have revealed parasite host specificity apparent in both host preference and spatiotemporal prevalence patterns. Most haemosporidia were capable of infecting both bird species, and yet, despite the ecological and taxonomic proximity of the hosts, parasite prevalence was frequently host dependent. What causes these host preferences (i.e., host-parasite dynamics vs. vector-driven) and the role they play in divergence between lineages (e.g., Leucocytozoon clades) could be fruitfully explored using a system such as ours.
Understanding the mechanisms underlying these patterns would benefit from knowledge of vector host preferences, the timing of vector emergence, and the role of parasite-parasite interactions within hosts. Considering the high prevalence of co-infection detected in this study, a future exploration that focuses on both multiple infections, as well as host-specific dynamics, could provide promising insight into the ecology and epidemiology of complex host-parasite communities.

ACK N OWLED G M ENTS
We thank Jessica Lewis, Ethan Hermer, Maya Mould, and all Moulis volunteers for help with capture and data collection. This work was funded by research grants from the Agence National pour la Recherche (ANR-JCJC "NetSelect"), the Human Frontiers Science Partnership (RGP0006/2015 "WildCog") to ASC, and a Royal Society International Exchange grant to CB (IE150476). This work is part of the Laboratoire d'Excellence (LABEX) entitled TULIP (ANR-10-LABX-41). We also thank the two anonymous reviewers whose comments helped improve and clarify this manuscript.

CO N FLI C T O F I NTE R E S T
None declared.

E TH I C A L A PPROVA L
Capture of birds and blood collection was conducted under animal care permits from the French bird ringing office (CRBPO; n°13619; PP576) held by Dr. Alexis Chaine.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in Dryad (https://doi.org/10.5061/dryad.b2rbn zsbf). Novel lineage sequences are available on GenBank (MN782320-MN782321) and have also been uploaded to MalAvi.