Ecological and evolutionary drivers of haemoplasma infection and bacterial genotype sharing in a Neotropical bat community.

Most emerging pathogens can infect multiple species, underlining the importance of understanding the ecological and evolutionary factors that allow some hosts to harbour greater infection prevalence and share pathogens with other species. However, our understanding of pathogen jumps is based primarily around viruses, despite bacteria accounting for the greatest proportion of zoonoses. Because bacterial pathogens in bats (order Chiroptera) can have conservation and human health consequences, studies that examine the ecological and evolutionary drivers of bacterial prevalence and barriers to pathogen sharing are crucially needed. Here were studied haemotropic Mycoplasma spp. (i.e., haemoplasmas) across a species-rich bat community in Belize over two years. Across 469 bats spanning 33 species, half of individuals and two-thirds of species were haemoplasma positive. Infection prevalence was higher for males and for species with larger body mass and colony sizes. Haemoplasmas displayed high genetic diversity (21 novel genotypes) and strong host specificity. Evolutionary patterns supported codivergence of bats and bacterial genotypes alongside phylogenetically constrained host shifts. Bat species centrality to the network of shared haemoplasma genotypes was phylogenetically clustered and unrelated to prevalence, further suggesting rare-but detectable-bacterial sharing between species. Our study highlights the importance of using fine phylogenetic scales when assessing host specificity and suggests phylogenetic similarity may play a key role in host shifts not only for viruses but also for bacteria. Such work more broadly contributes to increasing efforts to understand cross-species transmission and the epidemiological consequences of bacterial pathogens.

INTRODUCTION

Most pathogens that cause disease in humans, domestic animals and wildlife are capable of infecting multiple host species (Woolhouse, Taylor, & Haydon, 2001). However, predicting which hosts maintain pathogens and identifying their role in cross‐species transmission can be challenging, as many hosts can be infected but not play key roles in the reservoir community (Fenton, Streicker, Petchey, & Pedersen, 2015). Pathogen jumps between species depend on infection prevalence in the donor host, transmission opportunities between donor and recipient species, and suitability of the recipient host for pathogen replication (Plowright et al., 2017). Each of these steps can be shaped by ecological and evolutionary factors (VanderWaal & Ezenwa, 2016). For example, small‐bodied species can have greater competence, the ability to transmit new infections, than larger species (Downs, Schoenle, Han, Harrison, & Martin, 2019), and host switching is often constrained by phylogeny, owing to similarity in immunological barriers to pathogen replication between closely related species (Streicker et al., 2010). Identifying the ecological and evolutionary factors that allow some species to harbour greater prevalence and have facilitated pathogen sharing can improve our general understanding of disease emergence (Fountain‐Jones et al., 2018). Examining evolutionary associations between hosts and pathogens can further uncover factors favouring host shifts versus codivergence and assess the frequency of cross‐species transmission (Geoghegan, Duchêne, & Holmes, 2017).

Given the public health and agricultural burdens of many zoonotic pathogens such as avian influenza viruses, henipaviruses and lyssaviruses, many investigations of pathogen prevalence and emergence focus on viruses (Geoghegan et al., 2017; Olival et al., 2017). However, more zoonoses are caused by bacteria than other pathogen taxa (Han, Kramer, & Drake, 2016), and bacterial pathogens can negatively impact newly infected host species (e.g., Mycoplasma galliscepticum, a poultry pathogen, caused rapid population declines in wild house finches; Hochachka & Dhondt, 2000). More attention to bacteria and their propensity for host specificity versus generalism is accordingly important for understanding whether factors that govern cross‐species transmission of viruses can be extended to other pathogens (Bonneaud, Weinert, & Kuijper, 2019). Bacterial pathogens have been especially understudied for bats (Mühldorfer, 2013), in contrast to intensive studies of zoonotic viruses across the Chiroptera (Luis et al., 2015). However, many bacterial pathogens are probably important to both bat conservation and human health due to pathogenic effects on bats themselves as well as their zoonotic potential (Becker et al., 2018; Evans, Bown, Timofte, Simpson, & Birtles, 2009).

To determine the ecological and evolutionary drivers of bacterial prevalence and barriers to pathogen sharing, we focused on haemotropic Mycoplasma spp. (i.e., hemoplasmas) in a species‐rich bat community in Belize (Fenton et al., 2001; Herrera, Duncan, Clare, Fenton, & Simmons, 2018). The Neotropics have remarkable bat diversity owing to adaptive radiation in the Phyllostomidae (Gunnell & Simmons, 2012), producing a range of feeding strategies (e.g., frugivory, carnivory, sanguivory), body sizes and roost preferences (Monteiro & Nogueira, 2011). Haemoplasmas are facultative intracellular erythrocytic bacteria transmitted through direct contact (Cohen et al., 2018; Museux et al., 2009) and also possibly via arthropod vectors (Willi, Boretti, Meli, et al., 2007). Haemoplasmas can cause acute and chronic anaemia, especially for immunocompromised hosts; however, many animals develop inapparent infections and are asymptomatic (Messick, 2004). As Mycoplasma spp. lack many of the metabolic pathways associated with energy production and synthesis of cell components found in other bacteria, they are fully dependent on host cells (Citti & Blanchard, 2013). Haemoplasmas have therefore been described as mostly host specialists (Pitcher & Nicholas, 2005), yet interspecies and potentially zoonotic transmission can occur (Maggi et al., 2013; Willi, Boretti, Tasker, et al., 2007). Haemoplasmas are common and genetically diverse in bats (Di Cataldo, Kamani, Cevidanes, Msheliza, & Millán, 2020; Ikeda et al., 2017; Mascarelli et al., 2014; Volokhov, Becker, et al., 2017), which offers an ideal model system for identifying the ecological and evolutionary factors structuring bacterial infection risks within and between host species.

