Revealing the drivers of parasite community assembly: using avian haemosporidians to model global dynamics of parasite species turnover

non-migratory hosts). Our findings provide a robust foundation for the prediction of avian pathogen spread and the emergence of infectious diseases.


Introduction
Some species are globally widespread and occur across multiple ecological communities, whereas most biodiversity is geographically restricted, but what determines species turnover among distinct regions?Niche theory predicts species inhabit a subset of conditions and resources they can tolerate and/or in which they thrive.According to Hutchinson (1957), a species' niche is a multidimensional space in which all environmental conditions and resources required for its persistence are met.For parasites (i.e.organisms that exploit another to complete their life cycle), however, resource limitation is mainly related to the presence of suitable hosts and, in some cases, vectors (i.e.mobile blood-feeding invertebrates involved in the transmission of pathogens (Wilson et al. 2017, Mestre et al. 2020)).Indeed, host specificity (i.e.spectrum of hosts a parasite can infect) has been directly linked with their geographic and environmental ranges (Krasnov et al. 2005, de Angeli Dutra et al. 2021a), indicating generalist parasites are more likely to expand their range and disperse globally.Moreover, host specificity can vary across distinct climatic conditions (Fecchio et al. 2019a), suggesting parasite functional traits might shift according to habitat features.As a result, similar host communities living in similar environmental conditions should harbor similar parasite community composition (Fecchio et al. 2019b).
Nonetheless, parasite colonization and assembly are constrained by limitations other than just host species community composition.For instance, host community and parasite dissimilarity tend to increase with increasing geographic distance (Nekola and White 1999, Tuomisto et al. 2003, Dallas and Poisot 2018, Martins et al. 2021), due to dispersal restrictions of organisms and the presence of geographic barriers (Nekola and White 1999).In the case of parasites, the ability to infect highly mobile hosts (e.g.migratory species) could expand their distribution by increasing their dispersal ability and opportunities to invade new host communities and environments (de Angeli Dutra et al. 2021b, Poulin andde Angeli Dutra 2021).However, at the same time, this dispersal mechanism comes with a challenge: adaptation to new environmental conditions and resources (Poulin and de Angeli Dutra 2021).Hence, even if parasites reach a new region via host dispersal, parasite colonization could be prevented by inadequate climatic conditions and/or absence of suitable new hosts or vectors.For this reason, multiple factors can contribute to parasite community composition similarity or dissimilarity among distinct regions/communities across the globe.
Parasites comprise many of the most successful and prevalent organisms in nature, accounting for a third of the world's biodiversity (Clayton et al. 2015, Mestre et al. 2020).Nonetheless, due to the nature of parasite interactions with their hosts, the drivers of parasite community composition similarity and turnover are very particular (Mestre et al. 2020).Parasites can have major impacts on community structure and energy flow through direct and indirect cascading effects on host fitness, population dynamics, ecosystem functions and extinction risk (Dunne 2002, Kuris et al. 2008, Lafferty et al. 2008).For these reasons, elucidating the mechanisms underpinning parasite community composition similarity and turnover among regions is important to predict parasite species distribution worldwide, estimate parasite impact on regional host communities and forecast distributional shifts in response to anthropogenic changes.For instance, limited parasite turnover could indicate that different regions are prone to colonization by multiple distinct parasite species, whereas high turnover might reflect local parasite specialization.Therefore, we aimed to identify the drivers of parasite community composition similarity and turnover worldwide, using avian malaria and malaria-like parasites as a model system.
Avian malaria and malaria-like (haemosporidian) parasites are cosmopolitan vector-borne parasites which are among the most prevalent and diverse wildlife parasites, comprising almost 4000 lineages (Bensch et al. 2009).Avian haemosporidians are mainly represented by three distinct genera: Plasmodium, Haemoproteus and Leucocytozoon (Valkiūnas 2005).These parasites have been intensively studied in the past two decades and represent an ideal system to investigate the mechanisms shaping parasite turnover due to the comprehensive data freely available online on parasite distribution and diversity (Bensch et al. 2009) as well as host distribution and functional traits (Wilman et al. 2014, Dufour et al. 2020).Here, we test how environmental conditions (i.e.temperature and precipitation patterns) and host species distribution and functional traits (e.g.range size, territoriality, migratory status and body mass) relate to global parasite community dissimilarity while also controlling for geographic and host phylogenetic distances.We predict that 1) parasite species turnover increases as a function of increasing dissimilarity among communities (i.e.greater variation in environmental and host features) and 2) host species distribution is a key driver of parasite community composition.Our results identify the main factors shaping avian haemosporidian dissimilarity and turnover worldwide and quantify their relative importance.

