Aim: Parasites play a fundamental role in shaping ecological communities and influencing trophic interactions. Understanding the factors that drive parasite impacts on community structure and stability (i.e., resilience to disturbances) is crucial for predicting disease dynamics and implementing effective conservation strategies. In this study, using avian malaria and malaria-like parasites as a model system, we investigated the relationship between specificity, community stability, and parasite vulnerability and their association with host diversity and climate.

Location: Global

Time period: 2009-2023

Major taxa studied: Avian malaria and malaria-like parasites.

Methods: By compiling occurrence data from a global avian haemosporidian parasite database (MalAvi), we constructed a comprehensive dataset encompassing 60 communities. We utilized a phylogenetic model approach to predict missing host-parasite interactions, enhancing the accuracy of our analyses. Network analyses based on bipartite interactions provided measures of network specificity, stability, modularity, parasite competition, and vulnerability to extinction.

Results: We found that the high network specificity reduced community stability and decreased competition among parasites. Furthermore, we found that parasite vulnerability decreased with increasing community stability, highlighting the importance of community stability in host-parasite interactions for long-term parasite persistence. When exploring the influence of local host diversity and climate conditions on host-parasite community stability, we demonstrated that increasing host biodiversity and precipitation reduces parasite competition. Conversely, higher temperature raises competition among parasites.

Conclusion: These findings provide valuable insights into the mechanisms underlying parasite impacts on communities and the interplay between specificity, community stability, and environmental factors. Further, we reveal the role of climate in shaping host-parasite interactions. By unravelling the complexities of parasite-mediated interactions, our research substantially improves the current knowledge of the importance of specificity as a modulator of interactions in bipartite networks.

Data compilation       

