Prevalence and diversity of haemosporidian parasites in the yellow‐rumped warbler hybrid zone

Abstract Parasites can play a role in speciation, by exerting different selection pressures on different host lineages, leading to reproductive barriers in regions of possible interbreeding. Hybrid zones therefore offer an ideal system to study the effect of parasites on speciation. Here, we study a hybrid zone in the foothills of the Rocky Mountains where two yellow‐rumped warbler subspecies, Setophaga coronata coronata and S. c. auduboni, interbreed. There is partial reproductive isolation between them, but no evidence of strong assortative mating within the hybrid zone, suggesting the existence of a postzygotic selection against hybrids. Here, we test whether haemosporidian parasites might play a role in selecting against hybrids between S. c. coronata and S. c. auduboni. We screened birds from five transects across the hybrid zone for three phylogenetic groupings of avian haemosporidians Plasmodium, Haemoproteus and Leucocytozoon parasites and quantified intensity of infection. Contrary to our prediction, hybrids did not have higher haemosporidian parasite prevalence. Variation in Haemoproteus prevalence was best explained by an interaction between a birds’ hybrid index and elevation, while the probability of infection with Leucocytozoon parasites was only influenced by elevation. We also found no significant difference in the diversity of haemosporidian lineages between the warbler subspecies and their hybrids. Finally, intensity of infection by Haemoproteus increased significantly with elevation, but was not significantly linked to birds’ hybrid index. In conclusion, our data suggest that haemosporidian parasites do not seem to play a major role in selecting against hybrids in this system.

To investigate parasite-mediated speciation, tension zones, a particular type of hybrid zone, are relevant. Tension zones are characterized by a relatively narrow width and stabilized by a balance between selection against the hybrids within the zone and dispersal of parental types into the zone (Barton & Hewitt, 1985).
The sources of selection against hybrids in most tension zones are unknown, although the role of parasites in this context has not been extensively studied (Alexandrino et al., 2005;Singhal & Moritz, 2012).
One well-characterized hybrid zone where hybrid fitness is unclear is between myrtle warblers (Setophaga coronata coronata) and Audubon's warblers (S. c. auduboni) (Figure 1). These two subspecies interbreed in a narrow region across the Canadian Rocky Mountains (Hubbard, 1969). Brelsford and Irwin (2009) tested birds in the hybrid zone for evidence of assortative mating, and found a pairing pattern consistent with random or very weak assortative mating. Brelsford and Irwin (2009) also found that selection against hybrids is necessary to maintain the observed linkage disequilibrium and cline width.
The mechanisms of this inferred selection are still unclear. Here, we test the hypothesis that haemosporidian parasites may play a role in the selection against S. coronata × auduboni hybrids.
Patterns of prevalence and diversity of haemosporidian parasites can vary according to several different biotic and abiotic factors, yet the underlying factors are still poorly understood. Scordato and Kardish (2014) showed that host species is one of the main predictors of prevalence and diversity, and a better predictor than geography.
Here, we investigated the potential role of haemosporidian parasites in selecting against S. c. coronata × auduboni hybrids.
The data set we used is extensive in sample size and geographic area, composed of five independent transects, which allowed us to test hypotheses regarding the geographical structuring of haemosporidian distribution. By amplifying and sequencing a fragment of haemosporidian cytochrome b in blood samples from warblers in the hybrid zone, we assessed the infection status and, when possible, identified the parasite lineages present. Using a quantitative PCR protocol, we measured the intensity of infection (parasitaemia) by the most common lineage in the hybrid zone.
By these means, we determined whether hybrids had higher haemosporidian prevalence and parasitaemia than pure S. c. coronata and S. c. auduboni, which could be expected if hybrids were less adapted to resist haemosporidian infection. As hybrids may inherit specialist parasites from both parental subspecies, we also tested for higher haemosporidian lineage diversity in hybrids.
However, differences in infections may be affected by environmental parameters as well as bird ancestry, so we also assessed how geographical variables, specifically elevation, influenced haemosporidian prevalence in myrtle and Audubon's warblers.