Datasets
We obtained datasets from multiple openly available sources.Haemosporidian (i.e.Plasmodium, Haemoproteus and Leucocytozoon) data were extracted from MalAvi (http://130.235.244.92/Malavi/)(Bensch et al. 2009) using the function extract_table from the 'malaviR' package in R (www.r-project.org).MalAvi compiles records of haemosporidian parasites for each site sampled; sites with fewer than 10 records were excluded from the analyses.This first dataset contains information on the haemosporidian lineage identity, infected host species and site (with geographic coordinates) where each parasite lineage was recorded (Fig. 1).Bird distribution data were acquired from BirdLife International (www.birdlife.org/)(BirdLife International and Handbook of the Birds of the World (2020)) and bird functional traits (i.e.body mass, range size and territoriality) were obtained from Open Traits datasets (https://opentraits.org/datasets.html)(Wilman et al. 2014).We used data published by Dufour et al. (2020) to determine bird migratory category.The authors classify migratory birds in three categories depending on the distance the species migrate (short, variable or long) or as resident species.Multiple datasets containing host functional traits data were used to allow us to explore a wider range of functional traits in our analyses.Lastly, environmental data were extracted from Wordclim (https://worldclim.org/)using the function getData from the 'raster' package in R and resolution equal to 10 km.This last dataset is composed of 19 distinct environmental features relating to temperature and precipitation measures.We classified each locality where haemosporidians were recorded based on their geographic coordinates by clustering localities into geographic cell grids of 5 × 5 degrees.We then considered those geographic cell grids as distinct geographic units, each characterized by the occurrence of particular haemosporidian lineages, bird species and their traits and environmental conditions.

Data preparation
Initially, in order to combine all three datasets, we used principal component analyses (PCA) to reduce the dimensions of our data.We performed two PCAs: one for environmental data using all 19 variables extracted from Worldclim and a second PCA for all the functional traits.We then extracted the individual values for the first axis (53% and 72% of variance explained, respectively), which were each incorporated into the bird distribution raster file as a new layer.In order to transform all functional traits variables into a numeric format, we used the function dummy.data.framefrom the 'dummies' package.This function converts a categorical variable and returns a numerical matrix.Then, we transformed each PCA into a raster format and combined them into a unique raster file.Then, we employed this new three-layer raster object as our predictor data in a unique full data model.The three predictor datasets (bird distribution, bird functional traits and environmental data) clustered in a raster file were converted into a site-pair data format using the function formatsitepair from the 'gdm' package.This function inputs geographic and other predictor data to create a site-pair table, which is the format required to run generalized dissimilarity models.This new set of site-pair data contains distance (i.e.dissimilarity) and weight columns (i.e.importance of predictors) for each variable and site.These data were then included in a single generalized dissimilarity model that evaluated 1) bird functional traits, 2) bird distribution, 3) environmental features and 4) pairwise geographic distances among sites, as predictors of haemosporidian community composition variation.