In June 2023 we extracted haemosporidian occurrence data from MalAvi. MalAvi (http://130.235.244.92/Malavi/) is an online, freely available global database on avian haemosporidian parasites (Plasmodium, Haemoproteus, and Leucocytozoon) (Bensch et al., 2009). The database comprises >15,000 entries, each representing a recorded malaria infection in an individual bird, with information on parasite lineage identity, bird species identity, and geographical location. Each entry in our dataset consisted of an individual link of an identified malaria lineage in an identified host species from a given geographical locality. However, the absence of a parasite lineage from a particular host species may reflect inadequate sampling, rather than a true host-parasite incompatibility. 

Network analyses and specificity

All analyses were performed in the R environment version 4.3 (R Core Team, 2023). Each community was composed of all parasite lineages and host species occurring within a given geographical locality including the potential missing links. Communities were defined as cell grids of equal (5x5 degrees, total = 2592) dimensions harboring 20 or more distinct parasite lineages and at least 10 distinct host species. Grids harboring fewer than 20 distinct parasite lineages, or less than 10 host species were excluded from our analyses. It is important to note that the majority of the grids in our study represent regions with scarce sampling (e.g., Africa and Asia) and/or regions with no haemosporidian parasites (e.g., oceans and Antarctica) as they cover the entire globe. We chose the cell grid size that provides the optimal compromise between the number of cells that met the inclusion criteria mentioned above, and the spatial extent that could reasonably be expected to allow for species interactions. We acknowledge that the area covered by a 5X5 degrees cell encompasses diverse habitats and that not all bird species will occupy the entire cell; thus, it is at the upper end of what may be considered a community. The use of geographical grids to delimit communities is consistent with previous studies on avian haemosporidian parasites (de Angeli Dutra et al., 2023) and well established in ecology (Bousquin, 2021; Rangel et al., 2018).

As detailed below, network specificity (H₂’) and stability (measured as nestedness and modularity), parasite inter-specific competition, and vulnerability to extinction (measured as niche overlap and extinction slope, respectively) were calculated using the package “bipartite” in R (Dormann et al., 2008). Network modularity was calculated by comparing observed modularity with the partition that maximizes modularity (Beckett, 2016).

Network specificity (H₂’) was estimated as the degree of specialization or partitioning among both host and parasite species within a bipartite network based on their interaction patterns (Blüthgen et al., 2006). This metric is useful to compare the degree of specialization across distinct interaction webs as it considers both parties of a bipartite network. Likewise,  nestedness is also a network-level index (i.e., considers both parties, hosts, and parasites, of a bipartite network). It reflects the degree to which specialist parasites interact with a subset of the hosts used by progressively more generalist parasites. Here, we calculated nestedness based on decreasing fill (or DF) and paired overlap (Almeida-Neto et al., 2008). Since nested networks are cohesive and robust and nestedness is a network-level measurement (Bascompte, 2007), we used this metric to estimate community stability. Additionally, in certain cases, modularity can be associated with higher stability (Grilli et al., 2016), therefore we also evaluated the role of network specificity on modularity. Here, we calculated modularity using label propagation and multi-step agglomeration to attempt to maximize modularity in binary bipartite networks applying the LPAwb+ algorithm (Beckett, 2016).  

On the other hand, niche overlap and extinction slope metrics only take into account one party of the bipartite networks: the parasites. Niche overlap measures the similarity in resource utilization between two or more species within a community (Dormann et al., 2008). In our analyses, it measures the degree to which parasites share host species. As a result, communities with low overlap in host utilization by parasites are subject to generally lower levels of interspecific competition. The extinction slope estimates the impact of removing a species from a bipartite network. It measures the coextinction of parasites as a consequence of repeated removals of host species from the network (Dormann et al., 2008). In the selected algorithm, the host species infected by fewer parasites are removed first. The present algorithm was selected since those species are expected to have the weakest impact on community stability and extinction slope is a metric of parasite vulnerability to their hosts’ extinction.  

All network metric values were standardized among distinct communities using null models and Z-scores (i.e., a numerical measurement that describes the position of a given value according to the mean) to account for the fact that the different community networks differ widely in the number of nodes (diversity of hosts and parasites) and the number of interactions (biases in sampling effort and/or prevalence). This method corrects for those biases as we used null models to create a thousand different and reshuffled networks. Our null models created shuffled binary matrices, assigning interspecific interactions according to species-specific probabilities based on abundance of host species while ensuring that all hosts and parasites still had at least one interaction. This procedure was repeated a thousand times for each host-parasite network. All those new matrices were used to calculate randomized network metrics (e.g., nestedness, niche overlap, extinction slope, and specificity). After we obtained a thousand random values of each network metric for each host-parasite network, we used them to calculate Z-scores. Therefore, the final network metric values obtained and used in the analyses represent how the original value compares to randomly generated values, i.e., if the parasites compete as much, are as likely to get extinct or whether our networks are more, less, or equally specific and stable, as they would be by chance alone based solely on their number of nodes and interactions. 

Local biodiversity and climate conditions

Local host taxonomic and functional diversity and climatic conditions were calculated for each community evaluated. Host taxonomic diversity was calculated based on avian occurrence data extracted from the BirdLife International database. Host functional traits for each bird species were extracted from the Open Traits Network database (https://opentraits.org/datasets.html) (Wilman et al., 2014), migratory status followed classification by Dufour et al., 2020. We used Hill numbers to estimate both diversity metrics using the package hillR in R (Chao et al., 2014). Hill numbers are a common metric of diversity used in community ecology. Taxonomic diversity is measured as Shannon entropy. In our study, it refers to the level of uncertainty in identifying the host species of a randomly selected parasite. On the other hand, functional diversity was calculated as the sum of the pairwise distances (i.e., variety in host functional attributes evaluated: migratory status, territoriality, body weight, diet and habitat) among the host species that parasites infect in each cell grid (Chao et al., 2014).

We extracted climate data obtained from Worldclim (https://worldclim.org/) at a resolution of 10km for the whole grid area. To do so, we used the "getData" function from the "raster" package in R climate data (Hijmans, 2023). Values obtained across each grid were averaged so that our climatic metrics took into consideration conditions for the entire cell. The Worldclim dataset consists of 19 distinct environmental attributes, encompassing measurements related to temperature and precipitation. To reduce the dimensions of our climatic data, we ran two Principal Component Analysis (PCAs) separating the variables for temperature variables (Bio1-11) and precipitation variables (Bio12-19). PCAs for temperature and precipitation were performed separately to provide a clearer interpretation of the principal components for each variable. This is important in understanding the dominant patterns of variability in temperature and precipitation independently. We added the first axis of both PCAs (which explained 78% and 60% of the variance, respectively) as explanatory variables in the Bayesian models.