Determinants of haemosporidian single and co-infections risks in western paleartic birds

Co-infections with multiple pathogens are common in the wild and may act as a strong selective pressure on both host and parasite evolution. Yet, contrary to single infection, the factors that shape co-infection risk are largely under-investigated. Here, we explored the extent to which bird ecology and phylogeny impact single and co-infection probabilities by haemosporidian parasites using large datasets from museum collections and a Bayesian phylogenetic modelling framework. While both phylogeny and species attributes (e.g. size of the geographic range, life-history strategy, migration) were relevant predictors of co-infection risk, these factors were less pertinent in predicting the probability of being single infected. Our study suggests that co-infection risk is under a stronger deterministic control than single-infection risk. These results underscore the combined influence of host evolutionary history and species attributes in determining single and co-infection pattern providing new avenues regarding our ability to predict infection risk in the wild.

). Sampling was conducted on salvaged birds that were obtained 129 between 1990 and 2019. and consisted of tissues (muscle and liver) stored in 85% ethanol at    Table S1). The first axis explained 62.7% of the variability 156 and represented a gradient going from fast (negative values) to slow (positive values) life- The trophic niche of bird species was estimated using 35 variables describing the diet 159 during the breeding season (Pearman et al. 2014). Specifically, we considered 14 variables 160 characterizing diet (e.g. seeds, fruits, invertebrates, fishes), nine variables characterizing food 161 acquisition behavior (e.g. air pursuit, foliage glean, dig, probe), nine variables characterizing 162 the substrate from which food is taken (e.g. air, water surface, mud, canopy) and three variables 163 characterizing the daily foraging period (Appendix 1, Table S2). As in Pearman et al. (2014), 164 we also included body weight as a surrogate for total energy requirements. These variables were 165 scored as either 0 or 1, with the exception of body weight, which was scored as the average  The remaining traits (nest type and migration status) were extracted from Storchová & 172 Hořák (2018). Nest type was categorized as either "open" or "closed" while migration status 173 was categorized as "sedentary" (species living in the same area in both the breeding and the 174 non-breeding season), "migratory" (species migrates between breeding and non-breeding 175 season) and "facultative migrant" (species makes irregular shifts in breeding and/or 176 nonbreeding season). 178 Estimating species climatic realized niches (i.e. the set of suitable environmental conditions 179 accessible to the species and constrained by biotic interactions; Jackson and Overpeck (2000)) 180 requires data for the full geographical range of species together with the corresponding 181 environmental variables (Guisan et al., 2017). For each species, we estimated its climatic niche 182 by cross-referencing IUCN range maps (https://www.iucnredlist.org/; only considering the 183 resident and the breeding range) with climatic data. Climatic layers for 19 bioclimatic variables 184 were extracted from worldclim (https://www.worldclim.org/) at a 10' resolution (roughly 340 185 km² at the equator). From these variables, we performed a PCA from which we extracted the 186 two first axes, which explained 55% and 19% of the total variance, respectively. Owing to the 187 difficulty to characterize the climatic niche of migratory species, analyses presented below were 188 repeated with full and facultative migrants excluded (see Fig. S1, Fig. S2B).

189
To obtain the climatic niche of each species, we first projected the environmental values 190 corresponding to the geographic range coordinates falling inside the IUCN polygon within the 191 two-dimensional space defined by the two PCA axes. We then used a kernel density estimator 192 (KDE) to delineate species envelopes (Fig. 1) probability density that included 99% of points (to leave out environmentally atypical 198 occurrences). From species envelopes, we extracted its area (niche breadth) and computed its 199 centroid as the mean of point coordinates falling inside the delimited niche. We then extracted 200 the coordinates of the centroid on each of the two PCA axes (Fig. 1). To test the robustness of 201 our results we used two other algorithms to delineate species realized niches: convex hulls and 202 alpha hulls (see Fig. S2, Fig. S3). The area of the geographic range was directly extracted from 203 IUCN polygons.

204
Infection prevalence among avian phylogenetic tree 205 To illustrate the effects of host evolutionary history on single and co-infection probabilities, we 206 followed the methods of Barrow et al. (2019). Specifically, we used 1000 trees from 207 BirdTree.org (backbone tree from Hackett et al. 2008) to generate a consensus tree using the 208 ls.consensus function (package phytools v. 0.6-00, Revell 2020). The prevalence of infections 209 (i.e. the proportion of birds infected by haemosporidian parasites altogether and by each 210 haemosporidian genus separately) was calculated for each species and mapped on the 211 phylogenetic tree using the contMap function (package phytools). The prevalence of infection 212 was visualized over the whole phylogeny using species for which there were at least 5 samples. 214 Overall, we computed ten predictors: the first PCA axis summarizing bird life-history strategies 215 (i.e. position along the slow-fast life history continuum), the two OA axes summarizing bird 216 trophic niches, the climatic niche breadth, the geographic range size, the coordinates of the 217 climatic niche centroid on each PCA axis, the migration status, the nest type and the museum 218 affiliation (to account for potential effect related to tissue collection and/or conservation). To 219 prevent collinearity between predictors, we checked that variables had Pearson's correlation 220 coefficients ρ ≤ |0.7| (Dormann et al. 2013). We found that OA1 and the geographic range size 221 were strongly correlated with bird life-history strategies (ρ = -0.81) and the climatic niche 222 breadth (ρ = 0.76), respectively. These variables were thus removed from the analysis.

223
However, we also tested the effect of these variables in separate models where they replaced 224 the variables they were correlated with. 225 We built two types of phylogenetic multilevel models using Bayesian multilevel models 226 from the brms package (Bürkner et al. 2021). First, we used a phylogenetic multinomial model 227 with a "Categorical" error distribution to study the effects of species attributes and phylogeny 228 on the probability to be single and co-infected. Second, we used three phylogenetic generalized 229 linear multilevel models with a "Bernoulli" error distribution to study the effects of species 230 attributes and phylogeny on the probability to be infected by each parasite genus, separately.

231
For all models, host species and phylogeny were treated as random effects. Continuous 232 predictors were standardized to z-scores (mean=0, variance=1) to improve model convergence 233 and parameter estimation. We used the default priors of the brms package (i.e. weakly 234 informative priors) and ran three chains with 11,000 iterations. The first 1,000 iterations were 235 considered as burn-in and were thus discarded. Chains were thinned every 10 iterations. To 236 account for phylogenetic uncertainty, we randomly sampled 100 trees from the set of trees 237 extracted from BirdTree.org and ran the analyses on each tree. We then combined all models 238 using the combine_models function. Note that this approach critically depends on the 239 assumption that all trees are equally likely. We followed the methods of Barrow et al. (2019)

245
Unless otherwise stated, all results were consistent regardless of the method used to delineate 246 species envelopes as well as when migratory species were removed.

247
Single and co-infection prevalence 248 We detected 486 infected individuals (35.7%), among which 101 (7.4%) were co-infected by 249 at least two different parasite genera (20.8% of infected birds). Focusing on species for which 250 there were at least 5 samples, we observed a large variation in the prevalence of infections 251 across the phylogeny (Fig. 2). Species with the highest prevalence (>90% of infected 252 individuals) included Asio otus, Falco subbuteo and three passerines: Corvus corone, 253 Coccothraustes coccothraustes and Turdus philomelos (Fig. 2). In contrast, eleven species 254 belonging to eight different genera were never infected (Fig. 2). The rate of co-infection greatly 255 varied between species, ranging from zero individual in 32 species to more than 50% 256 individuals in seven species (Fig. 2). The most frequent parasite genus was Leucocytozoon 257 (51.03% of infected birds) followed by Plasmodium (38.65%) and Haemoproteus (26.38%).

258
The infection prevalence of each haemosporidian genus also varied greatly across the 259 phylogeny ( Fig. S4) with some families exhibiting higher infection prevalence than others. For 260 instance, Leucocytozoon infection prevalence was relatively high for Corvidae (78%) whereas 261 Plasmodium infection prevalence was the highest for Turdidae (86%).

262
Predictors of haemosporidian single and co-infection 263 The structure (trophic niche, OA1) and the height of the foraging environment (OA2), the area 264 of the climatic niche (niche breadth) and the location of the niche centroid were not pertinent 265 to explain the probability to be single or co-infected (i.e. the 95% confidence interval [CI] 266 overlapped with 0, Fig. 3). We note however that while this absence of effect was generally 267 robust to the type of algorithm used to delineate climatic niches (Fig. S2) and the removal of 268 migratory species, we detected a tendency toward an effect of niche breadth on co-infection 269 probability when excluding migratory species (Fig. S1, Fig. S2)). Migratory species (full and

283
None of the tested variables appeared to be a strong predictor of Plasmodium infection 284 probability (Fig. 4, Fig. S3). On the other hand, Haemoproteus infection probability tended to 285 be impacted by both the height of the foraging environment and species life-history strategies 286 thus indicating a lower probability of infection for slow-living species and species foraging on 287 the ground (Fig. 4). Species with open nests also tended to have higher Haemoproteus 288 prevalence than species with closed nests (Fig. 4). Regarding Leucocytozoon, infection 289 probability tended to increase with the geographic range size (Fig. 4) but to be lower for 290 facultative migrants and species with closed nests (Fig. 4)

294
In this study, we evaluated the extent to which host phylogeny and species attributes related to 295 climatic niche properties and other ecological and life-history traits influenced the probability 296 of haemosporidian infection and co-infection in western Palearctic birds. We found that (1) 297 while some attributes (e.g. migration) influenced both the probability of being single and co-298 infected, others (e.g. geographic range size) were only pertinent in explaining variation in co-299 infection probability and (2) that phylogeny is a far more important predictor of co-infection 300 probability than of single infection probability. Overall, the effect size of all predictors and the 301 proportion of variance explained by our models were significantly larger for co-infection than 302 for single infection probability. Altogether, these findings suggest that co-infection probability 303 is under a stronger deterministic control than single infection which may rather be influenced 304 by stochastic processes (e.g. random encounter rate with parasites) or by other factors not 305 included in our study (e.g. hosts' features).

306
Two parameters associated with host ecology clearly influenced single and co-infection 307 probabilities. First, we found an effect of nesting behavior, with the probability to be single and

387
Beyond the implications of this study for predicting host infection status, it demonstrates 388 that tissue sample collections from museum specimens can be used to investigate evolutionary 389 and ecological hypotheses at a relatively low cost and with no impact on wild populations.

390
However, this also implies certain limitations. For example, our samples come from dead or 391 injured individuals and our results may be biased towards weak individuals, who may have less 392 resistance to pathogens (particularly for individuals found dead without injuries). In addition, 393 we focused on haemosporidian parasites and we cannot rule out the possibility that birds have 394 been infected by other blood parasites genus (e.g. Trypanosoma, Babesia) or by macroparasites 395 such as helminths (e.g. Trematodes, Nematodes, Cestodes). We also concentrated on co- By taking advantage of a large set of samples associated with museum specimens and using a 405 state-of-the-art bayesian phylogenetic modelling framework, we identified relevant predictors 406 of among-species differences in haemosporidian (co-) infection risks for western palearctic 407 birds. Interestingly, we showed that our ability to predict co-infection risk was much higher 408 than single infection risk, suggesting that random processes (e.g. encounter rate) may be more 409 prevalent for the latter than for the former, where deterministic processes (e.g. species and parasites will not only help us understand why some species are more susceptible to (co-419 )infection than other, but will also be of importance to address urgent public health problems 420 regarding the emergence and evolution of infectious diseases. Ecol. Resour., 9, 1353Resour., 9, -1358 Blonder           Loxia curvirostra