Generalized dissimilarity modeling
In order to determine the drivers of variation in haemosporidian species community composition worldwide, we employed generalized dissimilarity models (GDM) (Ferrier et al. 2007, Mokany et al. 2022) using the function gdm from the 'gdm' package in R. In this model, parasite community composition dissimilarity among regions was the response variable; and geographic distance, bird distribution and the first axes of the PCAs comprising both environmental conditions and bird functional traits were used as explanatory variables.After fitting the GDM, we estimated the contribution of each of the four predictor variables toward explaining the deviation of parasite dissimilarity using variation partitioning analyses.For this, we used the function gdm.partition.deviance.Since there is a limitation of three variables for this function we ran the analyses twice.First, we evaluated the variation explained by host functional traits, host distribution and environmental features while excluding geographic distance.Second, with the same procedure, we explored the variation explained by bird distribution and functional traits clustered together, environmental features and geographic distance.Afterwards, we assessed the significance of each predictor for each model using the function gdm.varImp.This function estimates the deviance explained from the fitted model and compares it to the one derived from randomized models (n = 50).The given p-value represents the proportion of randomized models that explained better the deviance than the actual predictor data (Mokany et al. 2022).Finally, we evaluated the fit of each GDM using the function gdm.crossvaliation.Model fitting and individual predictors were evaluated visually and by summarizing the GDM.The height of the splines (i.e.function composed of pieces of simple functions joined to provide a suitable degree of smoothness to the GDM) illustrated the degree of parasite composition turnover among regions.All analyses were repeated for each parasite genus (i.e.Plasmodium, Haemoproteus and Leucocytozoon) separately.
One additional model incorporating beta phylogenetic diversity of hosts as a covariable was run.To estimate host beta phylogenetic diversity among regions we created a host occurrence matrix, downloaded a file containing 10 000 global bird phylogenetic trees (AllBirdsHackett1.tre)from https://birdtree.org/ and randomly selected one phylogenetic tree.Then, we calculated an index of beta-diversity based on the phylogenetic distances among species present in each region and used a PCA to reduce the dimensions of our data.However, since the addition of host phylogenetic relationships did not contribute to further explain the variance in our data, they were not incorporated in most of our analyses (Results and Supporting information).

Results
Haemosporidian richness varied between five and 276 distinct lineages per regions (i.e.defined as a geographic cell grid of 5 × 5 degrees, n = 100) and the maximum number of parasite records per region was 934.In general, we observed a mean of 36 distinct parasite lineages and 71 records per region.Regional bird richness (i.e.number of bird species per region) varied from four to 349, with a mean of 97 bird species.We found that geographic distance is the main predictor of spatial shifts in parasite community composition (p < 0.01), followed by host functional traits (p < 0.01), environmental features (p < 0.01) and, lastly, host geographic distribution (p = 0.02) (Fig. 2, Table 1A-C, model significance p < 0.01).To assess all environmental features and host functional traits in a unique model, dimension reduction was performed for both variables using PCAs.Interestingly, parasite species turnover responds differently to these predictors.For example, parasite turnover is mostly associated with temperature features (Fig. 3A), as mean temperature and temperature seasonality values (i.e. the variables that make the greatest contribution to negative and positive values in the PC1, respectively).Our results demonstrate that regions with warmer, rainier and less seasonal climate (negative values in Fig. 3A) show greater turnover of parasite lineages among localities than regions with lower temperature and precipitation, and high temperature seasonality (positive values in Fig. 3A).Dissimilarity in parasite community composition among warmer and wetter regions accounts for most of the turnover associated with environmental features (Fig. 2).
On the other hand, moderate to high host community composition dissimilarity is increasingly associated with turnover in parasite composition (Fig. 2F).Hence, regions with greater turnover of birds among localities must also present higher turnover of haemosporidians (Fig. 2F), which suggests parasite turnover increases in response to turnover of their avian hosts.Further, the main functional traits associated with haemosporidian turnover among regions were the presence of resident birds and strong territoriality; these were the traits that contributed the most to positive values in the PC1 (Fig. 3B).For those regions, shifts in the proportion of territorial and migratory birds between regions are associated with the greatest values of haemosporidian turnover.Moreover, regions harboring mostly large-bodied and/or migratory birds (i.e. the traits that present the greatest contribution to negative values in the PC1) showed the lowest rates of parasite turnover when compared to regions where birds present similar functional traits values in the PC1 (Fig. 3B).This suggests regions where migratory birds are abundant share more parasites among themselves (lower turnover) than regions harboring mostly resident birds.Naturally, regions that are the most distinct for any of the features considered also harbor the most distinct parasite fauna.In addition, we also observed a large overlap in deviance explanation among the variables analyzed, which is evidenced by the small increase in the variation explained when combining variables (Fig. 4, Table 1A, B and Supporting information).Overall, similar patterns were observed when analysing each parasite genus separately (Supporting information, models significance p < 0.01); however, when evaluating Haemoproteus and Leucocytozoon separately, bird distribution was not a significant predictor of haemosporidian turnover.
It is important to note that the addition of host beta phylogenetic diversity did not increase the amount of variation explained by our data (Supporting information).For this reason, our main models were run without taking into account the effects of host phylogenetic relationships on haemosporidian turnover.Again, the two major variables associated with haemosporidian turnover were geographic distance followed by host functional traits (p < 0.01).