| Field sampling
The majority of samples (n = 196 S. c. coronata, n = 193 S. c. auduboni and n = 228 hybrids) used for this study were initially collected by Brelsford and Irwin (2009). Briefly, we captured warblers defending their territories on their breeding grounds (that is, mostly males), along five transects across the hybrid zone and additional allopatric sites, in Alberta and British Columbia ( Figure 2). We took several morphometric measurements and scored five plumage color traits that differ between S. c. coronata and S. c. auduboni, following Hubbard's (1969) hybrid index: 0 for a coronata-like trait, 2 for an auduboni-like trait, and 1 for intermediate states. The hybrid index is the mean of the scores of the different traits, listed in Supporting Information Table S1. We excluded the sixth trait used by Hubbard (1969), tail pattern, from analysis due to concerns over its repeatability. We also F I G U R E 1 Migrating yellow-rumped warblers (photograph credit: David P. L. Toews) F I G U R E 2 Map of sampling sites (black dots). Letters indicate the five transects excluded auricular color from analysis in females, because in both subspecies female auricular patches were brown rather than black or gray.
In addition, we also sampled yellow-rumped warblers during autumn migration in 2015 (n = 131), between August 28th and September 12th near Barrier Lake (51.023591, −115.060657, elevation: ca. 1,400 m) in the Kananaskis region of Alberta, Canada. The aim of this new sampling was to catch hatch-year birds starting their first migration to compare their parasite lineages composition to that of birds that have already been on wintering ground. This site is close to the geographic center of the hybrid zone, but birds captured during migration may have come from distant breeding locations.
We captured migrating yellow-rumped warblers using mist nets with song and call playback, and determined the age of the birds (hatchyear/after hatch-year) by examining skull ossification (Norris, 1961).
We also took morphometric measurements. The use of the full hybrid index (Hubbard, 1969) is more difficult when the birds do not display their breeding plumage, so in this case, we used a genomic hybrid index (see below). Blood samples were obtained by brachial venipuncture and stored in Queen's lysis buffer (Seutin, White, & Boag, 1991).

| Genomic hybrid index
We extracted genomic DNA from each sample either using a standard phenol-chloroform extraction or with a Qiagen ® DNeasy kits (Qiagen, Valencia, CA) and determined the final concentration of each extraction using Qubit Fluorometric Calibration (QFC; Invitrogen, Carlsbad, CA). To estimate a high-resolution genomic hybrid index, we used a double digest restriction association DNA sequencing (ddRAD) protocol following Peterson, Weber, Kay, Fisher, and Hoekstra (2012) with the modifications outlined in Campagna, Gronau, Silveira, Siepel, and Lovette (2015). We sequenced the two lanes of an Illumina HiSeq 2000 (150 bp, singleend) at the Cornell University Life Sciences Core Laboratories Center (Ithaca, NY).
We demultiplexed sequencing reads within each index group using the barcode-splitting program Sabre (https://github.com/najoshi/sabre), allowing for one mismatch in the barcode plus enzyme cut-site sequence. We used BOWTIE2 (Langmead & Salzberg, 2012) to map each of the individual reads to a build of the myrtle warbler genome (Toews, Brelsford, Grossen, Milá, & Irwin, 2016), using the "very sensitive local" set of alignment presets. For SNP discovery and variant calling, we used the UnifiedGenotyper in GATK (DePristo et al., 2011), and used GATK and VCFtools (Danecek et al., 2011) to apply the quality filters outlined in Toews, Taylor et al. (2016). We coded genotypes with a Phred-scaled quality lower than 20 as missing data and excluded loci with more than 30% missing data and/or a minor allele frequency of less than 1%.
To estimate hybrid ancestry, we used the program STRUCTURE (version 2.3.4, Falush, Stephens, & Pritchard, 2003). For this analysis, we focused on markers that were spaced at least 10 kb apart using the "thin" function in VCFTools (Danecek et al., 2011). We then ran STRUCTURE on this subset of SNPs (n = 4,661 loci) with K = 2 for 100,000 MCMC steps following a burn-in of 100,000 iterations.

| Determination of haemosporidian parasites
Infections by haemosporidian parasites were diagnosed by performing nested PCR as described in Jenkins, Delhaye, and Christe (2015), modified from Hellgren, Waldenström, and Bensch (2004). Briefly, the first PCR round was conducted using HaemNF1 and HaemNR3 primer pair in order to amplify a 617 bp conserved region of the haemosporidian (Plasmodium, Haemoproteus, and Leucocytozoon) cytochrome b (cytb) gene. We then amplified a 479 bp region from 1 μl of the product of the first PCR round, using the Plasmodiumand Haemoproteus-specific HaemF and HaemR2 primers pair and the Leucocytozoon-specific HaemFL and HaemR2L primers pair.
PCR products were then visualized after agarose gel electrophoresis (2%). Positive PCR products were sequenced in both directions by Sanger sequencing (Microsynth AG, Balgach, Switzerland).
Double-peak(s) in a chromatogram meant that several haemosporidian lineages were present in the sample ); these mixed infections were excluded from diversity analyses.
Sequences were blasted against the GenBank ® and MalAvi databases . Sequences that did not match at 100% of identity with any deposited sequence were named SETCOR03 to SETCOR19.

