Original Research ARTICLE
Host Phylogeny, Geographic Overlap, and Roost Sharing Shape Parasite Communities in European Bats
- 1Graduate Degree Program in Ecology, Colorado State University, Fort Collins, CO, United States
- 2Department of Biology, Colorado State University, Fort Collins, CO, United States
- 3Division of Vector-Borne Diseases, Centers for Disease Control and Prevention, Fort Collins, CO, United States
- 4Centre for Zoonoses and Environmental Microbiology, Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
- 5Department of Parasitology and Parasitic Diseases, University of Agricultural Sciences and Veterinary Medicine Cluj-Napoca, Cluj-Napoca, Romania
- 6Department of Zoology, Hungarian Natural History Museum, Budapest, Hungary
- 7Department of Evolutionary Zoology and Human Biology, University of Debrecen, Debrecen, Hungary
- 8Department of Research and Collections, Natural History Museum, University of Oslo, Oslo, Norway
- 9Department of Parasitology and Zoology, University of Veterinary Medicine, Budapest, Hungary
- 10Evolutionary Systems Research Group, Centre for Ecological Research, Hungarian Academy of Sciences, Tihany, Hungary
- 11Terrestrial Ecology Unit, Department of Biology, Ghent University, Ghent, Belgium
- 12Department of Animal Ecology and Physiology, Radboud University, Nijmegen, Netherlands
- 13Endless Forms Group, Naturalis Biodiversity Center, Leiden, Netherlands
How multitrophic relationships between wildlife communities and their ectoparasitic vectors interact to shape the diversity of vector-borne microorganisms is poorly understood. Nested levels of dependence among microbes, vectors, and vertebrate hosts may have complicated effects on both microbial community assembly and evolution. We examined Bartonella sequences from European bats and their ectoparasites with a combination of network analysis, Bayesian phylogenetics, tip-association and cophylogeny tests, and linear regression to understand the ecological and evolutionary processes that shape parasite communities. We detected seven bat–ectoparasite–Bartonella communities that can be differentiated based on bat families and roosting patterns. Tips of the Bartonella tree were significantly clustered by host taxonomy and geography. We also found significant evidence of evolutionary congruence between bat host and Bartonella phylogenies, indicating that bacterial species have evolved to infect related bat species. Exploring these ecological and evolutionary associations further, we found that sharing of Bartonella species among bat hosts was strongly associated with host phylogenetic distance and roost sharing and less strongly with geographic range overlap. Ectoparasite sharing between hosts was strongly predicted by host phylogenetic distance, roost sharing, and geographic overlap but had no additive effect on Bartonella sharing. Finally, historical Bartonella host-switching was more frequent for closely related bats after accounting for sampling bias among bat species. This study helps to disentangle the complex ecology and evolution of Bartonella bacteria in bat species and their arthropod vectors. Our work provides insight into the important mechanisms that partition parasite communities among hosts, particularly the effect of host phylogeny and roost sharing, and could help to elucidate the evolutionary patterns of other diverse vector-borne microorganisms.
The enormous complexity of natural communities results from the large number of coexisting species and the diverse and unequal strengths of their interactions. Parasites, including macroparasites (e.g., worms and arthropods) and microparasites (e.g., bacteria and viruses), are an integral component of natural communities. Parasitism is a widespread life history strategy used by approximately one-third to over one-half of all species (Poulin, 2014; Morand, 2015). These parasitic organisms are under selective pressure to optimize their life history traits to efficiently colonize and reproduce in or on their hosts. This process of developing host specificity can be complicated, however, when there are several layers of parasitism. Such is the case for vector-borne microorganisms. For these organisms, selection can occur in both the host and the vector. When combined with host-associated selection on vectors, vector-borne microbes may exhibit complicated patterns of host and vector associations and phylogenetic differentiation. In light of this, classic models of parasite cospeciation and host-switching (de Vienne et al., 2013) must give way to novel approaches that examine the contributions of both hosts and vectors to the evolution and community assembly of vector-borne microorganisms.
Examining these processes is fundamental to understanding microbial diversity and surveillance of vector-borne microorganisms (Braks et al., 2011). Vectors vary in their host specificity, potentially leading to transmission of microorganisms to atypical hosts including humans. Vector-borne microorganisms account for a substantial proportion of emerging infectious diseases worldwide (Jones et al., 2008), and the zoonotic potential of mammalian viruses has been positively linked with being vector-borne (Olival et al., 2017). Knowledge of associations between microbes, hosts, and vectors will help to understand how humans become exposed to zoonotic agents and mitigate these risks. Thus, disentangling the ecological and evolutionary relationships between microorganisms and their hosts and vectors is important for managing the potential spillover of zoonotic agents to humans.
Disentangling these complex ecological and evolutionary processes requires sampling and analytical methods that integrate across trophic levels. If sampling is done across multiple ectoparasites and hosts, we can characterize the strength of host–parasite associations and identify host–vector–microbe communities using network-based approaches. Knowledge of these communities would directly facilitate disease management and the prevention of spillover events. For instance, hosts or vectors that have high infection prevalence or are connected with a large number of other nodes in the tripartite host–vector–microbe network may be targeted for pathogen surveillance or vector control. Looking at patterns of microbe sharing among hosts, we can highlight factors that constrain microbial host range using multiple regression, including host phylogenetic distance, vector sharing, geographic range overlap, and roost sharing as covariates (Streicker et al., 2010; Willoughby et al., 2017). Finally, we can examine how biases in historical microbial host-switching result in the observed congruence between host and microbial phylogenies (Charleston and Robertson, 2002).
To understand how complex host–vector–microbe communities are assembled and maintained in nature, we examined the associations of Bartonella spp. bacteria and ectoparasitic arthropods with their bat hosts using compiled data from nine European countries. We argue that Bartonella infections in bats and their ectoparasites represent an ideal system for understanding these complexities, first because Bartonella infections are prevalent and genetically diverse in many bat species studied to date (McKee et al., 2016; Stuckey et al., 2017b), providing rich data with which to analyze complex patterns. Second, bats are present on all continents except Antarctica and have traits that favor parasite transmission and geographic spread, including flight and long life spans. Many bat species are highly social and may form large colonies (Kerth, 2008), frequently co-roosting with other species, which could facilitate cross-species parasite transmission. Third, bats are a phylogenetically ancient lineage (Shi and Rabosky, 2015; Foley et al., 2016), allowing extended time for microbes and ectoparasites to develop host specificity. Finally, bats have many ectoparasites that vary in host specificity, ranging from highly specific wing mites (Bruyndonckx et al., 2009) to more generalist vectors like ticks (Hornok et al., 2016, 2017), which could have opposing effects on the evolution of host specificity in microorganisms they transmit. Such ectoparasite life history traits can interact with bat social systems in shaping microbial transmission (van Schaik et al., 2015). All these forces may combine to generate complex host–vector–microbe communities over evolutionary time but may be predictable given sufficient data and appropriate analytical methods. Moreover, bats are a highly threatened group of wildlife species, play central roles in ecosystems, and deliver valuable ecosystem services such as pollination and pest control (Boyles et al., 2011; Kunz et al., 2011); thus, ecological and evolutionary information on parasites could be informative for bat conservation and ecosystem sustainability (Whiteman and Parker, 2005; van Schaik et al., 2018). In addition to these factors, there are outstanding questions regarding the forces that drive Bartonella evolution in bats. Previous work has shown that the phylogeny of bat-associated Bartonella lineages is congruent with the bat host phylogeny (Lei and Olival, 2014), and Bartonella lineages tend to cluster by bat suborders and families (McKee et al., 2016). This previous work indicates that Bartonella species have developed some level of host specificity; however, the relative influence of ectoparasites and the biogeography of bat hosts on bat–Bartonella associations remain unclear.
Lastly, recent studies have highlighted the zoonotic potential of bat-associated Bartonella species. Bartonella spp. infections in humans and domestic animals can lead to symptoms ranging from mild fever to potentially life-threatening endocarditis (Chomel and Kasten, 2010). In one case of human endocarditis, the etiological agent was identified as a novel pathogenic species, Candidatus Bartonella mayotimonensis (Lin et al., 2010). This species and related strains have been identified in European and North American bat species (Veikkolainen et al., 2014; Lilley et al., 2017; Stuckey et al., 2017a; Urushadze et al., 2017). Bartonella spp. have been detected in numerous bat ectoparasites (Stuckey et al., 2017b; Hornok et al., 2019), some of which are known to occasionally attack humans (Jaenson et al., 1994; Estrada-Peña and Jongejan, 1999). A recent study also found serological evidence of a Bartonella species specific to fruit bats in humans in Nigeria where members of the community capture and sometimes eat bats (Bai et al., 2018). Given these emerging patterns, knowledge of the host and vector associations of bat-associated Bartonella species could have implications for managing spillover risk.
Our strategy to investigate how bat–ectoparasite–Bartonella communities are assembled and how they evolve involves a multifaceted analytical approach (Figure 1) that splits this problem into three fronts: (1) assessing the diversity of Bartonella species in European bats and their ectoparasites and the structure of bat–ectoparasite–Bartonella communities using network analysis, (2) understanding the evolutionary implications of these ecological patterns using tip-association and cophylogeny tests and phylogenetic measures of historical host-switching rates, and (3) linking patterns of Bartonella host specificity to ecological and evolutionary covariates using linear regression. We hypothesized that associations between bats, ectoparasites, and bacteria can be resolved into identifiable communities that separate by bat phylogeny, geographic overlap, and roosting patterns. Second, we expected that the phylogeny of Bartonella species will exhibit significant clustering by bat taxonomy and will have significant congruence with the phylogeny of bat species. Linking these patterns together, we hypothesized that host phylogenetic distance, ectoparasite sharing, geographic range overlap, and summer roost sharing are predictors of bacterial species assemblages and host-switching rates among bat species. This multifaceted approach aims to bridge ecological processes to observed evolutionary patterns to better understand the diversity and epizootiology of bartonellae in bats. Such an approach could be generalized to study and manage other microorganisms with complex, multihost dynamics.
Figure 1. Analysis framework diagram. Boxes link raw data (dark blue) to derived data (light blue), analytical methods (green), and results (orange) to connect all of the scientific analyses used in the study.
Materials and Methods
Study Sites and Specimen Collection
Bat ectoparasites were collected in the Netherlands, Belgium, Hungary, and Romania between 1993 and 2015. Sampling sites included roosting, swarming, and foraging areas, and all sampling occurred during the summer maternity and autumn mating phases when ectoparasites are more active (van Schaik and Kerth, 2017). In the Netherlands and Belgium, ectoparasite specimens were collected with forceps either directly from bats during inspections of bat boxes and night mist netting or from their roosts. In addition, ectoparasites were sampled during inspection of dead bats collected in the Netherlands between 1993 and 2011 and stored in the Naturalis Biodiversity Center, Leiden, the Netherlands. Bat flies from Hungary and Romania derived from a study by Sándor et al. (2018). All bats were morphologically identified to the species level. Initial identification of bat flies was based on morphological characteristics (Theodor and Moscona, 1954; Theodor, 1967). Ectoparasites were stored in 70% ethanol in separate vials prior to further analysis. The distribution of sampling sites is mapped in Figure 2, and the coordinates of sampling sites are listed in Table S13.
Figure 2. Map of ectoparasite sampling sites. Points summarize the total number of ectoparasite samples from the Netherlands, Belgium, Hungary, and Romania collected from each sampling site. Latitude, longitude, and ectoparasite counts for each sampling site are listed in Table S13.
Bat Species Tree, Geographic Range Overlap, and Roosting and Mating Data
A phylogenetic tree of bats (Figure S1) was obtained from the Open Tree of Life (http://www.opentreeoflife.org) from a previous study of bat taxonomy using multiple mitochondrial and nuclear loci (Shi and Rabosky, 2015). The tree was pruned to the 21 species in Table 1. Myotis oxygnathus was considered a synonym for My. blythii (Agnarsson et al., 2011; Balvín and Bartonička, 2014; Wilson and Reeder, 2015).
Table 1. Summary of Bartonella phylogenetic diversity for 21 European bat species included in our study.
The geographic ranges of each bat species were downloaded from the International Union for Conservation of Nature (IUCN) Red List website (http://www.iucnredlist.org) (IUCN, 2014). IUCN ranges are convenient data that are available for all the bat species in this study, and previous studies have successfully used these data for understanding the determinants of viral diversity in bats (Luis et al., 2013, 2015; Maganga et al., 2014; Webber et al., 2017; Willoughby et al., 2017). Shape files were imported into R using the “readShapeSpatial” function in the “maptools” package (Bivand et al., 2017; R Core Team, 2018). Individual range maps (Figure S2) and a map of overlapping ranges (Figure 3) were generated by drawing shape files over the “worldHires” map from the “maps” package (Becker et al., 2016). We then calculated pairwise percent geographic range overlap between each bat species as described previously (McKee et al., 2016); see the Supplementary Material for more details. Data on roosting patterns and mating systems of bats (Tables S1 and S2) were collected from books by Dietz and Kiefer (2014) and Niethammer and Krapp (2001,2004).
Figure 3. Map of geographic ranges for the 21 European bat species studied. Transparent layers were mapped on top of one another to highlight regions with dense range overlap, primarily in Europe and the Middle East. Some species have additional overlap in Central and East Asia. Individual maps are shown in Figure S2.
Ectoparasite DNA Extraction and Barcoding
DNA from bat flies, mites, fleas, and bat bugs was extracted with ammonium hydroxide as described previously (Wielinga et al., 2006). DNA from bat ticks was extracted using the Qiagen DNeasy Blood & Tissue Kit according to the manufacturer's protocol for the purification of total DNA from ticks (Qiagen, Venlo, the Netherlands). Confirmation of ectoparasite identification was performed by sequencing a 658-base-pair fragment of the cytochrome oxidase subunit I (COI) using primers LCO1490 and HCO2198 (Folmer et al., 1994); see the Supplementary Material for details. For species identification, both strands of PCR products were Sanger sequenced (BaseClear, Leiden, the Netherlands) using the same forward and reverse primers as in the conventional PCR. Trimming and manual cleaning of COI sequences were performed in BioNumerics v7.1 (Applied Math, Belgium). A phylogeny was inferred using the GTR+Γ model with 25 distinct rate categories with 1,000 bootstrap replicates using RAxML (Stamatakis, 2014). Based on this phylogeny, individual associations between a Bartonella sequence and a vector species were corrected if the phylogenetic position of the COI sequence conflicted with the morphological identification, replacing the morphological identification with the phylogenetic identification.
Review of Ectoparasite Host Range
We recognize that the relationships between bat hosts and ectoparasites in our dataset may not capture the full ectoparasite host range due to infrequent associations between ectoparasites and some bat hosts. To capture some of this additional variation, we performed a literature review (Table S3) of the host range of the 17 ectoparasite species for which we had associated Bartonella data. The search was implemented in Google Scholar, Web of Science, and GenBank using the ectoparasite species epithet in the search terms. Additional publications were obtained based on citations therein and from previous reviews of bat fly host associations (Szentiványi et al., 2016). This review gathered a total of 302 publications from 1835 to 2018, 212 of which yielded information relevant to the ectoparasite and bat species in the current study. We counted the number of publications for which an ectoparasite species was noted as occurring on a bat species, adding the current study toward each total only if an ectoparasite was sampled from a bat species in our collection. Bat–ectoparasite associations noted in the studies that did not record the full species epithet for both ectoparasite or host species were excluded.
Bartonella Amplification and Sequencing
Ectoparasites were tested individually for the presence of Bartonella spp. with a conventional PCR assay targeting the citrate synthase gene (gltA) using primers designed by Norman et al. (1995); the Supplementary Material contains additional details on Bartonella detection protocols. Previous studies have found gltA sequences to be sufficiently diverse to distinguish among Bartonella species and some subspecies (La Scola et al., 2003). Additionally, gltA is the most common marker for genotyping Bartonella species (Kosoy et al., 2018); hence, it is useful for comparing Bartonella diversity across studies. For Bartonella species identification, both strands of PCR products were Sanger sequenced (BaseClear, Leiden, the Netherlands) using the same forward and reverse primers as in the conventional PCR. To minimize cross-contamination and false-positive results, positive (pool of Bartonella-positive ticks; Tijsse-Klasen et al., 2011) and negative (water only) controls were included in each batch tested by PCR. Furthermore, DNA extraction, PCR mix preparation, sample addition, and PCR-product analysis were performed in assigned separate labs. Trimming and manual cleaning of Bartonella sequences were performed in BioNumerics v7.1 (Applied Math, Belgium) together with Bartonella reference sequences available in GenBank.
Sequences were further confirmed as Bartonella through Basic Local Alignment Search Tool (BLAST) searches (https://blast.ncbi.nlm.nih.gov/Blast.cgi). Based on this initial screening, some sequences were identified as originating from Bartonella species not associated with bats or from bacteria of other genera and were removed before further analysis (Supplementary Material). For the purposes of ecological and evolutionary analysis, we have assumed that a Bartonella strain amplified from an ectoparasite species may also be carried by the bat species on which that ectoparasite was found. We argue that even incidental ectoparasitism on an atypical bat host may lead to transmission of bacteria and is thus important for understanding available parasite host range. The host associations of ectoparasite-derived sequences were validated for the subset of Bartonella species that have been characterized from bats in previous studies (Supplementary Material).
To aid in phylogenetic inference and the delineation of novel Bartonella species, additional Bartonella gltA sequences amplified from bats and ectoparasites were compiled from published articles and records in GenBank. This search was implemented in Google Scholar, Web of Science, and GenBank. Initial BLAST screening also indicated that some Bartonella sequences had close similarity to Bartonella sequences found in humans (Podsiadly et al., 2010; Lin et al., 2012) and stray dogs (Bai et al., 2010); thus, representative sequences from these studies were also included. Details on the origin of Bartonella sequences are listed in Table S4. Sequences were trimmed to a common length of 337 base pairs and aligned using the local and accurate L-INS-i method in MAFFT v7.187 (Katoh and Standley, 2013). The sequences were inspected for gaps and misalignments and were removed if they contained obvious errors.
Phylogenetic Analysis of Bartonella Sequences
Two datasets were created for phylogenetic analyses. First, we compiled the full set of Bartonella sequences from bats and their ectoparasites from Table S4 plus sequences from the current study (full dataset, n = 754), using three sequences from Brucella spp. as the outgroup. Second, this full set was restricted to Bartonella sequences from bats and their ectoparasites in Europe (European bat dataset, n = 456). While ectoparasites collected from roosts and bat boxes were tested for Bartonella spp., the DNA sequences derived from these samples were only used in the phylogenetic analysis of the full dataset. Since they did not contain information on the host species, they were excluded from the European dataset and were thus not used to assess Bartonella diversity among bat species or in the other tests detailed below (network analysis, cophylogeny, and regression). This restricted set contained sequences from the 21 bat species represented in Figure S1 and the 17 ectoparasites represented in our literature review. Sequences from My. oxygnathus were combined with those from My. blythii. Due to the potential confounding factor of recombination in phylogenetic inference (Posada and Crandall, 2002), we performed the pairwise homoplasy index test (Bruen et al., 2005) in SplitsTree v4.13.1 (Huson, 2005). Tests using both the full and European bat datasets indicated no significant evidence that recombination affects our phylogenetic inference (P = 0.41 and P = 0.25, respectively).
We selected models derived from both datasets to determine the best sequence evolution, speciation, and codon partitioning models. Following this procedure, a phylogenetic tree was generated for the full dataset in BEAST using the GTR+Γ+I sequence evolution model and the birth–death speciation model with incomplete sampling and rate partitions for each codon position. See the Supplementary Material for additional details on model selection and phylogeny generation.
Branches of the phylogenetic tree were collapsed according to probable Bartonella species. Classification of bacterial species is challenging and may not conform well to species concepts developed for eukaryotes (Konstantinidis et al., 2006; Fraser et al., 2007). La Scola et al. (2003) proposed that Bartonella species could be distinguished if gltA sequences differed by >4% identity. Konstantinidis et al. (2006) advocated for a more stringent approach wherein bacterial species are demarcated by >5% difference in sequence identity, which corresponds well with another standard of bacterial species, 70% DNA–DNA hybridization. We chose to follow this more conservative approach and collapsed branches into operational taxonomic units (OTUs) based on ≤5% identity among sequences. Using these Bartonella OTUs, we can understand their ecology, specifically their host and ectoparasite associations, factors that can aid in demarcating Bartonella species or species complexes (Kosoy, 2010; Kosoy et al., 2012). Additional post hoc comparisons of OTUs with recently published Bartonella sequences from insectivorous bats in China and western, central, and eastern Europe were performed (Han et al., 2017; Stuckey et al., 2017a; Corduneanu et al., 2018). Phylogenetic diversity (PD) of Bartonella sequences from each bat species was assessed by the number of OTUs found in the species and (Faith, 1992) PD index based on branch lengths in 100 posterior samples of the gltA tree. Faith's PD was calculated in R using the “picante” package (Kembel et al., 2014).
Network Analysis and Community Detection
Weights for edges linking bat, ectoparasite, and Bartonella nodes were initially assigned based on the number of citations linking ectoparasite species to bat species in the literature review (Table S3) or the number of Bartonella gltA sequences for a given Bartonella OTU linked to a host bat or host ectoparasite (Table S10). To account for sampling intensity on edge weights, we adjusted bat–ectoparasite edge weights, wab, by dividing the number of citations linking ectoparasite b to bat a, nab, by the sum of the total unique publications surveyed for bat a, xa, and the total unique publications surveyed for ectoparasite b, xb; thus, wab = nab/(xa+xb). Similarly, we adjusted weights for edges linking Bartonella OTUs to hosts (either bat species or ectoparasite species), wcd, by dividing the number of gltA sequences linking OTU d to host c, ncd, by the sum of the total gltA sequences obtained from host c, yc, and the total gltA sequences obtained for OTU d, yd; thus, wcd = ncd/(yc+yd). Therefore, edge weights were constrained to be between 0 and 0.5, with an actual range of 0.00431 to 0.429.
We performed community detection on the tripartite network using three algorithms available in the R “igraph” package: the information map method (Rosvall and Bergstrom, 2008), the Louvain method (Blondel et al., 2008), and the spin glass method (Reichardt and Bornholdt, 2006). The purpose of using multiple algorithms was to account for some uncertainty in the identification of communities (see the Supplementary Material for details on the algorithms). Communities were visualized using the “HiveR” package in R (Hanson et al., 2016).
To identify bat, ectoparasite, and Bartonella species that might be highly influential in the network, we examined nodes that were highly connected in the tripartite network based on calculation of their weighted degree. A node's weighted degree represents the sum of the edge weights connecting a node to other nodes. Nodes were selected as influential if they were in the top 25th percentile of weighted degree. We then examined how the selected nodes were connected to other nodes within the communities detected by the community detection algorithms.
Cophylogeny and Tip-Association Tests for Bats and Bartonella
Clustering of traits among tips of the Bartonella tree (European bat dataset) was tested using the Bayesian Tip-association Significance Testing (BaTS) program (Parker et al., 2008). We assessed the clustering of bat taxonomic traits (species, genera, families, and suborders) and countries sampled. Significance of clustering for each trait was assessed by comparing the calculated association index (AI) and parsimony score (PS) for 100 posterior samples of the Bartonella tree against null distributions generated from 1,000 randomizations of traits to tips along each sampled Bartonella tree.
Phylogenetic trees of bat species and Bartonella OTUs were assessed for evidence of evolutionary codivergence using two algorithms: the Procrustean Approach to Cophylogeny (PACo; Balbuena et al., 2013) and the ParaFit method (Legendre et al., 2002). Trees were imported into R using the “ape” package (Paradis et al., 2004, 2016) and then rescaled to have a maximum branch length of one by dividing all branch lengths by the longest branch in each tree. A binary association matrix linking bat species to Bartonella OTUs was assembled based on Table S10. We tested the pattern of codivergence using ParaFit using the “ParaFit” function in the “ape” package with 999 permutations and stored the P-values for the contributions of individual linkages (known as ParaFitLink1 or F1 statistics). Codivergence was tested in PACo using the “paco” package in R (Balbuena et al., 2016) with 1,000 permutations. PACo residuals and mean jackknife contributions for individual linkages were stored to compare with results from ParaFit. We performed the tests on the maximum clade credibility Bartonella tree and 100 sampled posterior trees to assess the effect of phylogenetic uncertainty in the global fit tests.
Bayesian Prediction of Bartonella Host-Switching Rates
Historical rates of host-switching by Bartonella lineages among bat host species were predicted using Bayesian ancestral state reconstruction in BEAST. The host-switching rates discussed here represent a latent biological process (movement of agents between species) that was not observed but can be estimated based on observed host-switching events in the phylogenetic tree. Since we are interested primarily in the process, we chose to look at estimated host-switching rates rather than observed events. With the European bat Bartonella sequence dataset, we used the GTR+Γ+I sequence evolution model and the birth–death speciation model with incomplete sampling and rate partitions for each codon position (see the Supplementary Material for more details). Rates with Bayes factors (BF) > 3 after the stochastic search variable selection were considered well-supported (Lemey et al., 2009). A graph representing Bartonella transitions among bat species was drawn based on well-supported rates using the “arcdiagram” package in R (Sanchez, 2013).
Regression analyses centered around two primary response datasets: a dissimilarity matrix calculated from the counts of Bartonella OTUs found in each of the 21 bat species or their associated ectoparasites and the predicted host-switching rates from the ancestral state reconstruction analysis (Figure 1). Bartonella and ectoparasite dissimilarity were calculated using the binomial and Cao indices (Cao et al., 1997; Anderson and Millar, 2004), which can handle variable sample sizes and can calculate dissimilarity between species that have no shared parasites (Oksanen et al., 2015), and the Spearman rank correlation, which was subtracted from one to transform it into a dissimilarity measure. These same indices were used to measure ectoparasite dissimilarity. Other indices, specifically Pearson correlation and Bray–Curtis dissimilarity, were explored but were not chosen because they violated assumptions of normality in residuals. Only those host-switching rates with BF > 3 were kept for this analysis to be confident in the estimated rate values. Additional details on data selection for regression can be found in the Supplementary Material.
There were five primary predictors considered in the regression analyses: dissimilarity in ectoparasite sharing between bat species, phylogenetic distance between bat species, bat geographic range overlap, summer roost sharing of bats, and a vector of the least sampled species from the bat–Bartonella association matrix (Figure 1). The phylogenetic distance matrix from Shi and Rabosky (2015) was initially scaled in terms of branch ages (in millions of years); hence, we rescaled the branch lengths to be between zero and one by dividing all branch lengths by the maximum length. The summer roosting patterns of bats are a binary variable indicating whether or not bats share roosts during the summer months. The vector of least sampled species was selected from the bat–Bartonella association matrix as the minimum row sum for each species pair; this was then log-transformed. We only used the vector of least sampled species in the regression of Bartonella host-switching rates because there appeared to be a sampling bias in the predicted rates such that better sampled species tended to have higher median rates (Pearson's R = 0.62, t = 4.15, df = 28, P = 0.00028). Our dissimilarity measures did not appear to have such a bias.
Before performing regressions, we rescaled all the data and predictors to standard normal distributions with the exception of roost sharing, which was retained as a binary predictor. We performed separate model selection procedures on three global candidate model sets including a global model with all predictors and all subsets of the global model. The first global model set used ectoparasite dissimilarity, host phylogenetic distance, geographic range overlap, and roost sharing to predict Bartonella dissimilarity. The second set used host phylogenetic distance, geographic range overlap, and roost sharing to predict ectoparasite dissimilarity. The third set used the samples from the least sampled host species, ectoparasite dissimilarity, host phylogenetic distance, geographic range overlap, and roost sharing to predict Bartonella host-switching rates. For all models containing Bartonella or ectoparasite dissimilarity as data or predictors, we performed model selection based on regressions using all three dissimilarity indices (Spearman correlation, binomial, and Cao indices). For the host-switching models, we performed model selection on both median and mean host-switching rates.
Models were fit using linear regression with normally distributed errors. We selected models based on iterative testing of predictors in the full model and ranked them according to the Akaike information criterion with a correction for finite sample sizes (AICc) using the “dredge” function in the “MuMIn” package in R (Barton, 2016). We chose the model with the smallest AICc unless another model was less than two AICc away from the top model (Burnham and Anderson, 2004), in which case we chose the simplest model based on the principle of parsimony. For both the global and the top models, we recorded adjusted R2, inspected residual plots and quantile–quantile plots, and performed a Shapiro–Wilk test (Shapiro and Wilk, 1965) to confirm normality of residuals. We recorded standardized main effect coefficients, t statistics, F statistics, and associated P-values for regression parameters. To assess model fit, we performed k-fold cross-validation with 10-folds using the “cv.lm” function in the “DAAG” package in R (Maindonald and Braun, 2015). To assess the relative importance of parameters in models, we recorded adjusted partial R2 and relative importance values. Relative importance was calculated with the “calc.relimp” function using the (Lindeman et al., 1980)method in the R package “relimpo” (Groemping and Matthias, 2013). Bootstrap confidence intervals for relative importance values were estimated from 1,000 replicates. Due to potential nonindependence of comparisons between pairs of species, we also used Mantel tests (Mantel, 1967) to examine correlations between responses and predictors. We calculated the Pearson correlation and compared it to a null distribution generated from 999 random combinations of cells from the two matrices.
Ectoparasites (n = 903) were collected in the Netherlands and Belgium from 268 individual bats belonging to 11 species (E. serotinus, Ny. noctula, Pi. nathusii, Pi. pipistrellus, Pi. pygmaeus, Pl. auritus, My. bechsteinii, My. dasycneme, My. daubentonii, My. mystacinus, and My. nattereri). In addition, 170 nycteribiid flies from 169 individual bats (Mn. schreibersii, My. bechsteinii, My. blythii, My. capaccinii, My. daubentonii, My. myotis, My. nattereri, R. blasii, R. euryale, R. ferrumequinum, and R. mehelyi) derived from a study by Sándor et al. (2018) in Hungary and Romania were included.
A total of 1,073 ectoparasites were collected across the four countries (Figure 2). Morphological and molecular identification revealed 15 ectoparasite species from 7 families: two ticks, Argas vespertilionis and Ixodes ariadnae (Ixodida: Argasidae, Ixodidae); one bat bug, Cimex pipistrelli (Hemiptera: Cimicidae); one bat flea, Ischnopsyllus variabilis (Siphonaptera: Ischnopsyllidae); seven bat flies, Basilia nana, B. nattereri, Nycteribia kolenatii, Nb. schmidlii, Penicillidia conspicua, Pn. dufourii, and Phthiridium biarticulatum (Diptera: Hippoboscoidea: Nycteribiidae); and two bat mites, Spinturnix andegavinus and S. plecotina (Mesostigmata: Spinturnicidae). Additional mite specimens (Mesostigmata: Macronyssidae, Spinturnicidae) that could not be identified to the species level by morphology were delineated as two distinct taxa by COI sequences (Supplementary Material; Figure S3). Ectoparasite specimen counts collected in the Netherlands, Belgium, Hungary, and Romania from each bat species are summarized in Table S11. Table S12 records individual ectoparasite species identifications, Bartonella testing results, and sampling sites.
Phylogenetic Relationships Between Bartonella Sequences
In total, 412 gltA sequences were obtained from the 1,073 ectoparasites (38%) collected from European bats in Belgium, Hungary, the Netherlands, and Romania (Table S12). After filtering out sequences that were not Bartonella, these 316 sequences were combined with the 438 reference sequences listed in Table S4, resulting in the full dataset of 754 Bartonella gltA sequences used for phylogenetic analysis and delineation of Bartonella OTUs. The subset of 456 gltA sequences comprising the European bat dataset represents data from nine countries: Belgium, Finland, Georgia, Hungary, the Netherlands, Poland, Romania, Slovenia, and the United Kingdom. The sequences from Finland, Georgia, Poland, Slovenia, and the United Kingdom derive from past observational studies of Bartonella infections in bats and their ectoparasites (Concannon et al., 2005; Morse et al., 2012; Veikkolainen et al., 2014; Lilley et al., 2015; Urushadze et al., 2017; Szubert-Kruszynska et al., 2018).
Based on our demarcation of 5% sequence divergence for separating Bartonella species, we observed 49 monophyletic clusters of gltA sequences identified as OTUs and 39 individual gltA sequences that are distinct from these OTUs out of the 754 gltA sequences analyzed, resulting in an estimate of at least 88 distinct Bartonella species found in bats worldwide, 20 of which are found in European bats and ectoparasites (Figure S4). All OTUs had strong posterior support (posterior node probability > 0.9); however, support for nodes connecting OTUs into larger clades decreased significantly for deeper nodes. A general pattern of separation between Bartonella OTUs and sequences found in New World bats (colored green in Figure S4) and Old World or European bats (colored blue and red in Figure S4, respectively) was observed, as in a previous analysis by McKee et al. (2016).
Bartonella diversity varied across European bat species, with a range of one to nine OTUs and a range of 0.44–4.05 for Faith's PD for E. nilssonii and Mn. schreibersii, respectively (Table 1). Faith's PD per species was significantly positively correlated with the log number of gltA sequences obtained for that species (Pearson's R = 0.88, t = 8.04, df = 19, P < 0.0001) and with the number of Bartonella OTUs observed per species (Pearson's R = 0.88, t = 7.96, df = 19, P < 0.0001). After accounting for this significant sampling effect on Bartonella diversity, there was no significant correlation between Faith's PD and the number of ectoparasites associated with each bat species (t = 0.6, df = 18, P = 0.56) or the number of OTUs and ectoparasites (t = 0.17, df = 18, P = 0.87).
Network Analysis and Community Assignment
The three community detection algorithms consistently identified seven communities: Min/Myo, VespA–VespE, and Rhi (Figure 4 and Table S14). There were only minor inconsistencies in the community assignment of a few species. Two algorithms (information map and spin glass) lumped Pl. auritus and S. plecotina into a distinct community, but the Louvain algorithm placed them with community VespC. The Louvain and spin glass algorithms grouped I. ariadnae with community VespB, whereas the information map algorithm grouped this species with community VespD.
Figure 4. Communities of Bartonella operational taxonomic units (OTUs), ectoparasite species, and bat species. Tripartite networks were drawn using separate axes for Bartonella, bats, and ectoparasite nodes. Edges connecting nodes were drawn based on Tables S3 and S10, with edge weights adjusted for sampling intensity. Communities were identified by three detection algorithms (information map, Louvain, and spin glass). Species membership in each community is recorded in Table S14.
The seven communities are broadly organized based on host phylogeny, geographic overlap, and roost sharing (Figures S1, S2 and Tables S1, S2). Min/Myo contains species from two related families of bats (Vespertilionidae and Miniopteridae) that roost together in caves or other cave-like structures predominantly in southern Europe (with the exception of My. myotis, which is more widespread in Europe; Figure S2). Rhi contains only Rhinolophus species that roost together in caves mostly in southern Europe, North Africa, and the Middle East (Figure S2 and Table S2). VespA contains three closely related Myotis species that roost in tree cavities during the summer and swarm at underground sites during autumn and hibernate there in winter. VespB contains vespertilionid bats that roost in tree cavities, buildings, and caves during the summer and winter, although with little overlap in roosting patterns among species. However, these species do share some other traits in common, including long-distance migration (My. dasycneme, Ny. noctula, and Pi. nathusii) and a relatively northern distribution within Europe. VespC contains two Eptesicus species that roost in buildings during the summer and winter. Communities VespD and VespE contained single bat species that had highly specific ectoparasite or Bartonella species and were thus segregated from other communities despite having phylogenetic similarity or similar roosting habits to other species in these communities.
For 11 of the Bartonella OTUs in this study, community assignments corresponded well with other sequence data collected from related bats in Africa, Asia, and Europe (Kosoy et al., 2010a; Lin et al., 2012; Morse et al., 2012; Anh et al., 2015; Wilkinson et al., 2016; Han et al., 2017; Lilley et al., 2017; Stuckey et al., 2017a; Corduneanu et al., 2018). Further details on comparisons between OTUs and other sequences from previous studies can be found in the Supplementary Material.
Fifteen nodes were identified as being influential based on their weighted degree in the bat–ectoparasite–Bartonella association network. Five bats (Mn. schreibersii, My. blythii, My. daubentonii, My. myotis, and R. ferrumequinum), nine ectoparasites (A. vespertilionis, B. nana, C. pipistrelli, Nb. kolenatii, Nb. schmidlii, Pn. conspicua, Pn. dufourii, Ph. biarticulatum, and S. myoti), and one Bartonella OTU (OTU19) were identified as influential because they fell in the top 25th percentile for weighted degree. The ecology of these species may explain their influence in the network. My. daubentonii males are known to form social colonies and therefore high ectoparasite densities have been observed on both sexes (Encarnação et al., 2012). Mn. schreibersii, My. blythii, My. myotis, and R. ferrumequinum form mixed roosts in caves (Table S2). All nine ectoparasite species are very promiscuous in their host associations (Table S3). We note that C. pipistrelli and A. vespertilionis are known to bite humans (Jaenson et al., 1994; Estrada-Peña and Jongejan, 1999; Whyte et al., 2001) and were grouped in the same community as OTU19 and OTU26, which have both been found to infect humans (Veikkolainen et al., 2014; Urushadze et al., 2017). For additional details on the identification of highly influential nodes and their community assignments, see the Supplementary Material.
Cophylogeny and Tip-Association Tests for Bats and Bartonella
Tip-association tests showed significant clustering for all taxonomic levels (species, genera, families, and suborders) and sampled countries with none of the observed AI or PS distributions overlapping with the null distributions (Table S15). Mean AI and PS values were smaller for genus, family, and suborder compared to the AI and PS values for sampled countries, indicating that the clustering of tips of the tree is better explained by host phylogeny than geography. The cophylogeny global fit analyses both found significant evidence of evolutionary congruence between the bat and Bartonella phylogenies (PACo sum of squared residuals = 21.69, P < 0.0001; ParaFit sum of squared residuals = 21.97, P = 0.001). All PACo and ParaFit tests using sampled posterior Bartonella trees showed significant congruence, demonstrating that the results are robust with respect to phylogenetic uncertainty in the Bartonella tree.
Many of the supported links from the PACo and ParaFit tests are between bat species and Bartonella OTUs in the same community (indicated by line colors in Figure 5) and had high network edge weights (indicated by line width in Figure 5). Thirty-one out of the 81 (38%) bat–Bartonella links were highly supported by either PACo or ParaFit (Figure S5 and Table S16). These links were considered highly supported if the upper limit of the PACo jackknife 95% confidence interval was below the mean of all the squared residuals or if the ParaFit F1 statistic was assigned a P-value < 0.01. Of these 31 highly supported links, 22 (71%) were between bat species and Bartonella OTUs identified as being in the same community (Table S16). This proportion of highly supported links in the same community was higher than the proportion of less supported or unsupported links in the same community (χ2 = 3.23, df = 1, P = 0.036). Additional details on individual host–parasite links can be found in the Supplementary Material.
Figure 5. Comparisons of bat species and Bartonella OTU phylogenies. (A) Procrustes superimposition plot with Bartonella OTUs (open circles) and bat species (open triangles). In this plot, the axes represent the principle components of the bat species phylogeny. Bat species are projected onto the two main axes explaining the most variation. Names of bat families and subfamilies are indicated for clusters of species. The Bartonella OTUs are then projected and rotated to fit the bat phylogeny by minimizing the residual distances for each bat-Bartonella association (connected by lines). (B) Cophylogeny plot with the bat species phylogeny on the left side and the Bartonella OTU phylogeny on the right side. Lines in the middle connect bat species to Bartonella OTUs based on sequence data collected from bats or associated ectoparasites. Species abbreviations for bat species are listed in Table 1. In both (A) and (B), bat-Bartonella links colored by the community if both the bat species and Bartonella were placed in the same community (as in Table S14). Line widths in both (A) and (B) are proportional to network edge weight and the line transparency depend on the link support from PACo and ParaFit analyses.
Bayesian Prediction of Bartonella Host Transition Rates
Using the symmetrical rate partition model, only 30 (14%) out of the possible 210 bat host species combinations had significant Bayes factors (BF > 3) from the stochastic search variable selection procedure (Figure 6 and Table S17). The median host-switching rates varied from 0.37 for E. nilssonii and E. serotinus to 1.83 for My. blythii and My. myotis. Large host-switching rates tended to have higher BF support (Pearson's R = 0.52, t = 3.23, df = 28, P = 0.0032), and rates tended to be biased toward bat species in the same family and subfamily (Figure 6). Of the 30 rates, 25 (83%) were between bat species in the same family. This is significantly higher than the expected proportion (126/210, 60%) based on the number of possible species combinations that are in the same family (χ2 = 5.17, df = 1, P = 0.012). The majority of host-switching rates (25/30, 90%) have species pairs in the same identified community (11/30, 37%) or in a community containing bats from the same family (16/30, 53%). Fourteen of the 30 host-switching rates (47%) had species pairs that co-roost during summer months (Table S2).
Figure 6. Graph of Bartonella transition rates among European bat species. (A) Transition rates (median of the posterior distribution) with significant Bayes factors (BF > 3) are drawn as edges connecting nodes representing bat species. (B) Bayes factors are plotted as 2 ln K, where K is the Bayes factor support. Nodes were colored according to bat families/subfamilies. Edge colors were scaled according to median transition rates in (A) or Bayes factor support in (B). Details of transition rates and Bayes factor calculations can be found in Table S17. Species abbreviations for bat species are listed in Table 1.
The top regression model according to AICc (Table S18) for Spearman Bartonella dissimilarity was the model with host phylogenetic distance, geographic range overlap, and roost sharing covariates (AICc = 496). The top model for the binomial Bartonella dissimilarity was the model containing only the host phylogenetic distance covariate (AICc = 562). The top model for the Cao dissimilarity contained the host phylogenetic distance and roost sharing covariates (AICc = 565).
All top models were statistically significant (Table S18) and estimated a statistically significant positive effect of host phylogenetic distance on Bartonella dissimilarity among bat species (Table S19). Since the three candidate models (Spearman, binomial, and Cao) all use different data for regression, we cannot compare them via AICc. Instead, we compared them according to the proportion of variance in Bartonella dissimilarity explained (adjusted R2), the regression mean squared error, and the cross-validation mean squared error. The model using the Spearman correlation explained more variation in Bartonella dissimilarity and had lower mean squared error than the binomial and Cao indices (Table S18). From this model, we can infer that bats are more likely to have dissimilar Bartonella assemblages if they are more distantly related to each other (t = 8.6, df = 206, P < 0.0001), bats have little overlap in their geographic ranges (t = −2.43, df = 206, P = 0.016), and bats do not roost together during the summer (t = −3.81, df = 206, P = 0.00019). Between these three covariates, host phylogenetic distance explains more variation (adjusted partial R2 = 0.26) than geographic range overlap (adjusted partial R2 = 0.023) or roost sharing (adjusted partial R2 = 0.061). This difference in explanatory power is significant because the 95% bootstrap confidence intervals for their relative importance do not overlap (Table S19). Individual Mantel tests confirmed that phylogenetic distance had the strongest correlation with Bartonella dissimilarity (Table S20).
The second candidate model set used host phylogenetic distance, geographic range overlap, and roost sharing to predict ectoparasite dissimilarity. The top model chosen for all dissimilarity indices (Spearman, binomial, and Cao) included all covariates. All three models were statistically significant (Table S18) and estimated the effect of host phylogenetic distance to be positive and the effects of geographic range overlap and roost sharing to be negative (Table S19). According to the Spearman model, bats are more likely to have more dissimilar ectoparasite assemblages if they are more distantly related to one another (t = 3.85, df = 206, P = 0.00016), they have less geographic range overlap (t = 3.04, df = 206, P = 0.0026), and they do not roost together in the summer (t = −3.98, df = 206, P < 0.0001). All three covariates explain similar amounts of variation, and they do not significantly differ in their explanatory power since their relative importance confidence intervals overlap (Table S19). Mantel tests confirmed these findings with correlations of similar magnitude across all the three predictors (Table S20). Therefore, since ectoparasite dissimilarity is explained by host phylogenetic distance, geographic range overlap, and roost sharing, we observe no additional effect of ectoparasite dissimilarity on Bartonella dissimilarity after accounting for these effects.
Finally, the third set of candidate models used the number of gltA sequences from the least sampled bat species plus all four other covariates. The first covariate was used to account for the observed sampling bias in host-switching rates discussed above. We evaluated models using both mean and median host-switching rates taken from the posterior distributions of rates from the Bayesian ancestral state reconstruction analysis. Following model selection by AICc, both the mean and median host-switching models contained covariates for the least sampled bat species and host phylogenetic distance. Therefore, all the top Spearman, binomial, and Cao models were identical. Both the mean and median models were significant (Table S18) and estimated a positive effect for sampling and a negative effect for host phylogenetic distance (Table S19). The median host-switching model explained more variation and had lower mean squared error than the mean model (Table S19). This model indicates that after accounting for a positive sampling bias (t = 4.59, df = 27, P < 0.0001), Bartonella infections are more likely to switch between bat species that are more phylogenetically related (t = −2.32, df = 27, P = 0.028).
For all top model sets (Spearman Bartonella dissimilarity, Spearman ectoparasite dissimilarity, and median host-switching rates), we confirmed the fit of regression models visually by plotting the standardized data against the standardized partial regression coefficients. The plots for Spearman dissimilarity (Figure S6) confirm the positive effect of host phylogenetic distance and the negative effects of geographic range overlap and roost sharing on both Bartonella dissimilarity and ectoparasite dissimilarity. We also confirm the negative effect of host phylogenetic distance and positive sampling effect on Bartonella host-switching rates (Figure S6).
Diversity and Structure of Bat–Ectoparasite–Bartonella Communities
Using a multifaceted analytical approach (Figure 1), we explored the complex nature of Bartonella associations with bats and their ectoparasites in nine European countries. Our first objective was to measure the diversity of Bartonella infections in European bats and their ectoparasites and to analyze the structure of bat–ectoparasite–Bartonella communities. We observed high PD among Bartonella infections in our samples (Table 1), identifying 20 Bartonella OTUs that likely represent distinct species (Figure S4). Network analysis provided support for our hypothesis that associations between bats, ectoparasites, and bacteria could be resolved into identifiable communities. We detected seven distinct bat–ectoparasite–Bartonella communities (Figure 4) that separate by bat phylogeny, geographic ranges, and roosting patterns.
However, the separation of parasites among Myotis spp. bats in our study illustrates the complexity of ecological factors that shape host–parasite associations. Bartonella species found in Myotis spp. bats and their ectoparasites were partitioned into four distinct communities: Min/Myo, VespA, VespB, and VespE (Figure 4 and Table S14). All of these communities were separate from the Rhinolophus-associated Rhi community. This is notable especially for the Min/Myo community because there is overlap in the geographic range and habitat usage (caves) among the bat species in the Rhi and Min/Myo communities. This suggests that distinct host–parasite associations can form despite hosts living in sympatry. Community VespA separates from Min/Myo because the hosts roost mainly in trees during the summer and only use caves during the winter. Species in community VespB use a variety of roosts with little overlap; some are long-distance migrants, and others have relatively northern distribution within Europe. The parasite communities of Myotis spp. bats have thus been shaped by a mixture of different factors that reduce exchange of parasites among hosts, providing motivation for analyzing patterns of Bartonella sharing using a linear regression approach as we have done.
Evolutionary Patterns of Bat–Bartonella Associations
The second objective of our study was to understand the phylogenetic patterns generated from host–parasite associations. We expected that the phylogeny of Bartonella species would exhibit significant clustering by bat taxonomy and would have significant congruence with the phylogeny of bat species. As predicted, our tip-association and cophylogeny tests demonstrated that Bartonella lineages strongly cluster by host taxonomy, and the structure of the Bartonella phylogeny is congruent with the host phylogeny. However, these patterns are not entirely consistent with a pattern of strict cospeciation with bats. Rather, Bartonella lineages infecting bats appear to be polyphyletic, suggesting a more complex history of host-switching and possibly multiple introductions from other mammals deeper in the evolutionary tree (McKee et al., 2017; Urushadze et al., 2017; Frank et al., 2018). Cospeciation of hosts and parasites is a rare phenomenon in general (de Vienne et al., 2013), and several studies of bat viruses have shown that host-switching is the more dominant macroevolutionary force shaping microbial evolution than cospeciation (Cui et al., 2007; Mélade et al., 2016; Anthony et al., 2017). Previous research on Bartonella associations in bats and rodents showed that host-switching is more common than cospeciation (Lei and Olival, 2014). As we will discuss more below, it is more likely that the observed congruence of Bartonella and bat phylogenies is driven by phylogenetic bias in microbial host-switching, such that historical host shifts are more likely to happen between closely related species.
Predictors of Bartonella Sharing and Host-Switching Among Bat Species
Our last objective was to identify the ecological and evolutionary constraints that lead to Bartonella host specificity. We hypothesized that host phylogenetic distance, ectoparasite sharing, geographic range overlap, and roost sharing would be predictors of Bartonella sharing and host-switching rates among bat species. Our regression analysis (Figure S6; Table S18, Table S19) demonstrated that bats that are more phylogenetically related, overlap more in their geographic ranges, and share roosts are more likely to share Bartonella species, with phylogenetic distance being the most important predictor. Ectoparasite sharing between bats had no significant effect on Bartonella sharing after accounting for its own correlation with phylogenetic distance, geographic range overlap, and roost sharing. Finally, our analysis of historical host-switching rates showed that Bartonella lineages are biased to switching between phylogenetically related hosts.
Our regression results explaining variation in Bartonella sharing among bats are in agreement with previous work on bat viruses and other systems. Longdon et al. (2011) demonstrated using Drosophila sigma viruses that the host phylogeny explains most of the variation in viral replication among host species. Primates have more similar parasite communities if they are phylogenetically closely related and inhabit the same region (Davies and Pedersen, 2008). Streicker et al. (2010) found that the frequency of cross-species transmission (CST) of rabies virus between bat species (similar to our measure of Bartonella sharing) increases with decreasing phylogenetic distance and increasing geographic overlap. A later study analyzing these same data confirmed that host phylogenetic distance is a key determinant of rabies CST while other ecological covariates including roost structures, wing aspect ratio, wing loading, and body size were poor predictors (Faria et al., 2013). Luis et al. (2015) found that bat phylogeny and sympatry explained viral sharing in bats, with sympatry being the more important predictor. In addition, viral sharing communities of bats segregated by geographic regions. A recent global analysis of virus sharing in bats showed that after accounting for publication bias, bat species are more likely to share viruses if they have more geographic overlap and they roost in caves (Willoughby et al., 2017). Among cave-roosting bats, species shared viruses more frequently if they overlapped geographically and were documented as sharing roosts. These patterns indicate that host phylogeny and geographic overlap are general predictors of parasite communities in bats. These apparent biases in the arrangement of Bartonella and ectoparasite species among bat species likely contributed to our ability to identify communities of highly interacting bat, ectoparasite, and Bartonella species that tend to cluster by bat family and roosting patterns.
These results demonstrate how key ecological factors constrain the host range of Bartonella species in bats. Yet how do we explain the observed congruence between bat and Bartonella phylogenies? As noted above, we believe that this pattern is best explained by a phylogenetic bias in microbial host-switching, where host shifts occur more frequently between closely related species. This bias, if persistent over the evolution of Bartonella lineages, could produce a Bartonella tree that is largely congruent with the host tree without the need for strict cospeciation (Charleston and Robertson, 2002; de Vienne et al., 2013). This is in line with ecological fitting, which allows host colonization for ecological specialists prior to the evolution of novel capabilities for host exploitation (Araujo et al., 2015). We observed a negative relationship between historical Bartonella host-switching rates and host phylogenetic distance, but no significant relationship with ectoparasite sharing, geographic overlap, or roost sharing. These findings are similar to previous research on bat rabies showing that the host phylogeny is the strongest predictor of rabies host shifts (Streicker et al., 2010; Faria et al., 2013). Thus, biological constraints on parasite shifts among hosts are expected to be the dominant force that shapes host specificity over evolutionary time.
Influence of Ectoparasites on Bartonella Host Specificity
We noted above that ectoparasite sharing failed to explain additional variation in Bartonella sharing between bat species after accounting for the effects of host phylogenetic distance, geographic overlap, and roost sharing. In addition, ectoparasite sharing, geographic overlap, and roost sharing were not included as significant predictors of historical Bartonella host-switching rates. This suggests that the forces of host phylogenetic distance, geographic overlap, and ecological interactions may act on the assembly of ectoparasite and Bartonella communities independently and that host-associated, vector-borne microorganisms come to be vectored by the available ectoparasite communities associated with each host. This is supported by the fact that Bartonella lineages in European bats appear to be associated with a polyphyletic assemblage of arthropods, including ticks, mites, hemipteran bugs, fleas, and flies. Since ectoparasites exhibit differences in life history traits and among-host dispersal mechanisms (Giorgi et al., 2004; Dick and Patterson, 2006; Reckardt and Kerth, 2009), generalist ectoparasites may be more influential in spreading microorganisms among species inhabiting the same environment while specialist ectoparasites are important for the maintenance of microorganisms in separate host species. This broad vector usage could explain why Bartonella infections are so prevalent in bats.
The variation in ectoparasite life history traits and host specificity would also be expected to influence microbial host specificity, though in contrasting ways. While ectoparasites may develop their own specificity for particular host species, as observed in bat wing mites and bat flies (Bruyndonckx et al., 2009; Sándor et al., 2018), these associations between hosts and vectors will predominantly either compound or counteract the isolation already occurring as microbes develop associations with host species. In the first case, the effects of ectoparasite host specificity on microbial evolution may not be statistically separable from the overriding effect of microbial host specificity, as we noted above. On the other hand, generalist vectors could be seen as simply adding noise to the associations of specialist microorganisms through accidental associations with atypical hosts. As long as noise does not totally obscure the predominant host–microbe associations, then it will be possible to measure the host-specific signal, either statistically or through genetic data. For example, Withenshaw et al. (2016) were able to show that distinct genetic variants of Bartonella species separately infect two sympatric rodent species despite the presence of generalist flea vectors. Since the two host species have differences in microhabitat usage and activity patterns, they would rarely have opportunity to exchange fleas. This can lead to covert microbial host specificity even when vectors are host generalists. These patterns suggest that what separates vector-borne microorganisms from directly or environmentally transmitted microorganisms in terms of their host specificity is the added layer of vector host specificity, which will either inflate host specificity already present in the microbe or dilute host-specific patterns through noisy associations. In the case of Bartonella communities in bats, host specificity is clear despite the presence of generalist or polyxenous vectors (e.g., A. vespertilionis, C. pipistrelli, and S. myoti).
Despite the patterns noted above, we should not rule out the possibility that ectoparasites may contribute to the evolution of vector-borne microorganisms like Bartonella. The genus Bartonella appears to have evolved from insect gut symbionts that transitioned to a parasitic lifestyle after adapting to blood-feeding arthropods (Segers et al., 2017). Coevolutionary processes in early Bartonella lineages associated with blood-feeding arthropods may have influenced later patterns in Bartonella associations with mammalian groups. Gene exchange between Bartonella and other arthropod symbionts could also influence the formation of distinct phylogenetic lineages deep in the evolutionary tree (Zhu et al., 2014). The phylogeny of ectoparasite groups may help to understand these processes more fully, but considering the polyphyly of arthropod groups carrying Bartonella observed in our study, it may be more practical to study the influence of ectoparasite host specificity within particular arthropod groups (e.g., bat flies) on Bartonella host specificity and macroevolution. Future studies could use sequences from multiple genetic loci in hosts, vectors, and Bartonella to generate time-calibrated phylogenies and compare the influence of evolutionary processes on structural and temporal patterns within trees across trophic scales.
It is also possible that host and ectoparasite phylogenetic structure in geographically separate populations may influence microevolutionary patterns in associated microbes. Historical processes, such as the postglacial recolonization of regions of Europe by bats (Flanders et al., 2009; Dool et al., 2013), or patterns of host and ectoparasite dispersal across distant locations (Bruyndonckx et al., 2009; Witsenburg et al., 2015; van Schaik et al., 2018) may lead to the formation of distinct microbial lineages. These types of analyses demand additional genetic data to detect fine distinctions between related lineages and are thus beyond the scope of this current work but would be fruitful avenues for future research on vector-borne microorganisms like Bartonella.
Based on our results and previous work, we suggest that while geographic overlap, ecological interactions (e.g., roost sharing), and ectoparasite sharing provide the necessary conditions for Bartonella transmission between hosts, the success of transmission and perhaps an eventual host shift will ultimately depend on biological compatibility between the host and the microbe, which can be predicted by the phylogenetic distance between hosts (Pedersen and Davies, 2009). Strong patterns of host specificity in microbial communities can still be observed even when generalist ectoparasites are present, and while specialist vectors may be present in the community, their effects on microbial host specificity would be expected to be more pronounced for a generalist microorganism and not a specialist, wherein host specificity of vector and microorganism would not be statistically independent.
Bat-Associated Bartonella Species as Zoonoses
Beyond the scientific insights produced by this study, there may be additional practical value in our results. Several Bartonella species are known to be human pathogens, and new cases of zoonotic bartonellosis are consistently being described (Roux et al., 2000; Kosoy et al., 2003, 2010b; Iralu et al., 2006; Chomel and Kasten, 2010; Bai et al., 2012; Kandelaki et al., 2016; Vayssier-Taussat et al., 2016). While we recognize that not all of the 20 putative Bartonella species described in this study may have the potential to infect humans, OTU19 and OTU26 have been previously found to infect humans (Veikkolainen et al., 2014; Urushadze et al., 2017). These particular OTUs are also strongly linked with two ectoparasite species, C. pipistrelli and A. vespertilionis, which are known to sporadically bite humans (Estrada-Peña and Jongejan, 1999; Whyte et al., 2001).
A recent report also detected antibodies to a Bartonella strain specific to Egyptian fruit bats (Rousettus aegyptiacus) in eight people from a community in Nigeria (Bai et al., 2018). Members of the community enter caves and capture bats for consumption or sale as part of an annual festival. Thus, in addition to exposure to bat ectoparasites that may bite humans, these practices may provide alternative routes for human exposure to bat-associated Bartonella species. These include direct exposure through handling of bats (including possible bites or scratches), contamination of open wounds with blood or bat excreta during the capture process, or contamination of wounds from the excreta of bat ectoparasites. Bartonella DNA has been previously detected in bat guano, urine, and saliva (Veikkolainen et al., 2014; Banskar et al., 2016; Dietrich et al., 2017; Becker et al., 2018); however, more studies are needed to confirm that viable bacteria are present in the excreta of bats or their ectoparasites. Nonetheless, there is accumulating evidence that bat-associated Bartonella species may present an infection risk to human populations.
While assessing the zoonotic risks from the environment aids the prevention and control of human diseases, loss and disturbance of roost sites, among other reasons out of sensationalized fear of zoonotic diseases, is one of the threats for endangered bat species (López-Baucells et al., 2018). Disturbance of bat roosts in such sites could even be counterproductive if bat ectoparasites seek alternative food sources, potentially leading to infections of atypical host species such as humans or pets. Therefore, humans should avoid unnecessary contact with bats or their ectoparasites. Bat colonies in attics or walls of buildings do not pose a health risk as long as they are excluded from parts of the house occupied by humans (Tuttle, 2005). For individuals who may encounter bats within caves or other habitats, including tourists, scientists, and bat or guano harvesters, it is advisable to minimize their disturbance of the animals and to limit contact with bats, bat ectoparasites, and their excreta through protective equipment.
The analyses we have performed greatly expand our understanding of the complex ecology and evolution of Bartonella in bats. Our approach, which integrates across levels of parasitism and explores ecological and evolutionary patterns, could be applied to Bartonella in bats outside of Europe, to other Bartonella species associated with different mammalian orders, and possibly to other complex vector-borne diseases. However, we recognize that this study is observational and correlative and has limitations to the data that must be acknowledged.
Regarding our phylogenetic analysis, we needed to confirm that the gltA gene serves as an accurate marker for assessing the evolutionary history of the Bartonella lineages being studied. As detailed in the Supplementary Material, we are satisfied that the lineages identified in the gltA tree have a similar topology to a tree that uses additional markers, that recombination within and among loci has not significantly distorted evolutionary patterns, and that the topology of our tree is similar to trees assembled using another approach, neighbor-joining. Furthermore, by using 100 posterior sampled Bartonella trees in our tip-association and global fit tests, we demonstrated that our results are robust with respect to phylogenetic uncertainty in the parasite tree.
We recognize that our estimates of Bartonella diversity (number of OTUs and Faith's index) in European bats are limited by the extent of current sampling, with a significant positive correlation between these measures and the number of gltA sequences obtained per species. Thus, we probably do not yet have an accurate survey of Bartonella diversity in some of our bat species with few gltA sequences. Additional sequencing of Bartonella strains from some poorly sampled host species (e.g., E. nilssonii, My. capaccinii, My. mystacinus, Pi. pygmaeus) could clarify their membership within communities. Additional sampling of these species would also allow researchers to perform rarefaction analyses to determine if bat species actually differ in the number of Bartonella species that infect them. After accounting for this sampling effect, correlational analyses could attempt to explain this variation in Bartonella diversity using bat traits, as has been done successfully for viral diversity (Luis et al., 2013; Gay et al., 2014; Maganga et al., 2014; Webber et al., 2017; Willoughby et al., 2017).
Moreover, the Bartonella sequence data we compiled for this study were derived from numerous studies that relied on convenience sampling from numerous sources (mist netting, bat boxes, roosts, and dead bats) over varying time periods. No information on host species richness, host density, or community dynamics of host species within roosts was collected. Better structured and controlled sampling strategies could capture such data and would broaden our understanding of parasite diversity and persistence within bat populations.
We also acknowledge that phylogenetic distance only approximates the process of microbial host adaptation. Bartonella traits like the presence of particular secretion system genes and effector proteins (Harms et al., 2017) or bat traits linked to immune function, such as major histocompatibility complex or toll-like receptor alleles (Baker et al., 2013), may be better at predicting whether bats will share Bartonella species. Some bat species may also be more tolerant or resistant to certain Bartonella species, as has been shown for the fungus that causes white nose syndrome (Frank et al., 2014). This could help explain additional variation in Bartonella communities among phylogenetically related and sympatric bat species. Such traits have not been explored in the Bartonella or bat species from our study, but would be useful for understanding how bat species differ in their Bartonella prevalence and diversity.
Finally, we used geographic range overlap as a proxy for spatiotemporal proximity in our regression analysis. This assumes that if species share a roosting preference, the species will interact in a way that could lead to parasite exchange. Some species may interact more frequently within roosts than others, such as by using the same microhabitat or by overlapping within roosts during the same season; thus, there would be a greater potential for ectoparasite and microbial transmission. Fluctuations in colony size and bat community composition between summer maternity colonies and winter roosts would also be expected to influence parasite exchange (Dietrich et al., 2018). Our sampling approach did not capture these important variables and should be considered in future studies.
Our analysis of Bartonella infections in European bats and their ectoparasites has increased our understanding of how Bartonella species are segregated among sympatric hosts and how these ecological associations influence the evolution of Bartonella lineages. We find that while the host phylogeny primarily explains how bats share Bartonella and how Bartonella lineages evolve, the complete evolutionary history of Bartonella in bats may involve additional processes, including multiple introductions from other mammalian orders and biogeographic separation of hosts. Additional sampling and phylogenetic analysis of Bartonella from bats and other mammals and their ectoparasites could help to shed more light on these processes. As we have demonstrated, Bartonella is a productive system for studying complex host–parasite associations. The methods that we have used in this study and our findings regarding the important processes that constrain parasite evolution may be applicable to other vector-borne microorganisms, such as those carried by ticks. For example, studying Anaplasma spp. infections in different carnivore or ungulate species taking into account their phylogeny, geographic overlap, and tick associations could provide information about spillover potential to humans and domestic animals (Dugat et al., 2015; André, 2018). Investigating multitrophic interactions would also help to understand the processes behind the segregation of Borrelia genospecies among vertebrate hosts (Estrada-Peña et al., 2016). Comparison among these various systems may reveal general patterns in the ecology and evolution of host–vector–parasite associations.
Representative sequences from all ectoparasites and Bartonella OTUs have been submitted to GenBank. Accession numbers for Bartonella gltA sequences are MK140184–MK140336 and accession numbers for ectoparasites are MK140016–MK140183. Phylogenetic trees, R code, and additional data sheets are available on GitHub (https://github.com/clifmckee/eurobats).
No live bat was harmed for this study. Authorization for bat capture was provided by the Netherlands Enterprise Agency over the entire sampling period (ontheffing Flora en Faunawet permits FF75a/2008/033, FF75a2012/37a, and FF75a/2016/048 assigned to the Dutch Mammal Society) with permissions of all site owners; the National Inspectorate for Environment, Nature and Water (Hungary); the Underground Heritage Commission (Romania); the Flemish Forest and Nature Agency (ANB/BL-FF/V15-00095); and the Wallonian Department of Nature and Forests (SPW-DNF-DGo3 permit réf déro 2013/RS/n°15 and dero 2014/RS/n°12). Bat banding license numbers are 59/2003, 305/2015, 46/2016, and TMF-493/3/2005. Bats were handled according to the current laws of animal welfare regulation (L206/2004), and the Research Bioethics Commission of USAMV CN approved all working protocols for data collection in Romania.
HS and AK conceived the study. AK, AS, and TG performed the sample collection in the field. AS and GF organized parts of the sample collection and contributed to the study design. DD, A-JH, MF, GF, and TG contributed important samples to the study. MF identified bat flies. AK performed laboratory and molecular analyses. CM performed the phylogenetic and statistical analyses. CM and AK wrote the first draft of the manuscript. All authors read and approved the final version of the manuscript.
Part of this study was financially supported by the Dutch Ministry of Health, Welfare and Sport. This research was supported by grant PN-II-RU-TE-2014-4-1389, the grant In the light of evolution: theories and solutions (GINOP-2.3.2-15-2016-00057), and the János Bolyai Research Scholarship of the Hungarian Academy of Science (AS and GF). No specific funds, beyond the annual budget, were acquired from the US Centers for Disease Control and Prevention Division of Vector-Borne Diseases to support this study. The findings and conclusions in this report are those of the authors and do not necessarily represent the official position of CDC.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We thank René Janssen, Wout Willems, Thijs Molenaar, Jan Boshamer, Pierrette Nyssen, and members of the Dutch Bat Working Group (VLEN) for their assistance with field sampling. We thank Pepijn Kamminga for sharing samples from Naturalis Biodiversity Center. We thank the members of the Webb and Kosoy laboratories for their suggestions and support during this study. We also thank Michael Ryan, Sarah Gallup, and DeAna Nasseth for their comments on early drafts of the manuscript.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2019.00069/full#supplementary-material
Agnarsson, I., Zambrana-Torrelio, C. M., Flores-Saldana, N. P., and May-Collado, L. J. (2011). A time-calibrated species-level phylogeny of bats (Chiroptera, Mammalia). PLoS Curr. 3:RRN1212. doi: 10.1371/currents.RRN1212
Anderson, M. J., and Millar, R. B. (2004). Spatial variation and effects of habitat on temperate reef fish assemblages in northeastern New Zealand. J. Exp. Mar. Bio. Ecol. 305, 191–221. doi: 10.1016/j.jembe.2003.12.011
André, M. R. (2018). Diversity of Anaplasma. and Ehrlichia./Neoehrlichia. agents in terrestrial wild carnivores worldwide: implications for human and domestic animal health and wildlife conservation. Front. Vet. Sci. 5:293. doi: 10.3389/fvets.2018.00293
Anh, P. H., Van Cuong, N., Son, N. T., Tue, N. T., Kosoy, M. Y., Woolhouse, M. E. J., et al. (2015). Diversity of Bartonella. spp. in bats, southern Vietnam. Emerg. Infect. Dis. 21, 1266–1267. doi: 10.3201/eid2107.141760
Araujo, S. B. L., Braga, M. P., Brooks, D. R., Agosta, S. J., Hoberg, E. P., von Hartenthal, F. W., et al. (2015). Understanding host-switching by ecological fitting. PLoS ONE 10:e0139225. doi: 10.1371/journal.pone.0139225
Bai, Y., Kosoy, M. Y., Boonmar, S., Sawatwong, P., Sangmaneedet, S., and Peruski, L. F. (2010). Enrichment culture and molecular identification of diverse Bartonella. species in stray dogs. Vet. Microbiol. 146, 314–319. doi: 10.1016/j.vetmic.2010.05.017
Bai, Y., Kosoy, M. Y., Diaz, M. H., Winchell, J., Baggett, H., Maloney, S. A., et al. (2012). Bartonella vinsonii subsp. arupensis. in humans, Thailand. Emerg. Infect. Dis. 18, 989–991. doi: 10.3201/eid1806.111750
Bai, Y., Osinubi, M. O. V., Osikowicz, L., McKee, C., Vora, N. M., Rizzo, M. R., et al. (2018). Human exposure to novel Bartonella. species from contact with fruit bats. Emerg. Infect. Dis. 24, 2317–2323. doi: 10.3201/eid2412.181204
Balbuena, J. A., Poisot, T., Hutchinson, M., and Cagua, F. (2016). PACo: Procrustes Application to Cophylogenetic Analysis. Available online at: https://cran.r-project.org/package=paco
Banskar, S., Bhute, S. S., Suryavanshi, M. V., Punekar, S., and Shouche, Y. S. (2016). Microbiome analysis reveals the abundance of bacterial pathogens in Rousettus leschenaultii. guano. Sci. Rep. 6:36948. doi: 10.1038/srep36948
Barton, K. (2016). MuMIn: Multi-Model Inference. Available online at: http://cran.r-project.org/package=MuMIn
Becker, D. J., Bergner, L. M., Bentz, A. B., Orton, R. J., Altizer, S., and Streicker, D. G. (2018). Genetic diversity, infection prevalence, and possible transmission routes of Bartonella spp. in vampire bats. PLoS Negl. Trop. Dis. 12:e0006786. doi: 10.1371/journal.pntd.0006786
Becker, R. A., Wilks, A. R., Brownrigg, R., Minka, T. P., and Deckmyn, A. (2016). Maps: Draw Geographical Maps. Available online at: https://cran.r-project.org/package=maps
Bivand, R., Lewin-Koh, N., Pebesma, E., Archer, E., Baddeley, A., Bearman, N., et al. (2017). Maptools: Tools for Reading and Handling Spatial Objects. Available online at: https://cran.r-project.org/package=maptools
Blondel, V. D., Guillaume, J.-L., Lambiotte, R., and Lefebvre, E. (2008). Fast unfolding of communities in large networks. J. Stat. Mech. Theory Exp. 2008:P10008. doi: 10.1088/1742-5468/2008/10/P10008
Braks, M., van der Giessen, J., Kretzschmar, M., van Pelt, W., Scholte, E.-J., Reusken, C., et al. (2011). Towards an integrated approach in surveillance of vector-borne diseases in Europe. Parasit. Vect. 4:192. doi: 10.1186/1756-3305-4-192
Bruyndonckx, N., Dubey, S., Ruedi, M., and Christe, P. (2009). Molecular cophylogenetic relationships between European bats and their ectoparasitic mites (Acari, Spinturnicidae.). Mol. Phylogenet. Evol. 51, 227–237. doi: 10.1016/j.ympev.2009.02.005
Charleston, M. A., and Robertson, D. L. (2002). Preferential host switching by primate lentiviruses can account for phylogenetic similarity with the primate phylogeny. Syst. Biol. 51, 528–535. doi: 10.1080/10635150290069940
Concannon, R., Wynn-Owen, K., Simpson, V. R., and Birtles, R. J. (2005). Molecular characterization of haemoparasites infecting bats (Microchiroptera.) in Cornwall, UK. Parasitology 131, 489–496. doi: 10.1017/S0031182005008097
Corduneanu, A., Sándor, A. D., Ionică, A. M., Hornok, S., Leitner, N., Bagó, Z., et al. (2018). Bartonella. DNA in heart tissues of bats in central and eastern Europe and a review of phylogenetic relations of bat-associated bartonellae. Parasit. Vectors 11:489. doi: 10.1186/s13071-018-3070-7
Cui, J., Han, N., Streicker, D. G., Li, G., Tang, X., Shi, Z., et al. (2007). Evolutionary relationships between bat coronaviruses and their hosts. Emerg. Infect. Dis. 13, 1526–1532. doi: 10.3201/eid1310.070448
Davies, T. J., and Pedersen, A. B. (2008). Phylogeny and geography predict pathogen community similarity in wild primates and humans. Proc. R. Soc. B Biol. Sci. 275, 1695–1701. doi: 10.1098/rspb.2008.0284
de Vienne, D., Refrégier, G., López-Villavicencio, M., Tellier, A., Hood, M., and Giraud, T. (2013). Cospeciation vs. host-shift speciation: methods for testing, evidence from natural associations and relation to coevolution. N. Phytol. 198, 347–385. doi: 10.1111/nph.12150
Dick, C. W., and Patterson, B. D. (2006). “Bat flies: obligate ectoparasites of bats,” in Micromammals and Macroparasites: From Evolutionary Ecology to Management. eds S. Morand, B. R. Krasnov, and R. Poulin (New York, NY: Springer), 179–194.
Dietrich, M., Kearney, T., Seamark, E. C. J., and Markotter, W. (2017). The excreted microbiota of bats: evidence of niche specialisation based on multiple body habitats. FEMS Microbiol. Lett. 364:fnw284. doi: 10.1093/femsle/fnw284
Dietrich, M., Kearney, T., Seamark, E. C. J., Paweska, J. T., and Markotter, W. (2018). Synchronized shift of oral, faecal and urinary microbiotas in bats and natural infection dynamics during seasonal reproduction. R. Soc. Open Sci. 5:180041. doi: 10.1098/rsos.180041
Dool, S. E., Puechmaille, S. J., Dietz, C., Juste, J., Ibáñez, C., Hulva, P., et al. (2013). Phylogeography and postglacial recolonization of Europe by Rhinolophus hipposideros: evidence from multiple genetic markers. Mol. Ecol. 22, 4055–4070. doi: 10.1111/mec.12373
Dugat, T., Lagrée, A.-C., Maillard, R., Boulouis, H.-J., and Haddad, N. (2015). Opening the black box of Anaplasma phagocytophilum. diversity: current situation and future perspectives. Front. Cell. Infect. Microbiol. 5, 1–18. doi: 10.3389/fcimb.2015.00061
Encarnação, J. A., Baulechner, D., and Becker, N. I. (2012). Seasonal variations of wing mite infestations in male Daubenton's bats (Myotis daubentonii.) in comparison to female and juvenile bats. Acta Chiropterol. 14, 153–159. doi: 10.3161/150811012X654367
Estrada-Peña, A., and Jongejan, F. (1999). Ticks feeding on humans: a review of records on human-biting Ixodoidea with special reference to pathogen transmission. Exp. Appl. Acarol. 23, 685–715. doi: 10.1023/A:1006241108739
Estrada-Peña, A., Sprong, H., Cabezas-Cruz, A., de la Fuente, J., Ramo, A., and Coipan, E. C. (2016). Nested coevolutionary networks shape the ecological relationships of ticks, hosts, and the Lyme disease bacteria of the Borrelia burgdorferi. (s.l). complex. Parasit. Vectors. 9:517. doi: 10.1186/s13071-016-1803-z
Faria, N. R., Suchard, M., Rambaut, A., Streicker, D. G., and Lemey, P. (2013). Simultaneously reconstructing viral cross-species transmission history and identifying the underlying constraints. Philos. Trans. R. Soc. B 368:20120196. doi: 10.1098/rstb.2012.0196
Flanders, J., Jones, G., Benda, P., Dietz, C., Zhang, S., Li, G., et al. (2009). Phylogeography of the greater horseshoe bat, Rhinolophus ferrumequinum: contrasting results from mitochondrial and microsatellite data. Mol. Ecol. 18, 306–318. doi: 10.1111/j.1365-294X.2008.04021.x
Folmer, O., Black, M., Hoeh, W., Lutz, R., and Vrijenhoek, R. (1994). DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol. Mar. Biol. Biotechnol. 3, 294–299.
Frank, C. L., Michalski, A., McDonough, A. A., Rahimian, M., Rudd, R. J., and Herzog, C. (2014). The resistance of a North American bat species (Eptesicus fuscus.) to white-nose syndrome (WNS). PLoS ONE 9:e113958. doi: 10.1371/journal.pone.0113958
Frank, H. K., Boyd, S. D., and Hadly, E. A. (2018). Global fingerprint of humans on the distribution of Bartonella. bacteria in mammals. PLoS Negl. Trop. Dis. 12:e0006865. doi: 10.1371/journal.pntd.0006865
Gay, N., Olival, K. J., Bumrungsri, S., Siriaroonrat, B., Bourgarel, M., and Morand, S. (2014). Parasite and viral species richness of Southeast Asian bats: fragmentation of area distribution matters. Int. J. Parasitol. Parasites Wildl. 3, 161–170. doi: 10.1016/j.ijppaw.2014.06.003
Giorgi, M. S., Arlettaz, R., Guillaume, F., Nusslé, S., Ossola, C., Vogel, P., et al. (2004). Causal mechanisms underlying host specificity in bat ectoparasites. Oecologia 138, 648–654. doi: 10.1007/s00442-003-1475-1
Groemping, U., and Matthias, L. (2013). Relaimpo: Relative Importance of Regressors in Linear Models. Available online at: https://cran.r-project.org/package=relimpo
Han, H.-J., Wen, H., Zhao, L., Liu, J.-W., Luo, L.-M., Zhou, C.-M., et al. (2017). Novel Bartonella. species in insectivorous bats, northern China. PLoS ONE. 12:e0167915. doi: 10.1371/journal.pone.0167915
Hanson, B. A., Memisevic, V., and Chung, J. (2016). HiveR: 2D and 3D Hive Plots for R. Available online at: https://cran.r-project.org/package=HiveR
Harms, A., Segers, F. H. I. D., Quebatte, M., Mistl, C., Manfredi, P., Körner, J., et al. (2017). Evolutionary dynamics of pathoadaptation revealed by three independent acquisitions of the VirB/D4 type IV secretion system in Bartonella. Genome Biol. Evol. 9, 761–776. doi: 10.1093/gbe/evx042
Hornok, S., Szoke, K., Görföl, T., Földvári, G., Tu, V. T., Takács, N., et al. (2017). Molecular investigations of the bat tick Argas vespertilionis. (Ixodida: Argasidae) and Babesia vesperuginis (Apicomplexa: Piroplasmida.) reflect “bat connection” between Central Europe and Central Asia. Exp. Appl. Acarol. 72, 69–77. doi: 10.1007/s10493-017-0140-z
Hornok, S., Szoke, K., Kováts, D., Estók, P., Görföl, T., Boldogh, S. A., et al. (2016). DNA of piroplasms of ruminants and dogs in ixodid bat ticks. PLoS ONE. 11:e0167735. doi: 10.1371/journal.pone.0167735
Hornok, S., Szoke, K., Meli, M. L., Sándor, A. D., Görföl, T., Estók, P., et al. (2019). Molecular detection of vector-borne bacteria in bat ticks (Acari: Ixodidae, Argasidae.) from eight countries of the Old and New Worlds. Parasit. Vectors 12:50. doi: 10.1186/s13071-019-3303-4
Iralu, J., Bai, Y., Crook, L., Tempest, B., Simpson, G., McKenzie, T., et al. (2006). Rodent-associated Bartonella. febrile illness, southwestern United States. Emerg. Infect. Dis. 12, 1081–1086. doi: 10.3201/eid1207.040397
IUCN (2014). The IUCN Red List of Threatened Species. Available online at: http://www.iucnredlist.org/
Jaenson, T. G. T., Tälleklint, L., Lundqvist, L., Olsen, B., Chirico, J., and Mejlon, H. (1994). Geographical distribution, host associations, and vector roles of ticks (Acari: Ixodidae, Argasidae.) in Sweden. J. Med. Entomol. 31, 240–256. doi: 10.1093/jmedent/31.2.240
Kandelaki, G., Malania, L., Bai, Y., Chakvetadze, N., Katsitadze, G., Imnadze, P., et al. (2016). Human lymphadenopathy caused by ratborne Bartonella., Tbilisi, Georgia. Emerg. Infect. Dis. 22, 544–546. doi: 10.3203/eid2203.151823
Kembel, S. W., Ackerly, D. D., Blomberg, S. P., Cornwell, W. K., Cowan, P. D., Helmus, M. R., et al. (2014). Picante: R Tools for Integrating Phylogenies and Ecology. Available online at: https://cran.r-project.org/package=picante
Kosoy, M., McKee, C., Albayrak, L., and Fofanov, Y. (2018). Genotyping of Bartonella. bacteria and their animal hosts: current status and perspectives. Parasitology 145, 543–562. doi: 10.1017/S0031182017001263
Kosoy, M. Y., Bai, Y., Sheff, K., Morway, C., Baggett, H., Maloney, S. A., et al. (2010b). Identification of Bartonella. infections in febrile human patients from Thailand and their potential animal reservoirs. Am. J. Trop. Med. Hyg. 82, 1140–1145. doi: 10.4269/ajtmh.2010.09-0778
Kosoy, M. Y., Hayman, D. T. S., and Chan, K.-S. (2012). Bartonella. bacteria in nature: where does population variability end and a species start? Infect. Genet. Evol. 12, 894–904. doi: 10.1016/j.meegid.2012.03.005
Kosoy, M. Y., Murray, M., Gilmore, R. D. Jr., Bai, Y., and Gage, K. L. (2003). Bartonella. strains from ground squirrels are identical to Bartonella washoensis. isolated from a human patient. J. Clin. Microbiol. 41, 645–650. doi: 10.1128/JCM.41.2.645-650.2003
La Scola, B., Zeaiter, Z., Khamis, A., and Raoult, D. (2003). Gene-sequence-based criteria for species definition in bacteriology: the Bartonella. paradigm. Trends Microbiol. 11, 318–321. doi: 10.1016/S0966-842X(03)00143-4
Lei, B. R., and Olival, K. J. (2014). Contrasting patterns in mammal-bacteria coevolution: Bartonella. and Leptospira. in bats and rodents. PLoS Negl. Trop. Dis. 8:e2738. doi: 10.1371/journal.pntd.0002738
Lilley, T. M., Veikkolainen, V., and Pulliainen, A. T. (2015). Molecular detection of Candidatus. Bartonella hemsundetiensis in bats. Vector-Borne Zoonotic Dis. 15, 706–708. doi: 10.1089/vbz.2015.1783
Lilley, T. M., Wilson, C. A., Bernard, R. F., Willcox, E. V., Vesterinen, E. J., Webber, Q. M. R., et al. (2017). Molecular detection of Candidatus. Bartonella mayotimonensis in North American bats. Vector-Borne Zoonotic Dis. 17, 243–246. doi: 10.1089/vbz.2016.2080
Lin, E. Y., Tsigrelis, C., Baddour, L. M., Lepidi, H., Rolain, J.-M., Patel, R., et al. (2010). Candidatus. Bartonella mayotimonensis and endocarditis. Emerg. Infect. Dis. 16, 500–503. doi: 10.3201/eid1603.081673
Lin, J.-W., Hsu, Y.-M., Chomel, B. B., Lin, L.-K., Pei, J.-C., Wu, S.-H., et al. (2012). Identification of novel Bartonella. spp. in bats and evidence of Asian gray shrew as a new potential reservoir of Bartonella. Vet. Microbiol. 156, 119–126. doi: 10.1016/j.vetmic.2011.09.031
Longdon, B., Hadfield, J. D., Webster, C. L., Obbard, D. J., and Jiggins, F. M. (2011). Host phylogeny determines viral persistence and replication in novel hosts. PLoS Pathog. 7:e1002260. doi: 10.1371/journal.ppat.1002260
López-Baucells, A., Rocha, R., and Fernández-Llamazares, Á. (2018). When bats go viral: negative framings in virological research imperil bat conservation. Mamm. Rev. 48, 62–66. doi: 10.1111/mam.12110
Luis, A. D., Hayman, D. T. S., O'Shea, T. J., Cryan, P. M., Gilbert, A. T., Pulliam, J. R. C., et al. (2013). A comparison of bats and rodents as reservoirs of zoonotic viruses: are bats special? Proc. R. Soc. B 280:e20122753. doi: 10.1098/rspb.2012.2753
Luis, A. D., O'Shea, T. J., Hayman, D. T. S., Wood, J. L. N., Cunningham, A. A., Gilbert, A. T., et al. (2015). Network analysis of host–virus communities in bats and rodents reveals determinants of cross-species transmission. Ecol. Lett. 18, 1153–1162. doi: 10.1111/ele.12491
Maganga, G. D., Bourgarel, M., Vallo, P., Dallo, T. D., Ngoagouni, C., Drexler, J. F., et al. (2014). Bat distribution size or shape as determinant of viral richness in African bats. PLoS ONE 9:e100172. doi: 10.1371/journal.pone.0100172
Maindonald, J. H., and Braun, W. J. (2015). DAAG: Data Analysis and Graphics Data and Functions. Available online at: https://cran.r-project.org/package=DAAG
McKee, C. D., Hayman, D. T. S., Kosoy, M. Y., and Webb, C. T. (2016). Phylogenetic and geographic patterns of bartonella host shifts among bat species. Infect. Genet. Evol. 44, 382–394. doi: 10.1016/j.meegid.2016.07.033
McKee, C. D., Kosoy, M. Y., Bai, Y., Osikowicz, L. M., Franka, R., Gilbert, A. T., et al. (2017). Diversity and phylogenetic relationships among Bartonella. strains from Thai bats. PLoS ONE 12:e0181696. doi: 10.1371/journal.pone.0181696
Mélade, J., Wieseke, N., Ramasindrazana, B., Flores, O., Lagadec, E., Gomard, Y., et al. (2016). An eco-epidemiological study of Morbilli-related paramyxovirus infection in Madagascar bats reveals host-switching as the dominant macro-evolutionary mechanism. Sci. Rep. 6:23752. doi: 10.1038/srep23752
Morand, S. (2015). (macro-) Evolutionary ecology of parasite diversity: from determinants of parasite species richness to host diversification. Int. J. Parasitol. Parasites Wildl. 4, 80–87. doi: 10.1016/j.ijppaw.2015.01.001
Morse, S. F., Olival, K. J., Kosoy, M. Y., Billeter, S. A., Patterson, B. D., Dick, C. W., et al. (2012). Global distribution and genetic diversity of Bartonella. in bat flies (Hippoboscoidea, Streblidae, Nycteribiidae.). Infect. Genet. Evol. 12, 1717–1723. doi: 10.1016/j.meegid.2012.06.009
Norman, A. F., Regnery, R., Jameson, P., Greene, C., and Krause, D. C. (1995). Differentiation of Bartonella.-like isolates at the species level by PCR-restriction fragment length polymorphism in the citrate synthase gene. J. Clin. Microbiol. 33, 1797–1803.
Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., Minchin, P. R., et al. (2015). Vegan: Community Ecology Package. Available online at: http://cran.r-project.org/package=vegan
Olival, K. J., Hosseini, P. R., Zambrana-Torrelio, C., Ross, N., Bogich, T. L., and Daszak, P. (2017). Host and viral traits predict zoonotic spillover from mammals. Nature 546, 646–650. doi: 10.1038/nature22975
Paradis, E., Blomberg, S., Bolker, B., Claude, J., Cuong, H. S., Desper, R., et al. (2016). APE: Analyses of Phylogenetics and Evolution. Available online at: https://cran.r-project.org/package=ape
Parker, J., Rambaut, A., and Pybus, O. G. (2008). Correlating viral phenotypes with phylogeny: accounting for phylogenetic uncertainty. Infect. Genet. Evol. 8, 239–246. doi: 10.1016/j.meegid.2007.08.001
Podsiadly, E., Chmielewski, T., Karbowiak, G., Kedra, E., and Tylewska-Wierzbanowska, S. (2010). The occurrence of spotted fever rickettsioses and other tick-borne infections in forest workers in Poland. Vector-Borne Zoonotic Dis. 11, 985–989. doi: 10.1089/vbz.2010.0080
R Core Team (2018). R: A Language and Environment for Statistical Computing. Vienna: R Found. Stat. Comput. Available online at: http://www.r-project.org
Reckardt, K., and Kerth, G. (2009). Does the mode of transmission between hosts affect the host choice strategies of parasites? Implications from a field study on bat fly and wing mite infestation of Bechstein's bats. Oikos. 118, 183–190. doi: 10.1111/j.1600-0706.2008.16950.x
Roux, V., Eykyn, S. J., Wyllie, S., Raoult, D., and Roux, R. (2000). Bartonella vinsonii subsp. berkhoffii. as an agent of afebrile blood culture-negative endocarditis in a human. J. Clin. Microbiol. 38, 1698–1700.
Sanchez, G. (2013). Arcdiagram. Available online at: https://github.com/gastonstat/arcdiagram
Sándor, A. D., Földvári, M., Krawczyk, A. I., Sprong, H., Corduneanu, A., Barti, L., et al. (2018). Eco-epidemiology of novel Bartonella. genotypes from parasitic flies of insectivorous bats. Microb. Ecol. 76, 1076–1088. doi: 10.1007/s00248-018-1195-z
Segers, F. H. I. D., Kešnerová, L., Kosoy, M., and Engel, P. (2017). Genomic changes associated with the evolutionary transition of an insect gut symbiont into a blood-borne pathogen. ISME J. 11:1232–1244. doi: 10.1038/ismej.2016.201
Streicker, D. G., Turmelle, A. S., Vonhof, M. J., Kuzmin, I. V., McCracken, G. F., and Rupprecht, C. E. (2010). Host phylogeny constrains cross-species emergence and establishment of rabies virus in bats. Science 329, 676–679. doi: 10.1126/science.1188836
Stuckey, M. J., Boulouis, H.-J., Cliquet, F., Picard-Meyer, E., Servat, A., Aréchiga-Ceballos, N., et al. (2017a). Potentially zoonotic Bartonella. in bats from France and Spain. Emerg. Infect. Dis. 23, 539–541. doi: 10.3201/eid2303.160934
Stuckey, M. J., Chomel, B. B., de Fleurieu, E. C., Aguilar-Setién, A., Boulouis, H.-J., and Chang, C. (2017b). Bartonella, bats and bugs: a review. Comp. Immunol. Microbiol. Infect. Dis. 55, 20–29. doi: 10.1016/j.cimid.2017.09.001
Szubert-Kruszynska, A., Stanczak, J., Cieniuch, S., Podsiadły, E., Postawa, T., and Michalik, J. (2018). Bartonella. and Rickettsia. infections in haematophagous Spinturnix myoti. mites (Acari: Mesostigmata.) and their bat host, Myotis myotis. (Yangochiroptera: Vespertilionidae.), from Poland. Microb. Ecol. doi: 10.1007/s00248-018-1246-5. [Epub ahead of print].
Theodor, O. (1967). An Illustrated Catalogue of the Rothschild Collection of Nycteribiidae (Diptera) in the British Museum (Natural History). Publication No. 655, London: Trustees of the British Museum (Natural History).
Tijsse-Klasen, E., Fonville, M., Gassner, F., Nijhof, A. M., Hovius, E. K. E., Jongejan, F., et al. (2011). Absence of zoonotic Bartonella. species in questing ticks: first detection of Bartonella clarridgeiae and Rickettsia felis. in cat fleas in the Netherlands. Parasit. Vectors 4:61. doi: 10.1186/1756-3305-4-61
Urushadze, L., Bai, Y., Osikowicz, L. M., McKee, C. D., Kandaurov, A., Kuzmin, I. V., et al. (2017). Prevalence, diversity, and host associations of Bartonella. strains in bats from Georgia (Caucasus). PLoS Negl. Trop. Dis. 11:e0005428. doi: 10.1371/journal.pntd.0005428
van Schaik, J., Dekeukeleire, D., Gazaryan, S., Natradze, I., and Kerth, G. (2018). Comparative phylogeography of a vulnerable bat and its ectoparasite reveals dispersal of a non-mobile parasite among distinct evolutionarily significant units of the host. Conserv. Genet. 19, 481–494. doi: 10.1007/s10592-017-1024-9
van Schaik, J., Dekeukeleire, D., and Kerth, G. (2015). Host and parasite life history interplay to yield divergent population genetic structures in two ectoparasites living on the same bat species. Mol. Ecol. 24, 2324–2335. doi: 10.1111/mec.13171
van Schaik, J., and Kerth, G. (2017). Host social organization and mating system shape parasite transmission opportunities in three European bat species. Parasitol. Res. 116, 589–599. doi: 10.1007/s00436-016-5323-8
Vayssier-Taussat, M., Moutailler, S., Féménia, F., Raymond, P., Croce, O., La Scola, B., et al. (2016). Identification of novel zoonotic activity of Bartonella. spp., France. Emerg. Infect. Dis. 22, 457–462. doi: 10.3201/eid2203.150269
Veikkolainen, V., Vesterinen, E. J., Lilley, T. M., and Pulliainen, A. T. (2014). Bats as reservoir hosts of human bacterial pathogen, Bartonella mayotimonensis. Emerg. Infect. Dis. 20, 960–967. doi: 10.3201/eid2006.130956
Webber, Q. M. R., Fletcher, Q. E., and Willis, C. K. R. (2017). Viral richness is positively related to group size, but not mating system, in bats. Ecohealth. 14, 652–661. doi: 10.1007/s10393-017-1276-3
Wielinga, P. R., Gaasenbeek, C., Fonville, M., de Boer, A., de Vries, A., Dimmers, W., et al. (2006). Longitudinal analysis of tick densities and Borrelia., Anaplasma., and Ehrlichia. infections of Ixodes ricinus. ticks in different habitat areas in The Netherlands. Appl. Environ. Microbiol. 72, 7594–7601. doi: 10.1128/AEM.01851-06
Wilkinson, D. A., Duron, O., Cordonin, C., Gomard, Y., Ramasindrazana, B., Mavingui, P., et al. (2016). The bacteriome of bat flies (Nycteribiidae) from the Malagasy region: a community shaped by host ecology, bacterial transmission mode, and host-vector specificity. Appl. Environ. Microbiol. 82, 1778–1788. doi: 10.1128/AEM.03505-15
Wilson, D., and Reeder, D. M. (2015). Mammal Species of the World, 3rd Edn. Available online at: http://vertebrates.si.edu/msw/mswcfapp/msw/
Withenshaw, S. M., Devevey, G., Pedersen, A. B., and Fenton, A. (2016). Multihost Bartonella. parasites display covert host specificity even when transmitted by generalist vectors. J. Anim. Ecol. 85, 1442–1452. doi: 10.1111/1365-2656.12568
Witsenburg, F., Clément, L., López-Baucells, A., Palmeirim, J., Pavlinić, I., Scaravelli, D., et al. (2015). How a haemosporidian parasite of bats gets around: the genetic structure of a parasite, vector and host compared. Mol. Ecol 24, 926–940. doi: 10.1111/mec.13071
Keywords: Bartonella, ectoparasites, disease ecology, parasite communities, host-switching, network analysis
Citation: McKee CD, Krawczyk AI, Sándor AD, Görföl T, Földvári M, Földvári G, Dekeukeleire D, Haarsma A-J, Kosoy MY, Webb CT and Sprong H (2019) Host Phylogeny, Geographic Overlap, and Roost Sharing Shape Parasite Communities in European Bats. Front. Ecol. Evol. 7:69. doi: 10.3389/fevo.2019.00069
Received: 26 November 2018; Accepted: 22 February 2019;
Published: 26 March 2019.
Edited by:Kevin R. Theis, Wayne State University, United States
Reviewed by:Beza Ramasindrazana, Institut Pasteur de Madagascar, Madagascar
Kendra Phelps, EcoHealth Alliance, United States
Copyright © 2019 McKee, Krawczyk, Sándor, Görföl, Földvári, Földvári, Dekeukeleire, Haarsma, Kosoy, Webb and Sprong. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
†These authors have contributed equally to this work