Discussion
Species turnover and community assembly are ruled by three main filters, which include species' dispersal ability, abiotic tolerance and the pressures of interactions they face with other organisms (e.g.parasite-host interactions) (Vellend 2010).Here, we quantified the relative influence of these filters in driving parasite turnover worldwide using avian haemosporidians as a model.We found that geographic distance was the main driver of global haemosporidian turnover, indicating that dispersal ability may be the main constraint on parasite range.We also demonstrated that avian host functional diversity is a more important determinant of parasite turnover than host distribution (particularly in the case of Leucocytozoon parasites, Supporting information), since bird functional traits explained parasite turnover better than bird distribution in all models.Finally, we also revealed that climatic tolerance plays a minor role in shaping parasite species community composition worldwide, as environmental features (mainly temperature) explained only 19% of the variation in our main analyses.Nonetheless, our models explained only around 33-52% of total deviance in parasite community composition, suggesting there are many other filters associated with haemosporidian turnover, possibly related with vector distribution.Overall, we show that, despite the fact that abiotic and biotic filters seem highly correlated, dispersal ability and host resources are the principal forces shaping the spatial turnover of haemosporidian parasite species.Observed haemosporidian turnover as a function of GDM-predicted dissimilarity (i.e.expected turnover of haemosporidians).(C-F) I-splines illustrating the fitted relationship between GDM predictor variables (i.e.x-axis) and I-splines transformed variables (i.e.y-axis).Maximum height of each variable indicates the total degree of haemosporidian turnover predicted by that specific predictor.The shape of the splines illustrates the rate of turnover over each predictor gradient.Parasites rely on their hosts to complete their life cycle, with the presence of the latter in an area also associated with environmental features.Specifically, the range and turnover of parasites are also constrained by their hosts' tolerance to abiotic conditions (Mestre et al. 2020).Here, we demonstrated a great overlap in the deviance explained by host and environmental features.Since haemosporidians coevolve with their hosts (Pacheco et al. 2018, de Angeli Dutra et al. 2022), they are in turn constantly under selection to adapt to the environmental conditions those hosts inhabit or migrate into.Consequently, the linked effects of environmental and host features on parasite turnover might be the outcome of similar selective pressures for abiotic tolerance faced mutually by both haemosporidian parasites and their avian hosts through their evolutionary history.
Moreover, despite the large overlap in deviance explained by host functional traits and geographic distribution, the first emerges as a better predictor of parasite turnover.The two main functional traits associated with parasite turnover are the proportion of resident species in the local avian fauna and the degree of territoriality among those species.Regions inhabited by generally larger-bodied hosts present lower haemosporidian turnover.Host body mass is positively associated with both haemosporidian prevalence (Scheuerlein and Ricklefs 2004, Filion et al. 2020, Fecchio et al. 2021) and richness (Arriero and Møller 2008), therefore the presence of large-bodied hosts can shape local parasite species composition by providing a suitable resource for multiple haemosporidian species, enabling similar parasite lineages to colonize multiple communities.A high proportion of largebodied hosts in a broad region facilitates higher similarity of parasite lineages among localities within the region.On the other hand, migration can homogenize parasite composition through parasite dispersal (Møller and Szép 2011, Ellis et al. 2015, de Angeli Dutra et al. 2021b) and parasite-host network connectance (de Angeli Dutra et al. 2021c).A greater proportion of migrants -and, consequently, fewer resident birds -could also reduce parasite turnover because of the greater likelihood parasites will be transported among communities harboring multiple migratory host species.Thus, the mechanisms through which host functional traits structure parasite turnover vary, and distinct traits can have contrasting outcomes on parasite turnover.
Climate can also drive parasite turnover by shaping vector abundance and diversity (Benning et al. 2002, Atkinson et al. 2014, Clark et al. 2018, de La Torre et al. 2022).Here, we found that parasite turnover among regions with higher temperatures and precipitation rates, and low temperature seasonality, are the highest worldwide.This could be because haemosporidians use ectothermic hematophagous insects as vectors; as a result, these parasites must tolerate local temperatures to survive and reproduce (Valkiūnas 2005).Temperature is, therefore, the main driver of mosquito development due to its influence on metabolic rate (Chandrasegaran et al. 2020).
Local temperature drives the developmental rate of vectors of haemosporidian parasites (Mordecai et al. 2013) and unsuitable temperature conditions constrain parasite sporogonic development within their vectors (Lapointe et al. 2010).Thus, parasite turnover across localities in colder and drier regions with high temperature seasonality might be lower due to temperature placing constraints on parasite and vector development, limiting the diversity of parasites capable of persisting in those environments.Overall, environmental features have the potential to shape parasite turnover by constraining direct and indirect effects on parasite, host and vector development.
Previous research has evaluated factors driving haemosporidian turnover at local and regional scales (Williamson et al. 2019, McNew et al. 2021, de La Torre et al. 2022, Fecchio et al. 2022).Some found that geographic distance was not the principal driver of parasite turnover.Instead, researchers observed that host distribution (de La Torre et al. 2022) and environmental features (Williamson et al. 2019, McNew et al. 2021) were the main determinants of dissimilarity among parasite communities.Precipitation and elevation seem the most important drivers of haemosporidian parasite turnover in the Peruvian Andes (McNew et al.

2021
) and in North American sky islands mountain ranges (Williamson et al. 2019), respectively.On the other hand, de La Torre et al. ( 2022) also observed that host resources were the most important drivers of haemosporidian turnover in the Brazilian Atlantic Forest.However, contrary to the findings of this study, de La Torre et al. (2022) demonstrated that host community composition was a more relevant driver of parasite turnover than host functional traits.All these studies have explored haemosporidian turnover dynamics in a specific type of environment (e.g.Tropical Andes in Peru, sky island mountain ranges in the USA and Atlantic Rain Forest in Brazil); consequently, the variation in biotic and abiotic diversity expected among the sites evaluated is likely much lower than that seen at a global scale.This may explain the discrepancy between past research findings and the conclusions of the present study.Indeed, past research on large-scale parasite turnover has confirmed this pattern among haemosporidians and other parasite taxa (Dallas and Poisot 2018, Fecchio et al. 2019b, 2022, Martins et al. 2021).Therefore, geographic distance might only play an important role in shaping parasite turnover at larger scales or when there is a higher degree of heterogeneity among sites, whereas landscape and host features are the most important factors determining parasite turnover at smaller scales or among more homogeneous environmental conditions.
It is important to note, however, that the use and combination of multiple datasets comes with an increased risk of intrinsic errors, such as typos and divergence in species names (Nekola and Horsák 2022).Furthermore, data comprising bird distribution do not take into consideration the abundance of host species, which could be an important bias in our analyses.We also did not check for richness gradients directly; however, the addition of bird species distribution to the models indirectly takes this gradient into account, as bird diversity generally follows richness gradients (Jetz et al. 2012).Another relevant bias in our analyses is the concentration of haemosporidian samples in certain regions of the globe (e.g.Europe and Americas) whilst some other regions (e.g.Asia and Africa) lack data.Finally, due to absence of global data on the competent vectors of haemosporidians, our analyses could not consider vector distribution or their functional traits as well as parasite functional traits as potential determinants of haemosporidian turnover.
We demonstrated that, on a global spatial scale, avian hosts are the main predictor of the turnover of haemosporidian species across localities, after geographic distance (Table 1B).More specifically, we observed that haemosporidian species turnover among localities in regions harboring multiple migratory species is reduced, whereas this turnover is more pronounced in regions inhabited mostly by territorial and resident birds.We also demonstrated that temperature is the major climatic condition structuring haemosporidian turnover, with mean temperature values being the main factor associated with spatial shifts in haemosporidian composition.At the same time, since lower temperature values are associated with less pronounced haemosporidian turnover among localities within environmentally similar regions, temperate/high latitude regions should present more similar parasite community composition across localities than tropical regions.We also reveal that host functional traits, distribution and climate effects on haemosporidian parasite turnover are intrinsically inter-dependent.Finally, as geographic distance is the main predictor of haemosporidian species turnover, we conclude that dispersal limitations and ecological gradients may be major forces driving variation in parasite species composition globally, followed by resource (i.e.hosts) availability.