Many cross‐species comparisons of pathogen infection risks and sharing have used less‐diverse host communities (Johnson et al., 2012; VanderWaal, Atwill, Isbell, & McCowan, 2014) or global data sets of host–pathogen associations that can be limited by heterogeneous sampling effort and variation in pathogen detection methods (Dallas et al., 2019; Huang, Bininda‐Emonds, Stephens, Gittleman, & Altizer, 2014). Our focus on a widespread pathogen group in a highly diverse host community allowed us to capitalize on strong host trait variation while controlling for sampling effort and diagnostic methods (Becker, Crowley, Washburne, & Plowright, 2019; Han, Kramer, et al., 2016). Past work has also used host–pathogen networks to characterize contemporary or historical transmission at often coarse taxonomic scales (e.g., pathogen species complexes or genera; Blyton, Banks, Peakall, Lindenmayer, & Gordon, 2014; VanderWaal et al., 2014). However, as bat species can be infected by multiple haemoplasma genotypes, and because genotypes with ≥99% sequence identity of their 16S rRNA genes can represent different bacterial species (Volokhov, Hwang, Chizhikov, Danaceau, & Gottdenker, 2017; Volokhov, Simonyan, Davidson, & Chizhikov, 2012), focusing on genotypes provides finer‐scale resolution to determine the ecological and evolutionary features of species that facilitate pathogen sharing and to identify likely maintenance hosts of bacterial infections (Fountain‐Jones et al., 2018).

We investigated three questions about the relative contribution of ecological traits and evolutionary history to structuring infection patterns and pathogen sharing. First, what are the individual and ecological predictors of haemoplasma infection in a Neotropical bat community? Second, does the distribution of haemoplasma genotypes across the bat community map onto the bat phylogeny, as predicted by host–pathogen codivergence? Third, if genotype sharing among host species occurs, which host clades or traits best predict species ability to share pathogens? We predicted that ecological covariates such as ectoparasitism and large host colonies could increase bacteria risk through vector‐borne or density‐dependent transmission (McCallum, Barlow, & Hone, 2001; Willi, Boretti, Meli, et al., 2007). We also expected haemoplasma genotypes would be specific to particular host species but that more closely related bats would share haemoplasma genotypes, indicating phylogenetically restricted host shifts (Pitcher & Nicholas, 2005). However, ecological traits that increase the risk of pathogen exposure between species, such as occupying a greater diversity of roosting habitats, could also facilitate pathogen genotype sharing among less closely related hosts (McKee et al., 2019).

METHODS

Bat capture and sampling

From April 24 to May 6, 2017 and from April 23 to May 5, 2018, we sampled 469 bats from 33 species captured in two adjacent areas in the Orange Walk District of Belize: Lamanai Archaeological Reserve (LAR) and Ka’Kabish (KK; Figure 1). The LAR is bordered by the New River Lagoon, forest and agriculture, while KK is a remnant patch of forest surrounded by agriculture located 10 km away. We consider these sites to be independent, as the small home ranges of many Neotropical bat species and lack of continuous forest between LAR and KK have probably restricted most individual bat intersite movement (Jones, Hämsch, Page, Kalko, & O’mara, 2017; Loayza & Loiselle, 2008); however, some species may have greater connectivity (e.g., inter‐roost movement of Desmodus rotundus occurs rarely). At least 44 of the 70 bat species in Belize have been recorded in this region (Herrera et al., 2018; Reid, 1997). Bats were captured with mist nets primarily along flight paths and occasionally at the exits of roosts from 7 p.m. until 10 p.m. Harp traps were also set from 6 p.m. to 5 a.m. In KK, D. rotundus, Trachops cirrhosus and Chrotopterus auritus were captured at their shared roost. In the LAR, D. rotundus, Saccopteryx bilineata and Glossophaga soricina were sampled from a shared roost, although individuals were also sampled along flight paths across the greater site. More broadly, patterns of roost sharing of bat species in northern Belize remain elusive.

Figure 1

Study sites in northern Belize. The shaded inset shows the location of Orange Walk District. Borders show the boundaries of the LAR (Lamanai Archaeological Reserve) and KK (Ka’Kabish). White and brown shading indicates agricultural and urban development, while dark green shading represents intact forest. Satellite imagery was derived from Google Maps. Stacked bar plots show the relative abundance of each sampled bat family per study site [Colour figure can be viewed at wileyonlinelibrary.com]

Bats were placed in cloth bags until processing and were identified to species (and sex) based on morphology (Reid, 1997). Reproductive activity was indicated by the presence of scrotal testes in males and by evidence of pregnancy or lactation in females; across bat species, 69% of males and 42% of females were in reproductive condition. We also visually screened bats for the presence of ectoparasites (i.e., bat flies, ticks, bat bugs, mites; Ter Hofstede, Fenton, & Whitaker, 2004). We collected 3–30 µl of blood (volumes were dependent on bat mass) by lancing the propatagial vein with a sterile needle. Blood was collected with heparinized capillary tubes and stored on Whatman FTA cards to preserve bacterial DNA. Field procedures followed guidelines for the safe and humane handling of bats published by of the American Society of Mammalogists (Sikes, The Animal Care and Use Committee of the American Society of Mammalogists, 2016) and were approved by the Institutional Animal Care and Use Committees of the University of Georgia (A2014 04‐016‐Y3‐A5) and American Museum of Natural History (AMNHIACUC‐20170403 and AMNHIACUC‐20180123). Fieldwork was authorized by the Belize Forest Department under permits WL/2/1/17(16), WL/2/1/17(19) and WL/2/1/18(16). Sample size was similar between years (2017 = 202, 2018 = 267) but varied by site (LAR = 365, KK = 101). More species were sampled for blood at LAR (n = 33) than at KK (n = 17; Figure 1), reflecting site differences in species richness (Herrera et al., 2018). We sampled 1–139 individuals per bat species (the maximum being D. rotundus), with a mean of 14 individuals per bat species (Table S1).

DNA extraction, PCR amplification and amplicon sequencing