| Sequence analyses and phylogenetic reconstruction
In order to compute phylogenetic diversity metrics, we reconstructed a phylogeny of the sampled haemosporidian lineages.

| Parasitaemia measurements
We performed quantitative PCR on samples infected by hDEN-PEN02, the most prevalent haemosporidian lineage in our samples, to measure relative parasitaemia. In order to do so, we Host and parasite DNA concentrations (α) were calculated as: C t being the mean of the measured C t , I the intercept of the standard curve, and m the slope of the standard curve. Relative parasitaemia (R) was calculated as the ratio between parasite and host DNA: R was log-transformed in order to normalize the distribution.

| Analysis of the probability and intensity of infection
We used generalized linear mixed models (GLMMs) to examine the influence of hybrid index, elevation (as a proxy for associated environmental variables), and scaled mass index (SMI, as a proxy for bird body condition; Peig & Green, 2009) on probability of infection and coinfection in the hybrid zone, using site and transect of sampling as random factors, with the function "glmer" implemented in the R package "lme4" (Bates, Mächler, Bolker, & Walker, 2015).
We fitted different models for five different response variables (1) Plasmodium infections, (2) (1) and (2) We used linear mixed models to examine the variation in hDEN-PEN02 parasitaemia in the hybrid zone dependent on host hybrid index, elevation, scaled mass index, and using sites and transects of sampling as random factors (model 6), with the function "lmer" implemented in the R package "lme4" (Bates et al., 2015). We selected the best model according to the same procedure of the analyses of the infection probability.
As data used in models 1-6 were sampled along transects, we computed Moran's I (Gittleman & Kot, 1990) to test for spatial autocorrelation in the residuals of these models with the function "Moran.I" implemented in the R package "ape" (Paradis et al., 2004).
To assess the relative importance of age, sex, SMI, and hybrid index on infection probability in the migrating birds sampled in 2015, we fitted generalized linear models (GLMs) using these explanatory variables and the absence/presence of infection as a response variable with a binomial error structure, with the function "glm" implemented in the R package "lme4" (Bates et al., 2015). We tested the interactions between age and sex, as well as the hybrid index as a quadratic function, and selected the best model according to the same procedure as above. We fitted different models for (7) Plasmodium infections, (8) Haemoproteus infections, (9) Leucocytozoon infections, and (10) coinfection by Plasmodium and/ or Haemoproteus and Leucocytozoon. In models (7) and (8)

| Diversity analyses
Analyses regarding diversity of haemosporidian lineages were conducted on the whole sampling region and considering Plasmodium and Haemoproteus infections together, because the number of Plasmodium infections was low, and because Plasmodium and Haemoproteus groups cluster together in the phylogeny; Leucocytozoon infections were analyzed separately. We tested our hypothesis that genetically similar birds had more similar parasite lineages. In order to do this, we conducted a Mantel test using the "mantel.rtest" function implemented in the R package "ade4" (Dray  Dufour, 2007) with 10,000 permutations. We also did a partial Mantel test, correcting for geographic distances, using the "mantel.
partial" function implemented in the R package "vegan" (Oksanen et al., 2015) with the Pearson's method.
To compare the diversity of haemosporidian lineages between subspecies of yellow-rumped warblers and the hybrids, we calculated the Simpson's diversity index (D; Simpson, 1949)   Leucocytozoon infections could not be identified due to poor-quality sequences and were therefore excluded from further analysis.