Figure 1 .
Figure 1.Sampled localities included in the analysis.Color scale illustrates spatial variation in bird species richness worldwide.Crosses depict parasite distribution from a total of 100 regions and 207 localities (including offshore islands) extracted from the MalAvi database (http://130.235.244.92/Malavi/).

Figure 2
Figure 2. (A)Observed haemosporidian turnover as a function of GDM predicted ecological distance (i.e.calculated based on the dissimilarity of all variables included in the GDM -geographic distance, environmental features, bird functional traits and distribution).(B) Observed haemosporidian turnover as a function of GDM-predicted dissimilarity (i.e.expected turnover of haemosporidians).(C-F) I-splines illustrating the fitted relationship between GDM predictor variables (i.e.x-axis) and I-splines transformed variables (i.e.y-axis).Maximum height of each variable indicates the total degree of haemosporidian turnover predicted by that specific predictor.The shape of the splines illustrates the rate of turnover over each predictor gradient.

Figure 4
Figure 4. (A) Diagram illustrating the proportion of deviance in haemosporidian parasites turnover explained by host functional traits, environmental features and host distribution.(B) Diagram illustrating the proportion of deviance in haemosporidian parasites turnover explained by host functional traits and distribution clustered, environmental features and geographic distance.

Table 1 .
(A) Deviance in haemosporidian parasite composition variation explained by bird distribution (p = 0.02), bird functional traits (p < 0.01) and environmental features (p < 0.01).(B) Deviance in haemosporidian parasite composition variation explained by bird distribution and functional traits (clustered as 'Bird data'), environmental features and geographic distance (p < 0.01).(C) Variable importance (i.e.percent change in deviance explained in the full model) and p-value., 2023, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ecog.06634 by Ministry Of Health, Wiley Online Library on [02/05/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Figure 3. Contribution of each variable of principal component analyses (PCA) for (A) environmental features and (B) host functional traits variability worldwide.Font size represents the variable's contribution to the first axis of the PCA.Variables with negative values are illustrated in red whereas blue is used for variables that had positive values in the PCAs.Territ = territoriality strength, dist = migration distance.See Supporting information for contribution values for each variable.Bio2, Bio15, distance of migration short and territoriality none are not shown because their contribution was less than 1%.Vertical position of variables is arbitrary.Information regarding environmental features and host functional traits was extracted from WorldClim and OpenTraits and Dufour et al. (2020), respectively.16000587 16000587, 2023, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ecog.06634 by Ministry Of Health, Wiley Online Library on [02/05/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License