Mutualism or Parasitism? Using a Phylogenetic Approach to Characterize the Oxpecker-Ungulate Relationship

an apparently durable mutualism involving two lineages of vertebrates

The evolution of mutualistic interactions represents a major question in evolutionary biology, with the long-term stability of such interactions potentially eroded when one or both of the partners "cheats" by failing to provide benefits to the other (Axelrod and Hamilton 1981;Bronstein 1994;Noë et al. 2001). A major challenge in studying mutualism is that the benefits can be difficult to quantify, thus complicating efforts to assess whether mutualistic interactions are occurring. One solution to this problem is to use data on the preferences of one partner for the other partner in relation to factors expected under competing hypotheses. Oxpeckers represent a unique opportunity in this regard, as this group includes two species (red-billed oxpecker, Buphagus erythrorhynchus and yellow-billed oxpecker, B. africanus) that both exhibit strong preferences for specific ungulate species (Mooring and Mundy 1996a;Koenig 1997). In the field, some studies have found that oxpeckers play a role in removing ticks and other ectoparasites from ungulates (Bezuidenhout and Stutterheim 1980;Mooring and Mundy 1996a, b). Previous comparative research found that ungulate body mass can account for variation in oxpecker preferences (Mooring and Mundy 1996a;Koenig 1997). More recent work suggests, however, that oxpeckers may be acting as parasites, ignoring available ticks and instead consuming host tissues, such as dead skin, blood, and earwax (Weeks 1999(Weeks , 2000McElligott et al. 2004).
While a considerable amount of observational data confirm that oxpeckers feed on wounds, (i.e., Weeks 1999Weeks , 2000Plantan 2009), it is difficult to determine exactly what oxpeckers are eating strictly from behavioral observations. Stomach analyses by Moreau (1933) and van Someren (1958) demonstrated that oxpeckers feed extensively on ticks of all sizes and stages of development, along with blood-sucking flies. In Moreau's (1933) sample of 58 birds, for example, the average number of ticks per bird was over 41. Moreau (1933) further speculated that the blood engorged by the ticks, rather than the tick tissue per se, forms the main food of red-billed oxpeckers; hence he suggested that the birds might naturally feed on blood from the ungulates. Subsequent work also can be interpreted in multiple ways. Stutterheim et al. (1988) conducted an experiment showing that caged red-billed oxpeckers significantly reduced the tick loads on two oxen, while other studies can be interpreted as indicating that feeding on open wounds (as well as on bodily chaff such as dandruff and earwax) is the birds' preferred food (Bezuidenhout and Stutterheim 1980). Weeks' (1999Weeks' ( , 2000 more recent work supporting the importance of blood feeding and failure to find significant tick reduction on cattle prompted him to question the mutualism hypothesis. Thus, we know that oxpeckers eat ticks and feed on wounds opportunistically (Plantan 2009), but we still do not know whether they are benefiting the hosts (and thus mutualistic) or exploiting them.
To investigate the degree to which mutualism or parasitism characterizes the oxpeckerungulate association, we compiled data on host preferences for both oxpecker species and levels of tick parasitism on different ungulate host species. Under the mutualism hypothesis, oxpeckers should prefer ungulate hosts that harbor higher abundance of ticks, measured as the number of ticks per average individual ungulate or social group of ungulates. If oxpeckers function as parasites by feeding primarily on host tissues, we predicted that oxpecker preferences should covary negatively with ungulate hide thickness because thinner hides should be easier to pierce, resulting in more wounds and making it easier for the birds to access host flesh. Although longterm coevolutionary interactions between parasitic oxpeckers and their hosts might favor thicker hides in more preferred hosts, we expect to find that current preferences correlate with thinner host hides. Moreover, cross-species variation in hide thickness among mammalian herbivores may covary more strongly with body mass or levels of intraspecific competition (Jarman 1989), rather than with parasitic selective pressures from oxpeckers. We also investigated whether oxpeckers prefer larger-bodied ungulates, possibly because larger-bodied ungulates offer a larger or more stable environment for foraging or because larger ungulates have greater tick abundance (Ezenwa et al. 2006), and whether they prefer ungulates that live in more open environments, which would facilitate perching by the birds.

Data on tick infestations for ungulate hosts were obtained from the Global Mammal Parasite
Database (Nunn and Altizer 2005;Ezenwa et al. 2006). The data on tick abundance are measured as counts (rather than density) on individual hosts, and were compiled from publications detailing patterns of parasitism in wild ungulates from 1981 to 2002. Oxpecker preferences were measured as the mean across studies of the log-transformed host preference index calculated for each oxpecker-host combination, which was calculated as number of oxpeckers on a particular host species at a site divided by the number of hosts at that site and multiplied by 1000. Across the ungulate species in the dataset with measures from at least two localities, the association of host preference indices in pairs of localities were correlated (mean + SE Spearman rank correlation for red-billed oxpeckers: 0.57 + 0.07, n=27 comparisons; for yellow-billed oxpeckers: 0.71 + 0.06, n=21 comparisons; data from Koenig 1997 We aimed to investigate preferences in relation to tick abundance, which refers to the average number of ticks on individual hosts of each ungulate species, including those with zero ticks. The available literature, however, generally provides data on intensity (average number of ticks on infected individuals) and prevalence (percentage of animals with at least one tick, see Nunn and Altizer 2006). To obtain tick abundance per individual ungulate host species, we therefore multiplied intensity by prevalence. We identified 19 ungulate host species from Koenig (1997) with data on ticks, and of these, estimates of tick abundance were available for 14 species. To calculate group-level abundance we multiplied individual abundance by average group size for the ungulate species. The data are provided in the Supplementary Information.
Because the traits of different ungulate species are shared through common descent, the data on preferences and other traits are not necessarily statistically independent (Harvey and Pagel 1991;Nunn and Barton 2001). A phylogenetic supertree is available for the species in our dataset (Price et al. 2005;Bininda-Emonds et al. 2007). However, phylogenetic relationships and branch lengths are never known with certainty. We therefore implemented new methods to systematically control for phylogenetic uncertainty in comparative research (Pagel and Lutzoni 2002). We also used phylogenetic comparative methods that estimate phylogenetic signal and take the degree of phylogenetic signal into account when testing predictions (Pagel 1999).
Specifically, we used a Bayesian Markov chain Monte Carlo (MCMC) approach to infer multiple trees and then ran Bayesian MCMC comparative analyses using the resulting posterior probability distribution of trees (Huelsenbeck et al. 2000;Pagel and Lutzoni 2002;Ronquist 2004).
To generate the trees, we extracted data from GenBank for three mitochondrial genes (cytochrome B, 12S rRNA and 16S rRNA) for the 24 ungulate species with data on oxpecker preferences (see Supplementary Information). We used the program ModelTest (Posada and Crandall 1998) to identify the substitution model that best describes the data using the AIC criterion (with a maximum likelihood optimized tree as base tree for the likelihood calculations).
We found that a general time reversible model with gamma distributed rate variation among sites and a proportion of invariable sites (GTR+I+G) provided the best model for each of the three genes. To create the multiple sequence alignments (MSA), we used Muscle 3.7 (Edgar 2004).
Because alignment quality can have a substantial impact on the inferred tree (Ogden and Rosenberg 2006;Talavera and Castresana 2007), we used the program GBlocks (Castresana 2002) with the settings -b5=h, -t=d, and -b2=0.6 * number of sequences to exclude poorly aligned sites or sites with a high percentage of missing data (especially at the beginning and end of the MSA), which significantly improved alignment quality. We constrained two perissodactyls (Ceratotherium simum and Equus burchellii) to be sister species because, in analyses that did not constrain this node, the limited sequence data available for Equus burchellii produced higher levels of uncertainty, but this uncertainty did not reflect current understanding of ungulate phylogeny.
We ran three independent runs of 10,000,000 generations each in MrBayes version 3.1.2 (Ronquist and Huelsenbeck 2003), with the domestic cat (Felis catus) identified as the outgroup and sampling the chain every 2,000 generations. We used Metropolis coupled MCMC with six chains (one cold chain and five heated chains) to improve sampling of tree space. In all three independent runs, we excluded the first 500 sampled trees as burn-in, giving a total of 13,500 phylogenies (4,500 in each run) as an estimate of the posterior distribution of ungulate phylogeny for the species with oxpecker preference scores. We found high convergence among the three runs, with the potential scale reduction factor (PSRF) equaling (or very close to) 1 for all model parameters and an average deviation of split frequencies among the three runs equal to 0.0012 at the end of the MC analysis. For illustrative purposes, we summarized these topologies by constructing a 50% majority rule consensus tree (Fig. 1), but we ran analyses on individually sampled trees after pruning the outgroup from each tree.
We analyzed the log-transformed data across these trees using the MCMC regression model in the program BayesTraits (Pagel and Meade 2007). We ran three independent chains for 3,050,000 iterations, sampling every 100 generations and excluding the first 50,001 iterations as burn-in, with uniform priors on regression coefficients ranging from -100 to 100. One of the 13,500 trees was selected per iteration of the MCMC regression analysis. To ensure adequate sampling of parameter space, we adjusted the "ratedev" parameter so that the average acceptance rate fell between 0.2 and 0.4 (Pagel and Meade 2007). By examining the parameter values (including regression coefficients) across chains, we ensured that the MCMC sampling had converged on the posterior probability distribution, which it clearly had in all cases. Credible intervals were set as the distribution of the middle 95% of the distribution of sampled regression coefficients from the first MCMC run, and statistical support was inferred when the credible interval excluded zero. For some analyses, we provided the proportion of positive or negative coefficients to give further insight to support for the predictions.
We also used BayesTraits to estimate λ as a measure of phylogenetic signal (Freckleton et al. 2002). This parameter scales the off-diagonal elements of the variance-covariance matrix representing the phylogeny, where the off-diagonals are internal branch lengths. A value of 0 indicates no phylogenetic signal (i.e., a star phylogeny because all internal branches are collapsed to zero length). A value of 1 indicates that the pattern is consistent with a Brownian motion model of evolution using the available branch lengths. We estimated λ within the regression model, rather than for each variable separately; thus, estimates of λ are for residuals from the regression model (Lavin et al. 2008).

Results
Our phylogenetic analyses revealed a strong association between individual tick abundance and oxpecker preferences, and this was true for both species of oxpeckers ( Figure 2) and across phylogenies (Table 1). Our analyses further revealed that λ varied widely but is on average greater than zero and less than one. This finding argues strongly for using phylogeny-based methods when testing this prediction. We also found that both species of oxpeckers tend to prefer the same ungulate host species (regression of red-billed preferences on yellow-billed preferences: credible interval for slope estimate is 0.63 to 1.14; R 2 =0.74; credible interval for λ is 0.18 to 0.99).
Because many ungulates are gregarious, the relevant measure of tick abundance might be at the level of ungulate social groups rather than individuals, with oxpeckers selecting among "patches" of ungulates based on the number of ticks typically found in a social group. Grouplevel tick abundance yielded significant results for both species of oxpeckers, with group-level abundance strongly predicting the degree to which oxpeckers preferred different ungulate species (Table 1). We also investigated whether ungulates living in larger groups have higher tick abundance. The analyses provided limited support for this possibility, with regression coefficients positive in 94% of the samples from the MCMC chain, although the credible interval encompassed 0 (ranging from -0.19 to 1.42; R 2 =0.20; mean λ=0.66). (Table 1, see also Koenig 1997), and body mass and individual tick abundance associated positively among the species in our dataset (R 2 =0.34; credible interval on the regression coefficient: 0.04 to 0.85; credible interval on λ: 0.10 to 0.98). Using a multiple regression model in BayesTraits (Pagel and Meade 2007), we investigated the relative importance of body mass and tick abundance on oxpecker preferences. These analyses revealed that both body mass and tick abundance generally accounted for variation in oxpecker preferences, especially for red-billed oxpeckers (Table 2).

Both species of oxpeckers preferred larger-bodied hosts
Thus, the effect of tick abundance on oxpecker preferences was statistically independent of ungulate body mass.
Another potential confounding factor is habitat. It could be that ungulates living in more wooded or dense habitats would be less preferred hosts for oxpeckers, with overhanging vegetation tending to make it difficult for the birds to remain perched on ungulates. We found a slight tendency for both red-billed and yellow-billed preferences to covary with increasing habitat density, but the tendency was positive rather than negative as predicted (Table 1, only 7.1% and 15.9% of regression coefficients were negative, for red-billed and yellow-billed oxpeckers, respectively). In addition, hosts living in denser vegetation may be exposed to more ticks (e.g., see Mooring et al. 2004). We therefore tested for an association between habitat and individual tick abundance, but found that credible intervals included zero (credible interval: -0.18 to 0.52; R 2 =0.09; mean λ=0.50). In a multiple regression model, individual tick abundance was a strong predictor of oxpecker preferences for both species (credible interval for red-billed oxpeckers: 1.25 to 3.72; credible interval for yellow-billed oxpeckers: 0.82 to 3.54), while credible intervals for the effects of habitat included zero for both oxpecker species (credible interval for red-billed oxpeckers: -0.90 to 0.51; credible interval for yellow-billed oxpeckers: -0.86 to 0.78).
In contrast to results involving tick abundance, we found no support for a negative effect of hide thickness on oxpecker preferences, with credible intervals clearly encompassing zero for both species of oxpeckers (Table 1). It could be that tick abundance covaries with hide thickness, confounding this result. However, hide thickness did not predict individual abundance of ticks (credible interval: -1.23 to 1.65; R 2 =0.03; mean λ=0.41). We also investigated whether larger bodied ungulates have thicker hides, and found support for this possibility with a 95% credible interval that just barely excluded zero (credible interval: 0.001 to 0.86; R 2 =0.34; mean λ=0.51). Thus, the preference for larger-bodied hosts appears to be counter to expectations based on the parasitism hypothesis, which predicts oxpecker preferences for smaller-bodied ungulates with thinner hides.

Discussion
We found that oxpeckers exhibited stronger preferences for ungulates with higher tick burdens.
This association held when using preference data from both species of oxpeckers, and the association remained statistically significant after controlling for ungulate body mass. Moreover, by using recent advances in phylogenetic comparative methods, we were able to control for phylogenetic uncertainty in the evolutionary relationships and divergence times of the ungulate species in our sample. These analyses also revealed a moderate degree of phylogenetic signal in the evolutionary patterns that we examined. Collectively, these results strongly support the hypothesis that the oxpecker-ungulate relationship is mutualistic and provide a compelling argument for incorporating phylogeny when investigating the correlates of oxpecker preferences.
The mutualism and parasitism hypotheses are not mutually exclusive, however, and interpreted in light of recent field studies, we suggest that the mutualistic interaction between oxpeckers and their ungulate hosts may be sensitive to the underlying environmental context and can shift to parasitism under certain conditions. For example, although oxpeckers in one study frequently fed on cattle wounds, they tended to exploit existing lesions rather than create new ones (Weeks 2000). Thus, when ungulate hosts have readily available and abundant scratches, the benefits of feeding on wounds may exceed the benefits of foraging for ticks. Another study documented oxpecker preferences for feeding on wounds in captive black rhino (McElligott et al. 2004). In this case, the birds opened up new wounds, but the absence of ticks in the captive environment of this study may explain the switch to parasitic behavior (e.g., Plantan 2009). In addition, one might expect that greater risks of injury during intraspecific competition could favor thicker hides, but that oxpeckers may even prefer these hosts for tissue feeding on the occasions when wounds are present. This could weaken our expected negative prediction between hide thickness and oxpecker preferences.
Thus, oxpeckers may be more likely to ingest host tissue when they feed on non-preferred hosts, on captive hosts with lower tick burdens, or in species where wounding is common through intraspecific competition. Similarly, within a species, we might expect to see more parasitism on members of the sex that are most commonly involved in agonistic interactions. In addition, for the reasons just given, we expect that tests of the prediction involving oxpecker parasitism will provide more ambiguous comparative results, as compared to tests of the mutualism hypothesis. Despite these caveats, our results raise the intriguing possibility that ungulate hosts with lighter tick burdens, such as hartebeest (Alcelaphus buselaphus) and warthogs (Phacochoerus aethiopicus), may fall victim to "cheating" by oxpeckers more frequently. Our results therefore provide testable predictions for future field research.
Our study also provides further insights to previous research that links oxpecker preferences to host body mass (Koenig 1997). For red-billed oxpeckers, body mass-host preference associations appear to be due to covariation between body mass and individual tick abundance, with individual tick abundance being the key factor driving host preferences. By comparison, both individual tick abundance and body mass accounted for significant variation in yellow-billed oxpecker host preferences (see Table 1). For group-level patterns, group-level tick abundance was the primary predictor of red-billed oxpecker preferences, while body mass was the primary predictor of yellow-billed oxpecker preferences. This difference between the two oxpecker species may be driven by differences in their body size. Yellow-billed oxpeckers are larger, and previous work suggested that this size differential might account for slight differences in feeding behavior (Stutterheim et al. 1988). As the larger of the two oxpeckers, yellow-billed oxpeckers may require larger hosts for efficient foraging. Moreover, larger bodied ungulate hosts may support more individual birds, thus tending to increase the preference index that we used.
In summary, our analyses support the hypothesis that oxpeckers and their hosts exhibit a mutualistic relationship. Further strengthening this conclusion, our results are robust to uncertainty in the phylogenetic relationships of ungulate hosts and the underlying model of trait evolution, and the results are consistent in both oxpecker species. In contrast, we found no evidence in support of the parasitic hypothesis. Based on our findings and field evidence, we suggest that parasitic behavior by oxpeckers most likely occurs opportunistically within a generally mutualistic relationship. More generally, we provide new insights to the evolution of an apparently durable mutualism involving two lineages of vertebrates, and we show how new phylogenetic comparative methods can be used address fundamental questions involving mutualism.   The dependent variable was oxpecker preferences, with analyses conducted separately for each species. Ranges represent 95% credible intervals on the slope, and λ is a measure of phylogenetic signal (Freckleton et al. 2002). Credible intervals that exclude 0 are considered statistically meaningful.