| Analysis of the probability and intensity of infection
The results of the models to test elevation and hybrid index on probability of infection in the hybrid zone are reported in Table 1.
We found that the probability of infection by Haemoproteus varies according to an interaction between elevation and hybrid index (model 2: 2 1 = 6.36, p = 0.012): The probability of infection decreases with elevation in S. c. auduboni and hybrids, S. c. coronata are more likely to be infected at higher elevation ( Figure 3). This effect seems mainly driven by hDENPEN02 (model 3: 2 1 = 5.23, p = 0.022), a Haemoproteus lineage that was responsible for 96.1% (n = 100) of single Haemoproteus infections in the hybrid zone. Body condition also correlates the probability of infection by Haemoproteus  Figure S1).
A quadratic function of elevation also best predicted parasitaemia of birds infected by hDENPEN02 (model 6: 2 1 = 4.05, p = 0.044), and hybridization status did not influence it, nor the body condition (Table 1 and Figure 4). We found no spatial autocorrelation in the residuals of the models (Supporting Information Table S2).
In birds sampled in autumn 2015 during autumn migration, neither age, sex, nor body condition influenced the probability of infection. However, genomic hybrid index was significantly associated with the probability of infection by Haemoproteus: S. c. coronata was more infected than S. c. auduboni and hybrids (model 8: 2 1 = 5.52, p = 0.010; Table 1, Figure 5). Haemoproteus infections, 95.6% were hDENPEN02 and among identified Leucocytozoon infections, 52.9% were lCNEORN01, and 20.9% were lCB1, which match morphospecies Leucocytozoon majoris. We found 17 haplotypes that did not match at 100% of identity the sequences referenced in GenBank nor in MalAvi database of avian malaria lineages ). In total, eight lineages of the haemosporidian parasites that we found are shared between the two subspecies and the hybrids. One Plasmodium lineage, three