Genomic DNA was extracted from blood on FTA cards using QIAamp DNA Investigator Kits (Qiagen). We tested DNA for haemoplasmas using PCR (polymerase chain reaction) with primers and procedures described in prior analyses (Volokhov, Becker, et al., 2017). We included blank FTA punches as an extraction control, ultrapure water as a negative control and Candidatus Mycoplasma haemozalophi DNA as a positive control (Volokhov et al., 2011). Amplicons from PCR‐positive samples were purified by electrophoresis and extracted using the QIAquick Gel Extraction Kit (Qiagen).

To determine haemoplasma infection status, all 16S rRNA amplicons were directly sequenced by Macrogen (https://www.macrogenusa.com/). Amplicons were sequenced with the same primers used for PCR amplification and then with internal (walking) primers when needed (Volokhov, Becker, et al., 2017). Negative DNA samples were tested for amplification quality using the universal PCR primers targeting the mammal mitochondrial 16S rRNA gene (Volokhov, Kong, George, Anderson, & Chizhikov, 2008) or the mitochondrial cytochrome c oxidase subunit 1 (COI) gene (Clare, Lim, Engstrom, Eger, & Hebert, 2007); all haemoplasma‐negative DNA samples gave positive signal in the mitochondrial 16S rRNA gene‐ and/or the COI‐specific PCR assays. All amplified sequences were subjected to chimeric sequence analysis using decipher (Wright, Yilmaz, & Noguera, 2012) and uchime (Edgar, Haas, Clemente, Quince, & Knight, 2011). All haemoplasma sequences have been deposited in GenBank under accession nos. MH245119–MH245194 and MK353807–MK353892; four positive samples were identified as Bartonella spp. during sequencing and were considered haemoplasma negative in our analyses.

Bat phylogenetic data

We used the rotl and ape packages in R to extract a bat phylogeny from the Open Tree of Life and to calculate branch lengths with Grafen's method (Michonneau, Brown, & Winter, 2016; Paradis, Claude, & Strimmer, 2004). To assess haemoplasma genotype sharing as a function of host phylogenetic similarity, we derived pairwise phylogenetic distances between the 33 sampled bat species (Figure S1). Because more evolutionarily distant and distinct species could display less frequent bacterial genotype sharing owing to ecological and immunological barriers to pathogen exposure and replication (Huang, Drake, Gittleman, & Altizer, 2015), we used the picante package and our bat phylogeny to derive evolutionary distinctiveness (Kembel et al., 2010).

Host species trait data

We obtained species‐level data on host traits relevant to pathogen transmission from previously published sources (Table S2). We obtained fecundity (litter size, litters per year), body mass and diet from the Amniote Life History and EltonTraits databases (Myhrvold et al., 2015; Wilman et al., 2014). For foraging ecology, which could affect bacterial exposure (e.g., trophic interactions; Kellner et al., 2018), we defined three dietary guilds: frugivory (including nectarivory; n = 11), insectivory (n = 18) and carnivory (including sanguivory and piscivory, n = 4; González‐Salazar, Martínez‐Meyer, & López‐Santiago, 2014). We also considered the proportion of plant‐based items in the diet. We simplified foraging strata into aerial (n = 14), arboreal (n = 16, including scansorial), and ground‐ or aquatic‐level foraging (n = 3). We also expanded prior compilations of wing aspect ratios and roost preferences to serve as proxies for ecological overlap among species (Fenton et al., 2001; Herrera et al., 2018; Reid, 1997). Roost type was simplified to open (e.g., only foliage; n = 6) or closed (e.g., hollows, caves; n = 27), and roost flexibility was simplified to using one (n = 16) or multiple roost types (n = 17). We classified maximum colony sizes as small‐to‐medium (e.g., <100 individuals; n = 20) or large (e.g., hundreds to thousands; n = 13; Reid, 1997; Santana, Dial, Eiting, & Alfaro, 2011), as most values were reported in ranges. We did not record pairwise sympatry (e.g., Luis et al., 2015; McKee et al., 2019) given that all species occur in Belize (Figure S2). Yet because more widely distributed species could have more opportunities for pathogen sharing due to range overlap, we used the geosphere package and data from the International Union for Conservation of Nature to derive geographical range size (Baillie, Hilton‐Taylor, & Stuart, 2004; Hijmans, Williams, Vennes, & Hijmans, 2019). Missing species‐level traits were taken from other databases, primary literature or closely related species (Table S2).

Individual‐level analysis of bat infection status

We first used the prevalence package to estimate haemoplasma infection prevalence and its 95% confidence interval (CI; Wilson method). We then used phylogenetic generalized linear mixed models (GLMMs) to test if infection status varied by sex, reproductive status, year, site and ectoparasite presence while accounting for bat phylogenetic relatedness. We fit candidate GLMMs that considered all fixed effects as well as interactions between sex and reproduction and between site and year. We also considered a model that excluded site to account for possible nonindependence of LAR and KK alongside an intercept‐only model. As vampire bats were banded for a mark–recapture study (Volokhov, Becker, et al., 2017) and some were sampled between and within years (n = 14), we randomly selected one of each recapture. After removing recaptures and missing values (n = 323), we fit the phylogenetic GLMMs using the brms package, default priors, and infection status as a Bernoulli‐distributed response. We included random effects for bat species and phylogeny, the latter of which used the phylogenetic covariance matrix (Bürkner, 2017). We ran four chains for 20,000 iterations with a burn‐in period of 10,000, thinned every 10 steps, for a total of 4,000 samples. We compared GLMMs using the leave‐one‐out cross‐validation (LOOIC) and assessed fit with a Bayesian R 2, including the total modelled variance and that attributed to only the fixed effects (Gelman, Goodrich, Gabry, & Vehtari, 2019; Vehtari, Gelman, & Gabry, 2017). We then estimated fixed effects (means and 95% highest density intervals [HDIs]) from the posterior distributions of each predictor from the top GLMM. Lastly, we assessed the sensitivity of individual‐level analyses to over‐representation of D. rotundus by randomly subsampling this species by the maximum n for other bat species.

Species‐level analysis of haemoplasma prevalence

We next calculated infection prevalence per species, using the metafor package to estimate logit‐transformed proportions and sampling variances (Viechtbauer, 2010). We used the nlme package to estimate phylogenetic signal as Pagel's λ with a weighted phylogenetic generalized least squares (PGLS) model to account for within‐species variance (Garamszegi, 2014). We next used a graph‐partitioning algorithm, phylogenetic factorization, to flexibly identify clades with significantly different prevalence estimates at various taxonomic depths. We used the taxize package to obtain a bat taxonomy from the National Center for Biotechnology Information (Chamberlain & Szöcs, 2013) and used the phylofactor package to partition prevalence as a Bernoulli‐distributed response in a GLM (Washburne et al., 2019). We determined the number of significant bat clades using Holm's sequentially rejective test with a 5% family‐wise error rate.

To identify species trait correlates of prevalence, we fit 11 PGLS models (weighted by sampling variance) with body mass, annual fecundity (litters per year × pups per litter), dietary guild, quantitative diet, foraging strata, aspect ratio, roost type, roost flexibility, colony size, geographical range size and evolutionary distinctiveness as predictors. We also fit PGLS models with only sample size or an intercept. We compared models with Akaike's information criterion corrected for small sample sizes (AICc) and estimated R 2 (Burnham & Anderson, 2002).

Hemoplasma phylogenetic analyses and genotype assignment

We compared our 16S rRNA gene sequences to those in GenBank (Volokhov, Becker, et al., 2017; Volokhov et al., 2011). Briefly, we aligned sequences using clustal x, and inter‐ and intraspecies similarity values were generated using bioedit. Genetic distances were calculated with the Kimura 2‐parameter and Tamura–Nei models, and the phylogeny was constructed using mega x with the minimum evolution algorithm (Kumar, Stecher, Li, Knyaz, & Tamura, 2018).

We assigned haemoplasma genotypes to positive bats based on analysis of the 16S rRNA partial gene (860–1,000 bp) sequences in GenBank and their clustering on the phylogeny. Genotypes were designated as novel if (a) sequences differed from the closest haemoplasma sequences in GenBank by ≥1.5% and/or (b) if sequence similarity was <1.5% but genotype‐specific reproducible mutations (at least two per sequence) were observed between haemoplasma sequences from at least two independent bat samples and the nearest GenBank haemoplasma sequences. These genotype‐specific mutations were further used to differentiate closely related haemoplasma genotypes from our sample. We caution that genotype is not synonymous with species, as analysis of the 16S rRNA gene alone is insufficient for accurate species identification of Mycoplasma spp. (Volokhov et al., 2012). Future studies using genomics or housekeeping genes may identify independent but closely related haemoplasma species in our genotypes.

To assess if haemoplasma genotype assignments were associated with site and year, we used χ2 tests with p values generated through a Monte Carlo procedure. Prior to our phylogenetic and network analyses of genotype distributions across bat species (see below), we used another χ2 test to assess the association between haemoplasma genotype identify and bat host identity.

Evolutionary relationships between bats and haemoplasmas

To determine the degree to which bat haemoplasma genotypes display host specificity and to describe their evolutionary relationships with host species, we used our bat and haemoplasma phylogenies to construct a binary association matrix. To test the dependence of the haemoplasma phylogeny upon the bat phylogeny and thus assess evidence of evolutionary codivergence, we applied the Procrustes Approach to Cophylogeny (PACo) using distance matrices and the paco package (Hutchinson, Cagua, Balbuena, Stouffer, & Poisot, 2017). We used a jackknife procedure to estimate the degree to which each bat–genotype link supported a hypothesis of phylogenetic congruence; links were supported if their upper 95% confidence interval was below the mean of all squared jackknife residuals (Balbuena, Míguez‐Lozano, & Blasco‐Costa, 2013).

Haemoplasma genotype sharing among bat species

We used haemoplasma genotype assignments to create a network, with each node representing a bat species and edges representing shared genotypes among bat species pairs. We built an adjacency matrix using the igraph package and used the Louvain method to assess the structure of bat–haemoplasma communities within this network (Csardi & Nepusz, 2006). To test whether the distribution of haemoplasma genotypes across our Neotropical bat species is shaped by host phylogeny, we used two GLMs to predict counts of shared genotypes (Poisson errors) and the presence of sharing (binomial errors) by phylogenetic distance between bat species. We assessed statistical significance with a quadratic assignment procedure via the sna package (Butts, 2008).

We calculated two metrics of network centrality to quantify different aspects of how important a node (bat species) is to haemoplasma genotype sharing: degree and eigenvector centrality (Bell, Atkinson, & Carlson, 1999). Whereas degree indicates the number of other species with which a host shares bacterial genotypes (i.e., links per node), eigenvector centrality indicates the tendency for a host to share genotypes with species that also share more genotypes (i.e., connectivity). Eigenvector centrality is thus an extension of degree that can identify hubs of parasite sharing (Gómez, Nunn, & Verdú, 2013). These two metrics were moderately correlated (ρ = 0.59), with many non‐zero degree species displaying zero eigenvector centrality. To examine spatial and temporal patterns in host centrality, we built separate adjacency networks per each site and year. We fit separate GLMs to determine how haemoplasma sharing centrality was predicted by site, year and the two‐way interaction. Degree was modelled as a Poisson‐distributed response, while eigenvector centrality was logit‐transformed and used Gaussian errors. We next applied phylogenetic factorization to both metrics and weighted the algorithms by the square‐root sample size per species (Garamszegi, 2014). We then fit the same PGLS models used in our prevalence analysis to identify the most competitive trait predictors of bat species centrality to haemoplasma sharing. Lastly, to assess whether network centrality is associated with haemoplasma prevalence, we fit two weighted PGLS models with each centrality metric as a univariate predictor.

RESULTS

Haemoplasma infection status

We detected sequence‐confirmed haemoplasma infection in 239 of 469 individuals (51%; 95% CI: 46%–55%), with positive individuals in 23 of the 33 sampled bat species (Table S3). The most parsimonious phylogenetic GLMM explained 22% of the modelled variance in infection status and included only sex, reproductive status, ectoparasites and year (ΔLOOIC = 0.16, w = .27; Table S4), suggesting site to be uninformative. Males had higher odds of infection than females (odds ratio [OR] = 2.35, 95% HDI: 1.36–4.03), but risk was unrelated to ectoparasitism, reproduction or year (Figure 2a). This result was weakly sensitive to over‐representation by Desmodus rotundus, as effect size and the prevalence sex difference was weaker when we randomly subsampled this species (Figure 2). Across all models, fixed effects only explained up to 7% of the modelled variance, suggesting more variation explained by the species and phylogeny random effects.

Figure 2

Predictors of individual bat haemoplasma infection status. (a) Odds ratios and 95% HDIs from the most parsimonious phylogenetic GLMM (Table S4). Estimates that do not overlap with 1 (dashed line) are displayed in black. Reference levels for the odds ratios include bats sampled at LAR, females, reproductive bats, absence of ectoparasites and bats sampled in 2017. (b) Infection prevalence and 95% confidence intervals (Wilson method) stratified by sex. Results are shown for the full data set and after randomly subsampling Desmodus rotundus

Interspecies variation in haemoplasma prevalence

Across bat species, haemoplasma prevalence ranged from 0% to 100% (mean = 37%). We estimated Pagel's λ in logit‐transformed prevalence to be .39, indicating moderate phylogenetic signal. Similarly, phylogenetic factorization identified one bat clade with significantly lower prevalence compared to the paraphyletic remainder: the Emballonuridae (12% infected; Figure 3a). Our trait‐based analysis showed that relatively larger species (β = 1.48, p = .01, R 2 = .24) and those with larger colonies (β = 0.67, p = .06, R 2 = .20) had higher prevalence (Figure 3b; Table 1). Relatively heavier (≥20 g) and larger colony species included D. rotundus, Molossus nigricans and Pteronotus mesoamericanus, for which prevalence was greater than 58%. Although these three species were also heavily sampled, other well‐sampled species such as Sturnira parvidens and Carollia sowelli had lower prevalence, and sample size did not predict prevalence (Table 1).

Figure 3

Predictors of species‐level haemoplasma prevalence across the Belize bat community. (a) Clades with significantly different prevalence are highlighted. (b) Results from the top PGLS models predicting prevalence as a function of mass and colony size. Model fit and 95% confidence intervals are shown overlaid with data scaled by sample size; species from the clade identified through phylogenetic factorization are coloured as in (a). Species identified by phylogenetic factorization (Saccopteryx bilineata and Rhynchonycteris naso) and with larger body mass, colony size and haemoplasma prevalence (Desmodus rotundus, Molossus nigricans and Pteronotus mesoamericanus) are shown to the right (photographs by Sherri and Brock Fenton) [Colour figure can be viewed at wileyonlinelibrary.com]

Table 1

Competing weighted phylogenetic generalized least squares models predicting haemoplasma infection prevalence (logit‐transformed) across the Belize bat community

Model structure	k	ΔAICc	wi	R 2
Log body mass	2	0.00	.40	.24
Maximum colony size	2	1.71	.17	.20
Roost type	2	2.82	.10	.18
Foraging strata	3	2.88	.09	.24
Roost flexibility	2	3.56	.07	.16
Per cent plants in diet	2	4.30	.05	.14
Dietary guild	3	5.46	.03	.17
Log evolutionary distinctiveness	2	5.50	.02	.11
Log aspect ratio	2	5.61	.02	.10
Log annual fecundity	2	6.13	.02	.09
Square‐root geographical range size	2	6.27	.02	.08
1 (intercept only)	1	6.77	.01	.00
Sample size	2	8.38	.01	.02

Models are ranked by ΔAICc with the number of coefficients (k), Akaike weights (w) and a likelihood ratio test pseudo‐R 2.

Hemoplasma genotype diversity

Our phylogenetic analysis identified 29 Mycoplasma genotypes in the Belize bat community (Table 2), including three previously identified from vampire bats (VBG1–3; Volokhov, Becker, et al., 2017). All genotypes demonstrated minor levels of sequence variability (99.5%–100%; Table 2). Based on comparisons with sequences from GenBank, 21 of these genotypes represent novel haemoplasmas, five are closely related to nonhaemoplasma mycoplasmas, and many are phylogenetically related to previously identified Mycoplasma spp. and haemoplasmas identified from other bat species (e.g., Di Cataldo et al., 2020; Millán et al., 2019; Millán, López‐Roig, Delicado, Serra‐Cobo, & Esperón, 2015), bat ticks (e.g., Hornok et al., 2019), primates including humans (e.g., Alcorn et al., 2020; Bonato, Figueiredo, Gonçalves, Machado, & André, 2015; Hattori et al., 2020) and rodents (e.g., Gonçalves et al., 2015; Goto, Yasuda, Hayashimoto, & Ebukuro, 2010; Vieira et al., 2009). A more detailed description of these 29 bacterial genotypes is provided in the Figure S3.

Table 2

Haemoplasma genotypes identified from the Belize bat community. Genotypes are given with their bat host species, representative GenBank numbers and intragenotype variability

Genotype	Host species	Representative GenBank number	Mean intragenotype sequence similarity (%)
VBG1	Desmodus rotundus, Pteronotus fulvus c	KY932701	99.8
VBG2	Desmodus rotundus	KY932678	99.9
VBG3	Desmodus rotundus	KY932722	99.6
CS1 b	Carollia sowelli	MK353833	100
CS2 b	Carollia sowelli, C. perspicillata	MH245134	99.7
MR1 b	Molossus nigricans	MH245174	99.7
MR2 b	Molossus nigricans	MH245151	NA a
PPM b	Pteronotus mesoamericanus, P. fulvus c	MH245159	99.9
EF1 b	Eptesicus furinalis, Saccopteryx bilineata c , Glossophaga soricina c	MH245147	99.6
EF2 b	Eptesicus furinalis	MH245131	99.9
NM b	Natalus mexicanus c	MK353818	NA a
LE b	Lophostoma evotis	MK353892	99.9
TC1 b	Trachops cirrhosus	MH245145	99.8
TC2 b	Trachops cirrhosus	MK353860	99.8
APH1 b	Dermanura phaeotis, D. watsoni, A. lituratus c	MH245132	100
APH2 b	Artibeus jamaicensis, A. lituratus, A. intermedius	MH245187	99.9
APH3 b	Artibeus intermedius	MH245186	99.8
GLS b	Glossophaga soricina	MK353874	99.5
MYE b	Myotis elegans, Myotis pilosatibialis c	MK353840	100
MYK b	Myotis pilosatibialis c	MH245153	NA a
UB b	Uroderma convexum	MK353869	99.8
PLU b	Platyrrhinus helleri c , Uroderma convexum c	MK353883	99.6
SP b	Sturnira parvidens; A. lituratus c	MH245168	99.4
RHN b	Rhynchonycteris naso c	MK353871	NA a
M. moatsii‐like 1 d	Pteronotus mesoamericanus c	MK353864	NA a
M. moatsii‐like 2 d	Myotis pilosatibialis c	MK353862	NA a
M. moatsii‐like 3 d	Rhynchonycteris naso c	MH245146	NA a
M. lagogenitalium‐like d	Glossophaga soricina c	MH245140	NA a
M. muris‐like d	Saccopteryx bilineata c	MH245138	NA a

Intragenotype sequence variability could not be assessed, as only one sequence was identified.

Novel haemoplasma genotypes.

Genotypes were detected in only one individual of these bat species.

Non‐haemoplasma Mycoplasma genotypes.

After controlling for multiple comparisons, our 29 bacterial genotypes were associated with site (χ2 = 47.11, p < .01) and year (χ2 = 40.40, p < .01). Genotype composition was more diverse at LAR (Figure S4), and KK haemoplasmas were dominated by vampire bat genotypes (VBG1–3). Genotype composition was more idiosyncratic by study year. However, these 29 bacterial genotypes were most strongly associated with bat species (χ2 = 3,532, p < .01; Figure S4).

Bat–haemoplasma evolutionary relationships

Although some haemoplasma genotypes were shared between bat species (i.e., VBG1, CS2, PPM, EF1, AH1–2, MYE, PLU, SP; n = 9), most showed strong host specificity (n = 20; Table 2). Our coevolutionary analysis (PACo) supported congruence between the bat and haemoplasma phylogenies (  = 36.21, p = .02, n = 1,000; Figure 4), suggesting that haemoplasma evolution has mostly tracked bat speciation. However, PACo also demonstrated that only 56% of the 41 unique bat–haemoplasma links displayed significant evidence of coevolution (Figure S5), and these patterns were almost exclusively found within the Phyllostomidae (with the exception of Saccopteryx bilineata and its Mycoplasma muris‐like bacterial genotype). The other 18 bat–haemoplasma links therefore displayed evidence of phylogenetic incongruence and thus probable host shifts.

Figure 4

Evolutionary relationships between Belize bats and haemoplasma genotypes. The cophylogeny plot shows the bat phylogeny on the left and the haemoplasma genotype phylogeny on the right. We used the treespace package to collapse our complete haemoplasma phylogeny (Figure S3) to only the 29 bacterial genotypes (Jombart, Kendall, Almagro‐Garcia, & Colijn, 2017). Lines display bat–haemoplasma associations and are shaded by the inverse of the squared residuals from PACo (i.e., dark lines show small residuals more indicative of coevolution)

Hemoplasma genotype sharing networks

Within our bat–haemoplasma network, genotype sharing was restricted to five host communities, whereas six genotypes were each restricted to a single bat species (Figure 5a). GLMs showed that both the frequency and the presence of genotype sharing declined with phylogenetic distance between bat species (Poisson: p < .001, R 2 = .08; binomial: p < .001, R 2 = .51; Figure 5b).

Figure 5

Patterns of haemoplasma genotype sharing across the Belize bat community. (a) Nodes in the genotype network represent bat species (abbreviated by Latin binomials), and edges represent a shared genotype. Nodes are coloured by communities identified with the Louvain method and are scaled by the number of individuals per species. (b) Matrix showing pairwise haemoplasma genotype sharing, coloured by the number of genotypes shared between bat species [Colour figure can be viewed at wileyonlinelibrary.com]

Bats shared haemoplasma genotypes with zero to five other species (i.e., degree), and most hosts had eigenvector centrality of zero (Figure S6). Six bat species had nonzero eigenvector centrality (37%–100%), indicating that they shared more haemoplasma genotypes with other highly connected host species. Stratifying our haemoplasma network across sites and years showed that centrality measures varied by space but not time (Figure S7; Table S5). We observed no haemoplasma sharing at KK, which probably reflects lower host diversity (Herrera et al., 2018).

Phylogenetic factorization identified similar bat clades with significantly different centrality compared to the paraphyletic remainder (Figure 6a, b). For degree, the algorithm only identified Artibeus lituratus as being more central (  = 5) than other bats (  = 1.14). However, phylogenetic factorization identified three taxa in the subfamily Stenodermatinae that had significantly elevated eigenvector centrality: the genera Artibeus and Dermanura (  = .67 compared to  = .02 in all other bats), the species A. lituratus (  = 1 compared to  = .12), and the species S. parvidens (  = .37 compared to  = .15). Mirroring these results, phylogenetic signal was absent for degree (λ = 0) but high for eigenvector centrality (λ = .93).

Figure 6

Phylogenetic patterns in haemoplasma genotype networks for Belize bat species (a) degree and (b) eigenvector centrality. Clades showing significantly different centrality metrics are highlighted, and points are scaled by observed values. Results from the top PGLS models predicting both centrality metrics as a function of bat species traits (c and d). Model fit and 95% confidence intervals are shown overlaid with data scaled by sample size; species from the clades identified through phylogenetic factorization are coloured as in (a) and (b) [Colour figure can be viewed at wileyonlinelibrary.com]

Trait‐based analyses showed that degree centrality was best predicted by diet (Table S6); bat species feeding more heavily on fruit and nectar shared more bacterial genotypes with other species (β = 0.004, p < .001, R 2 = .20; Figure 6c). Similarly, eigenvector centrality was best predicted by bat colony size and diet (Table S7); highly central species had small colonies (β large = –1.93, p = .05, R 2 = .13) and fed more on plants (β = 0.03, p < .01, R 2 = .10; Figure 6d).

As a final analysis, we assessed whether network centrality (i.e., a bat species’ role in haemoplasma genotype sharing) predicted contemporary infection prevalence (Figure S8). However, we found no associations between species‐level infection prevalence and centrality as measured by degree (β = –0.13, R 2 = .03, p = .42) or eigenvector centrality (β = –0.20, R 2 < .01, p = .79).

DISCUSSION

By examining the prevalence and distribution of a common bacterial pathogen (haemoplasmas) in a diverse bat community, we expanded analysis of the ecological and evolutionary predictors of bat infection and pathogen sharing beyond viruses. Across the bat community, haemoplasma infection risk was weakly higher for males but was better predicted by phylogeny, with large‐bodied and large‐colony bat species showing greater prevalence. Haemoplasmas showed high diversity and mostly strict host associations, with strong congruence between the bat and haemoplasma phylogenies. Although codivergence was supported by our analyses, we also observed haemoplasma genotype sharing and evidence of historical host shifts between closely related bats. Species most central to this haemoplasma sharing network displayed taxonomic clustering and were disproportionately frugivores and nectarivores. Yet these highly central bat species did not also have the highest haemoplasma prevalence, reinforcing mostly infrequent bacterial sharing between species. Our work reveals phylogenetic patterns in haemoplasma infection in a diverse bat community while contributing to broader efforts to understand the host specificity of bacterial pathogens and their cross‐species transmission patterns in wildlife.

Whereas many bacterial pathogens, including haemoplasmas, are common in bats (Bai et al., 2011; Becker et al., 2018; Ikeda et al., 2017; Mascarelli et al., 2014; Millán et al., 2015; Volokhov, Becker, et al., 2017), the factors that confer high infection probability are poorly understood. In the Belize bat community, the odds of haemoplasma infection tended to be higher in males, suggesting male‐biased transmission as detected in feline and canine systems (Soto et al., 2017; Walker Vergara et al., 2016). Such patterns could stem from males mounting weaker immune responses than females (Kelly, Stoehr, Nunn, Smyth, & Prokop, 2018) or to male defence of multifemale roosts in many Neotropical bat species (Voigt, von Helversen, Michener, & Kunz, 2001). Direct transmission of haemoplasmas has been demonstrated in feline and rodent systems (Cohen et al., 2018; Museux et al., 2009) but only inferred in bats from metagenomic studies detecting these bacteria in saliva (Volokhov, Becker, et al., 2017). We found weak support for the hypothesis that ectoparasites play a role in infection risks (Hornok et al., 2019; Willi, Boretti, Meli, et al., 2007). The weak male bias in infection could cast further doubt on vector‐borne transmission, as females in some bat species have elevated ectoparasitism (Frank, Mendenhall, Judson, Daily, & Hadly, 2016); however, a secondary GLMM testing this hypothesis found generally weak support for a female bias in ectoparasitism in our system (Figure S9). Future work assessing ectoparasite burdens could better elucidate the roles of vectors in haemoplasma risk.

Across Neotropical bats sampled in Belize, we found phylogeny to be a better predictor of haemoplasma risk than individual traits, site or year. Phylogenetic factorization identified one clade, the Emballonuridae, with significantly lower prevalence than all other bats in the community. This moderate phylogenetic signal mirrors comparable effects of phylogeny for bat viruses (Guy, Thiagavel, Mideo, & Ratcliffe, 2019), similarly suggesting potential for innate differences in species susceptibility or pathogen exposure. Trait‐based analyses revealed that this taxonomic pattern was driven by heavier and large‐colony species having greater haemoplasma prevalence. Small‐bodied species could have low prevalence due to small blood volumes and low bacterial titres (Volokhov, Becker, et al., 2017). Alternatively, the positive, saturating relationship between body mass and bacterial prevalence could be driven by allometric patterns in competence (Downs et al., 2019), in contrast to weak or opposite relationships between mass and viral richness across bats (Guy et al., 2019; Han, Schmidt, et al., 2016). As larger‐bodied bat species can also be more abundant in Neotropical habitat fragments (Herrera et al., 2018), these results suggest land conversion could increase the frequency of bat species most capable of maintaining haemoplasma infection. Similarly, positive relationships between colony size and prevalence could support density‐dependent transmission of bacteria (McCallum et al., 2001), whereas mixed support has been found for bat viruses (Streicker et al., 2012; Webber, Fletcher, & Willis, 2017). Future work could test how community‐wide infection patterns vary across broader habitat gradients and use multiple bacteria to assess the generality of these trends.

Approximately two‐thirds of the Neotropical bat species sampled in Belize were infected by haemoplasmas, for which we observed high genetic diversity consistent with other studies of this pathogen in bats (Mascarelli et al., 2014; Millán et al., 2015; Volokhov, Becker, et al., 2017). However, these bacterial genotypes were mostly novel and only weakly related to haemoplasmas described elsewhere in Latin America (Ikeda et al., 2017; Millán et al., 2019), with the exception of those previously identified from vampire bats (Volokhov, Becker, et al., 2017). When considering the phylogenetic scale of genotypes, most haemoplasmas were host‐specific. Over half of our haemoplasma communities consisted of a single bat–genotype association, matching the degree of host specificity observed more generally for Mycoplasma spp. (Citti & Blanchard, 2013; Pitcher & Nicholas, 2005). When we did detect genotype sharing between species, this occurred mostly between closely related hosts (e.g., genotype PPM was detected in Pteronotus mesoamericanus and P. fulvus), indicating bat phylogenetic distance decreased the probability of bacterial transfer.

Analyses to characterize species centrality to the haemoplasma genotype sharing network showed that one species (Artibeus lituratus) and the subfamily Stenodermatinae played key roles. This clade, and especially the genera Artibeus and Dermanura (formerly all classified in Artibeus), was the only taxon with nonzero connectivity, and this pattern was reflected in fruit‐ and nectar‐based diets and small colonies being the primary predictors of centrality. The strictly frugivorous Stenodermatinae represents a recent divergence in the Phyllostomidae (Botero‐Castro et al., 2013), and high centrality of these species may indicate weaker phylogenetic barriers for bacterial transmission between hosts in this clade. Two other analyses reinforced infrequent and conserved haemoplasma sharing between species. First, phylogenetic patterns in prevalence were distinct from those in genotype sharing centrality (e.g., large‐colony species had higher prevalence but lower connectivity), and prevalence accordingly did not predict centrality. Second, we found general congruence between bat and haemoplasma phylogenies. Although this shows codivergence is a strong evolutionary force, congruence can also stem from preferential jumps to closely related hosts (De Vienne et al., 2013). Although we cannot rule out that some host shifts may be artefacts of the limited resolution of both the phylogenies, our evolutionary analyses and genotype sharing results imply that haemoplasma host shifts are possible yet rare.

By sampling a diverse assemblage of bacterial genotypes in an ecologically and evolutionary rich host community, our work has broader implications for our understanding of disease emergence. Many bacterial pathogens are thought to be generalists and relatively unlikely to specialize in a novel host (Pedersen, Altizer, Poss, Cunningham, & Nunn, 2005; Woolhouse & Gowtage‐Sequeria, 2005), in contrast to many viruses in which host shifts are more common owing to high mutation rates and short infectious periods (Geoghegan et al., 2017; Longdon, Brockhurst, Russell, Welch, & Jiggins, 2014). Recent theoretical work suggests host shift speciation may be less common for bacteria because of higher phenotypic plasticity (e.g., the ability to reside in diverse habitats) and a slower tempo of evolution (Bonneaud et al., 2019). Obligate reliance of Mycoplasma spp. on host cells and more chronic infections probably explains their propensity to specialize (Citti & Blanchard, 2013; Cohen et al., 2018). More broadly, however, using genetics to infer pathogen sharing, rather than coarser phylogenetic scales (e.g., species complexes or genera), is increasingly showing that many bacterial strains may be more host‐specific (Withenshaw, Devevey, Pedersen, & Fenton, 2016). The high specialism of bat haemoplasma genotypes thus underlines the importance of using finer phylogenetic scales in the study of infectious disease (Fountain‐Jones et al., 2018; Graham, Storch, & Machac, 2018).

Comparative analyses of viruses have suggested that phylogenetically conserved pathogen jumps between species may be a broader generality in the study of disease emergence (Albery, Eskew, Ross, & Olival, 2019; Luis et al., 2015; Streicker et al., 2010). With few exceptions, our results on haemoplasma genotype sharing between Neotropical bat species are generally consistent with this pattern for a bacterial pathogen. Two cases in which haemoplasmas were shared between more distantly related species included the VBG1 genotype in Desmodus rotundus and Pteronotus fulvus (Phyllostomidae and Mormoopidae) and the EF1 genotype in Glossophaga soricina and Saccopteryx bilineata (Phyllostomidae and Emballonuridae). For the latter, both bat species co‐roost in the LAR, which suggests an ecological context for pathogen exposure over current timescales. However, other genetic markers (e.g., rpoB, rpoC, gyrB) would be necessary to infer contemporary cross‐species transmission (Kämpfer & Glaeser, 2012; Volokhov et al., 2012), as analysis of the 16S rRNA gene alone is insufficient for haemoplasma species identification (Volokhov et al., 2012). If haemoplasmas are more likely to specialize rather than expand their range into new and unrelated species, genotype sharing between unrelated bats could represent more transient spillovers (Bonneaud et al., 2019). As specialized pathogens could be more transmissible than generalists (Garamszegi, 2006), species with high infection prevalence of specialist genotypes could be prioritized for bacterial surveillance.

In conclusion, our analysis of a diverse community of bats and their pathogen genotypes identifies several key ecological and evolutionary factors structuring bacterial infection within and between species and provides a starting point for contrasts with such patterns for viruses. Similar to bat viruses, we found moderate phylogenetic signal in haemoplasma prevalence. However, these phylogenetic patterns in prevalence were decoupled from those describing bat species centrality in sharing haemoplasmas, such that genotype sharing was generally restricted by bat phylogeny. These findings imply codivergence of bats and their bacterial pathogens alongside rare and phylogenetically constrained host shifts. Future work more broadly characterizing the ecological and evolutionary determinants of bacterial infections in diverse host communities will improve our understanding of cross‐species transmission beyond viruses and contribute to efforts to understand the epidemiological consequences of bacterial pathogens.

AUTHOR CONTRIBUTIONS

D.J.B., K.A.S. and A.M.B. collected samples; N.B.S. and M.B.F. coordinated fieldwork; R.K.P., S.A. and D.G.S. supported laboratory analyses; D.V.V. and V.E.C. conducted molecular and phylogenetic analyses; and A.D.W. and D.G.S. assisted with statistical analyses. D.J.B. analysed data, produced figures and wrote the manuscript. All authors provided critical review of the manuscript.

CONFLICT OF INTEREST

We have no competing interests.

Supplementary Material

Click here for additional data file.