| D ISCUSS I ON
In this study, we broadly sampled a naturally occurring hybrid zone between two yellow-rumped warbler subspecies to explore the role of parasites in potentially selecting against hybrids. We found that haemosporidian parasites-in diversity and prevalence-are unlikely to play a major role in selecting against S. Notes. Generalized linear mixed models (1)-(5) and generalized linear models (7)-(10) were fitted with a binomial error structure (logit link function).
Response variable values were 1 if infected and 0 if uninfected in (1)-(4) and (7)-(9), 1 if coinfected and 0 if single or uninfected in (5) and (10), and the logarithm of parasitaemia in (6). Models (1)-(6) had the transect and site of sampling as random factors. We performed likelihood ratio tests to remove nonsignificant variables in the following order: (a) interaction; and (b) quadratic terms, both kept if at least one of them was significant. If not significant, we report values obtained after removal of this term from the last model that contained it. Explanatory variables in bold are those that were kept in the final model. In every case, elevation, hybrid index and scaled mass index (SMI) were kept in the final model as our data are nonorthogonal, as suggested by Crawley (2012). *p-value <0.05; **p-value <0.01.
TA B L E 1 (Continued) to be more infected by haemosporidian parasites than parental subspecies. We found that prevalence, coinfection probability, parasitaemia and diversity of haemosporidian lineages were not higher in hybrids, as we had originally predicted. In migrating yellow-rumped warblers sampled in autumn, there was no effect of hybrid index on the probability of infection by Plasmodium or Leucocytozoon, but we found that S. c. coronata had a higher probability to be infected by Haemoproteus than S. c. auduboni and hybrids.
From our results, it seems that most haemosporidian lineages were shared between myrtle and Audubon's warblers. Indeed, three lineages were very abundant in S. c. coronata, S. c. auduboni, and their hybrids. Some other lineages were found only in one subspecies, or only in hybrids, or in one subspecies and in hybrids, but this was generally restricted to only one or two individuals. Additional sampling would be required to determine whether these lineages, when undocumented, are rare and specific, or if they are simply rare and exclusivity is due to stochastic sampling effects. According to theory on the evolution of specialization, we could expect specialist haemosporidian parasites to show a higher prevalence and virulence than generalists (Futuyma & Moreno, 1988). This is sometimes the case for virulence: For example, Garamszegi (2006) showed that specialist malaria parasites of primates had a higher parasitaemia (used as a proxy for virulence) than generalists. However, Hellgren, Pérez-Tris, and Bensch (2009) showed that overall more generalist Plasmodium and Haemoproteus parasites also were the most abundant in single subspecies. In addition, the haemosporidian lineages found in the sampled yellow-rumped warblers are relatively generalist and have been found in other species in the MalAvi database . For example, hDENPEN02 has been found in at least seven other species of Passeriformes in North America and lCNEORN01 in ten (Oakgrove et al., 2014;Outlaw & Ricklefs, 2009;Ricklefs & Fallon, 2002;Walther et al., 2016). Given these results, we suggest that it is unlikely that these parasites specialize and thereby exert differential levels of selection between myrtle and Audubon's warblers. In addition, the common lineages we found in the hybrid zone were present in both after hatch-year and hatch-year birds, which had not completed a full annual migration cycle yet. This means these young birds were infected on the breeding ground, as opposed to their wintering grounds, suggesting that infection occurs in the nest or soon after fledging. This is supported by a recent study on Setophaga coronata F I G U R E 4 Parasitaemia in birds infected by hDENPEN02 as a function of elevation (m) F I G U R E 3 Probability of infection by Haemoproteus in relation to elevation (in m) for each yellow-rumped warbler group. For the sake of clarity, we plot values predicted by a model that uses hybrid index of yellow-rumped warblers as categories instead of the continuous variable "hybrid index" presented in the results section. Except for this difference, the model has the same structure to model (2) in Table 1 auduboni in New Mexico (US) sky islands that found a high diversity of haemosporidian lineages but no bird infected by hDENPEN02 (Williamson et al., 2018).
In sampling of yellow-rumped warblers along five transects across the Rocky Mountains, we documented a consistent pattern: Elevation was an important predictor of prevalence, especially with respect to Leucocytozoon infections. Indeed, the probability of infection by Leucocytozoon decreased with elevation in pure individuals and hybrids.
In the context of studying factors shaping haemosporidian parasites distribution in general, especially regarding the importance of predicting biodiversity and shifts in community structures under climatic change, it would be of major value to determine whether this pattern is a result of an effect of elevation, geographic position and their correlated abiotic factors (e.g., temperature, solar radiation, humidity, snow cover) or an effect of correlated biotic factors (e.g., change in plant and animal communities). Future studies on elevational gradients of parasite distribution should test whether these effects specifically act on bird susceptibility, parasite distribution, or on vector distribution and preferences. We propose that the presence of suitable conditions for vectors was the main driver for the observed Leucocytozoon distribution.
Black flies (Simuliidae), the vectors of Leucocytozoon parasites, depend on running water bodies for their reproduction, and so streams and rivers are critical to their distribution (Crosskey, 1990). Temperature is also a factor that influences larval survival (Ross & Merritt, 1978). Many  This mechanism has been supported by the effect of food supplementation on Plasmodium parasitaemia in experimentally infected canaries (Cornet, Bichet, Larcombe, Faivre, & Sorci, 2014) and this effect might well be enhanced in the wild in interaction with other sources of stress.
This study also revealed a pattern, previously identified in great tits (Ots & Horak, 1998 but see Bennett, Caines, & Bishop, 1988), that birds infected by Haemoproteus have a higher body condition than uninfected birds. Indeed, we expected infected birds to have a lower body condition if haemosporidian infections were detrimental for their health. However, this pattern was observed only in the breeding birds sampled in spring, while in birds sampled in autumn 2015 at the beginning of their migration, there was no effect of the body condition on the probability of infection. It is possible that haemosporidian infections are more detrimental to yellow-rumped warblers of poor condition, and that only individuals of good condition are able to survive the winter with an infection. We expect the weakest birds to die especially during the winter because migration is potentially stressful and requires substantial energy reserves (Wikelski et al., 2003). In addition, birds might be exposed to other parasite species on their wintering ground, which would be another source of physiological stress. Evidence of the impact of haemosporidian infections on yellow-rumped warbler health and mortality would be useful to assess the validity of this hypothesis. It would also help us to conclude on the potential role of haemosporidian parasites in selection against S. c. coronata and S. c. auduboni hybrids, as even if hybrids were not more infected, they could suffer higher costs of being infected than their parental subspecies.
In conclusion, we found that haemosporidian parasites seem unlikely to play a major role in imposing a stronger selection pressure in hybrids of yellow-rumped warblers. S. c. coronata × S. c. auduboni hybrids seem to exhibit similar patterns to S. c. auduboni with regard to how elevation affects their infection probability. We also found that both subspecies and hybrids share most of their haemosporidian lineages, which is consistent with these lineages being generalists. Finally, it seems that elevation, or other correlated

DATA ACCE SS I B I LIT Y
Upon acceptance of this manuscript, all data supporting this study will be made available on Dryad. All sequences will be made available on GenBank and MalAvi.