A Fly on the Cave Wall: Parasite Genetics Reveal Fine-scale Dispersal Patterns of Bats

Dispersal and population structure are key components in quantifying demographic changes of a population, estimating the evolutionary history of species, and informing conservation and management strategies (Kokko and López-Sepulcre, 2006; Ronce, 2007; Clobert et al., 2009). Population connectivity can be difficult to quantify accurately in organisms where methods of direct monitoring (e.g., radio telemetry, GPS tagging) cannot be used and in situations where the genetic signal of dispersal is ambiguous or absent. The genetic signal of dispersal may be low in polyploid organisms; populations that have recently undergone rapid demographic change, have recently diverged, or have experienced hybridization; or when deep sampling of individuals or loci is not possible (Evanno et al., 2005; Latch et al., 2006; Moussy et al., 2013; Dufresne et al., 2014; Milligan et al., 2018). In addition, population genetics cannot be used to detect dispersal when there are pre- or post-zygotic barriers to gene flow between local and dispersed individuals. Parasites offer an independent source of information about their hosts that can help to overcome the difficulties in estimating dispersal, including instances where host dispersal is occurring, but gene flow is not (Whiteman and Parker, 2005; Galbreath and Hoberg, 2012). Previous research on host–parasite coevolution suggests that host-specific parasites, which have remained tightly linked with their hosts through evolutionary time, might provide valuable insights about their host (Nadler et al., 1990; Dybdahl and Lively, 1996; Whiteman and Parker, 2005; Galbreath and Hoberg, 2012; Block et al., 2015; Sweet et al., 2018). However, parasites vary in their specificity and life-history strategies, which may allow incongruities in parameters of population genetics between hosts and their parasites (Barrett et al., 2008; Balvín et al., 2013; Olival et al., 2013; van Schaik et al., 2015; Mazé-Guilmo et al., 2016; Talbot et al., 2017). In addition to extremely specific parasites, we propose that parasites that frequently move between individuals of the same host species (i.e., micropredators), and parasites that do not share an evolutionary history with their host but occur on a single host species within a geographic region also allow us to test hypotheses of host dispersal that we cannot examine using host genetics alone.

Bats (Order Chiroptera) are the only flying mammals, and their mobility makes them interesting models for examining colonization, speciation, and extinction in mammals. It can be difficult to directly measure dispersal in bats due to their size relative to the size of positional data loggers and their high mobility (Moussy et al., 2013). Recent genetic work in the Caribbean has indicated that bat dispersal between some islands is limited due to the distance between islands or secondary contact between cryptic species (Carstens et al., 2004; Fleming and Murray, 2009; Muscarella et al., 2011; Speer et al., 2017). In the case of Erophylla sezekorni, a pollen- and nectar-feeding bat in the family Phyllostomidae, previous work showed a strong pattern of isolation-by-distance between islands (Muscarella et al., 2011). However, in a small portion of its range, E. sezekorni shows population differentiation between 2 groups of islands that are separated by a relatively short distance (∼60 km; Fig. 1). These 2 groups of islands, the Little Bahama Bank and the Great Bahama Bank, are separated by 2 channels—the Northwest and Northeast Providence Channels (collectively referred to as NPC), which have acted as a barrier to dispersal in other bat species (Muscarella et al., 2011; Speer et al., 2017; Fig. 1). While Pleistocene glaciations repeatedly connected and separated islands within the Little Bahama Bank and within the Great Bahama Bank, the NPC are deep channels that remained permanent barriers between the banks through the Pleistocene (see review in Speer et al., 2015).

To examine whether dispersal of E. sezekorni across the NPC occurs rarely or occurs but does not always result in gene flow, we estimated the population structure of E. sezekorni and their associated bat flies, Trichobius frequens. Most bat flies are specific to a single species of bat (Dick and Dittmar, 2014) but do not share an evolutionary history with their hosts (Graciolli and de Carvalho, 2012). This may be due to regional variation in host associations and specificity. In Cuba, T. frequens is commonly collected from Phyllonycteris poeyi and Brachyphylla nana in addition to E. sezekorni. However, bat fly population connectivity is likely not influenced by dispersal of co-roosting bat species in The Bahamas, because P. poeyi and B. nana do not occur on these islands (Speer et al., 2015). In addition, previous research on the host specificity of bat flies indicates there is minimal host switching of bat fly species when multiple species of bat occupy the same roost (Dick and Patterson, 2007; Dick and Dittmar, 2014). To confirm this, we sampled co-roosting species of bats for bat flies and never detected Trichobius frequens on any bat other than E. sezekorni (see Methods).

As with all flies in the superfamily Hippoboscoidea, bat flies deposit partially mature larva on an environmental substrate to pupate and complete their development into adult flies. Female bat flies leave the bat host approximately every 10 days to deposit a larva, which will pupate for 3–4 wk and emerge as an adult fly (Dick and Patterson, 2007; Dick and Dittmar, 2014). This unique life cycle causes bat flies to move frequently between individuals of the same host species that co-occur in a roost. This allows us to test the hypothesis that bats have dispersed, but that this dispersal has not resulted in gene flow. If dispersed bats have occupied the same roost as local bats within a 3–4 wk window, bat flies will be shared between local and dispersed individuals. This would result in low bat population connectivity but high bat fly connectivity. As bats frequently carry more than 1 bat fly on their bodies, a dispersed bat would bring several dispersed flies, allowing us to test patterns of fine-scale bat dispersal that may be difficult to detect without extensive sampling or deep sequencing (i.e., rare or recent dispersal). Many bat fly species cannot fly, and those with functional wings cannot fly well, making it unlikely that bat flies could disperse between islands on their own. In addition, there is strong selection for bat flies to remain on their hosts, because they must feed on blood frequently (Dick and Patterson, 2006).

Parasites impact the structure and stability of the wildlife communities to which they belong (Lafferty et al., 2008), but little is known about the population genetics of wildlife parasites (Huyse et al., 2005). Only 2 previous studies have examined population genetics in bat flies, both conducted on nycteribiid bat flies (Olival et al., 2013; van Schaik et al., 2015). No study has examined the population genetics of the family Streblidae. It is not known how streblid bat fly populations or infestation intensity change through the year based on host life cycle, but both host and parasite life history are expected to impact parasite population genetics (van Schaik et al., 2015; van Schaik and Kerth, 2017). Female and young bats are more heavily parasitized by bat flies than adult male bats (Patterson et al., 2008; Frank et al., 2016), creating the opportunity for bat sex-biased dispersal and roosting patterns to impact bat fly population size and gene flow. As the environment continues to rapidly change, we expect to see greater destabilization of wildlife communities in response to the extinction of host species and secondary extinction of their associated parasites (Strona and Lafferty, 2016). Population connectivity of parasites will help us to better understand fine-scale genetic responses of wildlife communities to environmental change.

Previous research on the phylogeography of bat flies has made use of mitochondrial loci, but it has found low or no population differentiation in these markers (Wilson et al., 2007; Olival et al., 2013). However, many bacterial symbionts, including those that have been documented in bat flies, are known reproductive parasites that facilitate the spread of a single or few mtDNA haplotypes (Jiggins, 2003; Hurst and Jiggins, 2005; Lack et al., 2011). Thus, mtDNA may be a problematic marker with which to estimate population connectivity of bat flies, potentially due to the presence of these bacterial symbionts (Lack et al., 2011). Taking into account previous research that found a low signal of population differentiation in mitochondrial loci in bat flies (Wilson et al., 2007; Olival et al., 2013), we therefore developed the first microsatellites in the family Streblidae to characterize the population connectivity of the bat fly T. frequens. We provide the first examination of population genetics in the family Streblidae and the first examination of bat fly population genetics in the neotropics. This research will provide valuable information about gene flow in parasites, which is not well understood, and may document the consequences of changing environmental landscapes at a fine scale. We propose the use of parasites to test patterns of host dispersal that are difficult to detect using host genetics alone (e.g., dispersal that does not result in mating, or dispersal that is rare or recently started occurring).

We collected tissues and wing punches from a total of 130 individuals of E. sezekorni from 4 islands in The Bahamas: Abaco and Grand Bahama of the Little Bahama Bank and Eleuthera and Long Island of the Great Bahama Bank (Fig. 1; Suppl. Table S1). Bats were sampled from 2 caves on each island using hand nets or mist nets. Two wing punches were collected from each bat using sterile 4-mm biopsy punches. A few bats (<5 per island) were euthanized using halothane. Skin, skeleton, and organ tissue (liver, kidney, spleen, heart, and lungs) samples were curated at the Florida Museum of Natural History and National Museum of The Bahamas (Table S1). In total, 137 T. frequens individuals were collected from E. sezekorni at these localities using featherweight forceps. Upon capture, bats were placed in individual, clean cloth bags to prevent cross-contamination of parasites. All bat tissues and parasites were stored in 95% ethanol at room temperature for 2 wk and then transferred to −20 C. All animal-handling procedures followed guidelines of the University of Florida Institutional Animal Care and Use Committee (2201404427) and guidelines set by the American Society of Mammalogists (Sikes et al., 2011). Bat flies were morphologically identified following Peterson and Hůrka (1974) and Wenzel (1976) using a stereomicroscope (Leica S8 APO, Leica Microsystems, Inc., Buffalo Grove, Illinois). To confirm host specificity of bat flies, we also collected ectoparasites from other species of bats roosting in the same cave. In all but 2 roosts, E. sezekorni was the only species of bat present. Bat flies of Tadarida brasiliensis (84 bat flies from 22 bats) and Macrotus waterhousii (6 bat flies from 2 bats), which were found roosting in the same cave (but not the same bell holes) as E. sezekorni, were examined. Tadarida brasiliensis was parasitized exclusively by a new species of Trichobius, and Macrotus waterhousii was parasitized by Nycterophilia coxata. There was no evidence of host switching in Trichobius frequens. All genetic information is freely available from GenBank (accession numbers MK923981–924037; http://www.ncbi.nlm.nih.gov/bioproject/542721) and FigShare (https://doi.org/10.6084/m9.figshare.c.4561181).

DNA was extracted from one-half wing punch or ∼25 mg of organ tissue using the ZR Genomic DNA Tissue MicroPrep kit (Zymo Research, Irvine, California) and eluted with 50 μl of Elution Buffer. DNA was extracted from bat flies using the ZR Genomic Tissue and Insect MicroPrep kit or ZR Genomic DNA Tissue MicroPrep kit using bead or manual homogenization or puncture of flies followed by proteinase K digestion. Bat flies develop through their first 3 instars inside of the female fly, so extractions from a gravid female will include DNA from the larva. This could cause the presence of 3 microsatellite peaks following genotyping and make it impossible to distinguish the genotype of the mother from that of the offspring. To avoid this, the abdomen of gravid females was removed prior to DNA extraction to prevent amplification of larval DNA during PCR. The exoskeletons of punctured bat flies were retained following digestion of tissue with proteinase K for slide mounting.

Twelve microsatellite primer pairs have already been developed for E. sezekorni (Murray et al., 2008), and we used 8 of these loci. Prior to this study, there were no available microsatellite primers for bat flies of the family Streblidae. We used Roche 454 pyrosequencing and the Primer 3 pipeline to develop 9 polymorphic microsatellite markers for T. frequens (Untergasser et al., 2012; see Supplement for details and Suppl. Table S2). Because bacterial reproductive parasites have been previously implicated in selective sweeps of mtDNA in insects (Duron et al., 2008), including bat flies (Lack et al., 2011), we used the 454 data generated for development of microsatellite loci to data-mine all 16S rRNA bacterial reads using the pipeline mothur v.1.37.4 and identified bacteria to putative OTU using the BLAST 16SMicrobial database (Altschul et al., 1990; Schloss et al., 2009; NCBI Resource Coordinators, 2014). We used seven of these novel bat fly microsatellites for this study. We adapted primers to be used in a 3-primer approach using universal tags to fluorescently label PCR products (Schuelke, 2000; Blacket et al., 2012), following the dye-labeling scheme outlined in Ascunce et al. (2013). PCRs for both E. sezekorni and T. frequens were conducted in 20-μl or 15-μl reaction volumes containing 7.5–8 μl of 2X GoTaq Colorless Master Mix (Promega Corp., Madison, Wisconsin; containing Taq polymerase, 400 μM dNTPs, and 3 mM MgCl₂), final concentrations of 0.025–0.06 μM unlabeled forward primer, 0.25–0.6 μM unlabeled reverse primer, and 0.25–0.6 μM dye-labeled universal primer, diluted with sdH₂O and 1 μl of template DNA to final reaction volume. Thermocycler conditions followed Ascunce et al. (2013), but annealing temperatures followed recommendations by Murray et al. (2008; see Table S2 for T. frequens annealing temperatures). Individuals that failed initial PCR were redone using 2x Type-It Microsatellite PCR Master Mix (Qiagen, Inc., Hilden, Germany) with the same PCR and thermocycler protocols detailed above. All PCRs for loci TRFR 13 and TRFR 38 used Type-It Master Mix.

For genotyping, E. sezekorni PCRs were combined up to 6 plex, and T. frequens PCRs were 4-plexed, with 1.5 μl of FAM-labeled product, 1.5 μl of VIC-labeled product, 3.0 μl of NED-labeled product, and 3.0 μl of PET-labeled product, diluted with sdH₂O up to 50 μl. PCR products were genotyped using the ABI3730xl DNA Analyzer (Applied Biosystems, Foster City, California) at University of Florida ICBR using a GeneScan 600 labeled with Liz dye as a size standard. Raw peak data were visualized and scored using Geneious™ with the microsatellite plugin (Kearse et al., 2012).

We used MICROCHECKER to identify all potential null alleles, removed them from the data set, and estimated F_ST, G”_ST, and D_EST (see supplemental materials for details; Van Oosterhout et al., 2004). Patterns of population structure following removal of hypothesized null alleles matched findings of the complete data set (Tables S3, S4), indicating that homozygote excess in some alleles did not inaccurately bias the signal of differentiation between populations. We used all loci and alleles for the rest of the analyses. In addition, some primer pairs in both species were in linkage disequilibrium (ES17/ES43, p-value = 0.024; TRFR4/TRFR13, p-value = 0.025; and TRFR4/TRFR45 p-value = 0.014; Table I), but none was significant following sequential Bonferroni correction (Rice, 1989).

Cytochrome oxidase subunit II (COII), a mitochondrial gene, was sequenced for 45 T. frequens individuals collected from E. sezekorni and 12 Trichobius sp. 1 individuals collected from Tadarida brasiliensis using primers A-tLEU and B-tLYS (Maekawa et al., 1999). Tadarida brasiliensis co-occurs in roosts with E. sezekorni and has been previously shown to have high population structure across the NPC (Speer et al., 2015, 2017). We included sequences of Trichobius sp. 1 here as an additional examination of the use of mtDNA as a population genetics marker in bat flies. PCRs of COII were conducted in 15-μl reaction volumes containing 7.5 μl of 2x GoTaq Colorless Master Mix (Promega Corp., Madison, Wisconsin) and final concentrations of 0.33 μM of each the forward and reverse primers diluted with 1 μl of template DNA and sdH₂O to the final reaction volume. Thermocycler protocol followed Dittmar de la Cruz and Whiting (2003). Following successful amplification, PCR products were purified using 1 μl of ExoSAP-IT (Affymetrix, Inc., Cleveland, Ohio) per 5 μl of PCR product. Purified PCR products were sent to the University of Florida ICBR for Big Dye Terminator cycle sequencing reactions and analysis on an Applied Biosystems 3730xl Genetic Analyzer (Applied Biosystems). Raw chromatographs were visualized, trimmed, and coupled into consensus sequences using Geneious™ (Kearse et al., 2012). Consensus sequences were aligned using ClustalW with default parameters (Thompson et al., 2002).

We used the R package genepop v4.3 to calculate exact tests of Hardy-Weinberg equilibrium (HWE) and estimate linkage disequilibrium between loci (Rousset, 2008). For each test, we used the Markov chain algorithm with dememorization of 10,000 and 100 batches each consisting of 100,000 iterations. Tests were performed with individuals grouped by island (assuming 4 populations) and with individuals grouped by banks of islands (assuming 2 populations corresponding to the Little Bahama Bank and the Great Bahama Bank) to ensure that population assignment did not bias results.

We calculated G”_ST, Jost's D, and F_ST as measures of population differentiation using GenAlEx v.6.501 (Peakall and Smouse, 2012). Because microsatellites often have a high number of alleles, the maximum possible value of F_ST is not equal to 1, even when there are no shared alleles between 2 populations (Jost, 2008; Meirmans and Hedrick, 2011). All population differentiation statistics were calculated using AMOVA with 9,999 permutations of the data and 9,999 bootstrap replicates. We used the mantel test implemented in GenAlEx v.6.501 to test for the correlation between the natural log of geographic distance (km) and G”_ST/(1 − G”_ST) with 999 random permutations to determine significance. Geographic distance was measured as the shortest geodesic distance between sampling localities.

We used Structure v.2.3.4 (Pritchard et al., 2000) to cluster individuals into populations based on the genotypic and linkage disequilibria between alleles. To prevent incorrect inference of the number of populations based on uneven sampling (Puechmaille, 2016), we subsampled the data 10 times so that each island was subsampled to 27 individuals for E. sezekorni and 21 individuals for T. frequens. Structure was run on each subsampled data set for 500,000 generations following a burn-in of 75,000 under a model of admixture and correlated allele frequencies. The number of clusters was set to K = 1–8, with 10 iterations at each K value. Structure Harvester v.0.6.94 (Dent and vonHoldt, 2011) was used to compare mean lnP(K) and ΔK, from the Evanno method (Evanno et al., 2005), to identify the most likely K value. CLUMPAK was used to compare Structure outputs at each K value and to create figures of Structure outputs (Rosenberg, 2004; Jakobsson and Rosenberg, 2007; Kopelman et al., 2015).

BayesAss v3.0.3 (Wilson and Rannala, 2003) was used to estimate recent (within the last 2 generations) dispersal between all 4 sampled islands and when individuals were split into 2 populations by the NPC. Estimates were conducted on the 10 subsampled data sets used for Structure analysis to reduce effects of variable sample size. In T. frequens data sets, mixing parameters of the migration rate, allele frequencies, and inbreeding coefficients were set to 0.4, 0.6, and 0.7, respectively. For E. sezekorni, mixing parameters for migration rate, allele frequencies, and inbreeding coefficients were set to 0.3, 0.5, and 0.7 in the island-level data set and 0.2, 0.3, and 0.3 when individuals were separated by the NPC. BayesAss was run on all data sets for 50 million generations with a burn-in of 5 million, sampling every 100 generations. Adequacy of mixing, burn-in, and sampling was verified using Tracer v.1.6.0 (Rambaut et al., 2014).

Average gene flow was estimated using Bayesian inference implemented in Migrate-n v.3.6.6 (Beerli, 2006). Bayes factors were used to compare three models of migration (Beerli and Palczewski, 2010): (1) There is 1 panmictic population, (2) migration across the NPC occurs from the Little Bahama Bank to the Great Bahama Bank, and (3) migration across the NPC occurs from the Great Bahama Bank to the Little Bahama Bank. In models 2 and 3, individuals were assumed to be panmictic between islands of the same bank (i.e., between Abaco and Grand Bahama and between Long Island and Eleuthera). Data were coded as repeat numbers, and the Brownian motion model was used to approximate step-wise mutation of microsatellites. We used uniform prior distributions for θ and M (range 0 to 50, window size of 1). Each model was run for 1 billion generations following a burn-in of 25,000 generations, sampling every 1,000 generations and replicated 10 times. We used Metropolis-Hastings sampling and 4 heated chains, following the recommendations in the Migrate-n manual for comparison of Bayes factors.

Microsatellite loci for E. sezekorni yielded an average of 9 alleles per locus, ranging from 5 alleles at ES8 to 17 alleles at ES43 (Table I). Loci for T. frequens had an average of 5 alleles per locus, ranging from 4 alleles at TRFR45 to 8 alleles at TRFR40. Microsatellites genotyped from E. sezekorni had on average 5% missing data (ranging from no missing data for ES17 to 10% missing data in ES8), and genotypes of T. frequens had an average of 4% missing data (ranging from 1% missing data for TRFR4 and TRFR17 to 17% missing data for TRFR43).

Missing data were coded following program-specific guidelines and ignored during analyses. Erophylla sezekorni collected from The Little Bahama Bank had 6 private alleles, and those collected from the Great Bahama Bank had 16 private alleles. There were fewer private alleles detected in T. frequens from the Little Bahama Bank (4 private alleles) and the Great Bahama Bank (5 private alleles). When individuals were split into 4 populations by island, there were private alleles distinguishing the islands in E. sezekorni (Abaco: 2 private alleles, Grand Bahama: 3 private alleles, Eleuthera: 5 private alleles, Long Island: 4 private alleles) and T. frequens (Abaco: 1 private allele, Grand Bahama: 3 private alleles, Eleuthera: 2 private alleles, Long Island: 1 private allele). There was no evidence of differentiation between sampling localities on the same island. Exact tests of HWE suggested that some E. sezekorni loci (ES38, ES43, and ES46) and some T. frequens loci (TRFR4, TRFR13, TRRF37, TRFR43) violated null expectations. Violations of HWE are not surprising, given sampled individuals likely belong to small populations with low gene flow from external sources. However, deviations from HWE may also indicate the presence of null alleles that artificially inflate signatures of population differentiation (Van Oosterhout et al., 2004).

Both E. sezekorni and T. frequens showed population substructure across the NPC (Table II). G”_ST, D_EST, and F_ST showed significant differentiation between the Little Bahama Bank and Great Bahama Bank in both species compared to random, but E. sezekorni populations were much less differentiated than T. frequens populations (Table II). When individuals were split by island, neither species showed differentiation between Abaco and Grand Bahama (Little Bahama Bank; Table III). A Mantel test for matrix correspondence indicated significant correlation between genetic differentiation and geographic distance in both E. sezekorni (correlation coefficient = 0.391, p-value = 0.013) and T. frequens (correlation coefficient = 0.557, p-value = 0.005), in agreement with previous work on the population genetics of E. sezekorni in the broader Caribbean (Muscarella et al., 2011).

Population assignment through the clustering algorithm Structure suggested 2 populations in the bat fly T. frequens, but a single population in the host bat E. sezekorni (Fig. 2). Given the low level of substructuring evidenced by G”_ST, D_EST, and F_ST in E. sezekorni, it is not surprising that there was no evidence of population differentiation using Structure. In T. frequens, lnP(K) was lowest when K = 2 or K = 3, but K = 2 had the highest ΔK across all subsampled data sets (Suppl. Fig. S1). When K was set to 2, parasites north of the NPC generally had a distinct genotype from parasites south of the NPC, in support of estimates of population structuring across this oceanic barrier (Fig. 2; Table II). However, individuals from Eleuthera shared more of their genotype with individuals north of the NPC than did individuals from Long Island. In addition, when K was set to 3, all individuals had a portion of their genotype assigned to cluster 3, but individuals from Long Island tended to have a lower proportion of cluster 3-type alleles, in agreement with evidence for significant isolation-by-distance (Fig. 2). All subsampled data sets of E. sezekorni individuals indicated a single population that spanned the NPC. When K was set to 2–8, all individuals had a portion of their genotypes assigned to each cluster, suggesting a single panmictic population (Fig. 2; Fig. S2).

Estimates of recent dispersal (last 2 generations) support high proportions of non migrants (i.e., residents) in the Little and Great Bahama Banks in both E. sezekorni and T. frequens populations across all subsampled data sets (Fig. 3; Tables S5, S6). When individuals were split by island, recent dispersal was higher between islands on the same bank, and less dispersal occurred across the NPC (Table IV). In both species, low dispersal may be occurring from the Little Bahama Bank to the Great Bahama Bank, although the 95% credibility interval overlaps zero (Table IV; Fig. 3). Long Island showed the highest proportion of non-migrant individuals in both species and is the most geographically distant island. Considering the 95% credibility intervals, which are broad in most estimates, there are no significant differences between the recent dispersal rates of E. sezekorni and T. frequens.

Migrate-n provides estimates of gene flow as an average over the history of populations. If we assume gene flow is proportional to dispersal, these results can be extended to estimate dispersal. We ran Migrate-n for 3 scenarios of migration history and used Bayes factors to determine the most likely model (Beerli and Palczewski, 2010). In both E. sezekorni and T. frequens, models of directional dispersal had a lower marginal likelihood than the single panmictic population model. In E. sezekorni, the most highly supported model was directional dispersal from the Little Bahama Bank to the Great Bahama Bank. In contrast, the most highly supported model in T. frequens was dispersal in the opposite direction (Great Bahama Bank to Little Bahama Bank). Across all models, T. frequens had smaller estimated population size than E. sezekorni (Fig. 3). Long-term gene flow across the NPC was much higher in both species than recent gene flow (T. frequens = 8.10–22.60 and E. sezekorni = 3.63–13.33; Table IV; Fig. 3).

MtDNA sequences for T. frequens and Trichobius sp. 1 showed no variable sites within species. There were 59 fixed polymorphisms differentiating T. frequens from Trichobius sp. 1. Previous studies using mtDNA from Trichobius major also found an absence of variation in 2 other mtDNA genes, cytochrome oxidase I and NADH dehydrogenase 4 (Wilson et al., 2007). The lack of variation in mtDNA may be a consequence of symbiotic bacteria that act as reproductive parasites, like Arsenophonous, Rickettsia, or Wolbachia, all of which have been sequenced from bat flies (Lack et al., 2011; Wilkinson et al., 2016). Wolbachia was detected in pyrosequencing data of T. frequens used for microsatellite development and is known to cause mtDNA sweeps in other insects (Jiggins, 2003; Hurst and Jiggins, 2005; Duron et al., 2008; Table S5).

Our findings indicate low differentiation of E. sezekorni across the NPC. Erophylla sezekorni may not form a completely panmictic population across the NPC, although individuals exhibit much greater connectivity than has been evidenced in other bat species in The Bahamas (Muscarella et al., 2011; Speer et al., 2017). Trichobius frequens has significantly greater population substructuring than their host bats across the NPC but no difference in recent and long-term dispersal rates (Tables II, III; Figs. 2, 3). Both species showed low recent dispersal but high long-term dispersal (Fig. 3). By combining the evidence provided by E. sezekorni and T. frequens, we can confidently reject our hypothesis that bats are dispersing across the NPC but that this dispersal does not result in mating. Despite our initial skepticism that the NPC would act as a barrier to gene flow for E. sezekorni, it appears that recent gene flow across these oceanic channels is low. However, our hypothesis that bat flies act as a fine-scale marker of host bat dispersal is supported. Traditionally, only strict host specialists that occur on a single individual host or share an evolutionary history with a host species have been used to examine host population dynamics. This research indicates that parasites that behave like micropredators, but are reliant on a single host species for dispersal, may also provide valuable information about host population dynamics.

Contrary to our expectations, T. frequens exhibited smaller or comparable effective population size to its host (Fig. 3). This may be due to the greater, although non-significant, population differentiation of T. frequens between islands of the Great Bahama Bank compared to E. sezekorni (Tables III, IV; Fig. 2). Small effective population sizes lead to rapid differentiation in parasites (Huyse et al., 2005) and may make parasite populations more susceptible to the effects of inbreeding depression. Further research is needed to explore the prevalence of small effective population size in parasites and its implications for parasite extinction and host switching.

The findings that bat flies have faster generation times and may have lower effective population sizes than their hosts suggest population genetics of bat flies may be more reflective of fine-scale changes in dispersal of bats than population genetics of the bats themselves (Huyse et al., 2005). The greater population differentiation in T. frequens than E. sezekorni might be explained by 3 scenarios: Bat flies are not dispersing with their host bats (i.e., missing the boat), bat fly population genetics reflect a true pattern of a recent reduction in host dispersal, or bat flies reflect a true pattern of rare periods of high dispersal in their host bats.

The differences between estimated gene flow and direction of gene flow across the NPC for T. frequens and E. sezekorni may indicate that bat flies do not always leave the roost with their host. It has been reported that bat flies leave the host temporarily during their life cycles, including while bats have left the roost to forage (Dick and Patterson, 2007). However, many studies have reliably sampled bat flies from their hosts using mist nets to capture bats outside of the roost (ter Hofstede et al., 2004; Olival et al., 2013). Bat flies can only survive away from their hosts for between 5 and 25 hr (Dick and Patterson, 2006), so there is strong selection for bat flies to stay on the host. The fact that the most well-supported model of dispersal in each species supports dispersal in contradictory directions may indicate that dispersal is happening in both directions across the NPC. We were not able to test this bidirectional model of dispersal with our data set and would require more loci. However, if T. frequens does not disperse with its host, we would expect to see high differentiation between neighboring islands similar to the differentiation observed between banks of islands separated by the NPC. This is not the case; T. frequens showed no significant population structure between Abaco and Grand Bahama (Table III), very low differentiation between Eleuthera and Long Island (Table III), and high recent dispersal between islands of the same bank (Table IV). Based on this evidence, differences in population structure between E. sezekorni and T. frequens cannot be explained by bat flies staying in the roost rather than dispersing with their host.

Alternatively, population structuring of T. frequens may represent a recent decline in dispersal of E. sezekorni. Based on few studies, bat flies have a generation time of approximately 30 days (Dick and Patterson, 2006), while phyllostomid bats have a generation time of about 1 yr (Barclay and Harder, 2003). A decrease in dispersal of the host bat occurring within the last 2 bat generations could potentially have a large effect on bat fly population divergence, given their considerably shorter generation times. Comparison of long-term dispersal estimates to recent dispersal supports this hypothesis, because Migrate-n infers high levels of long-term population connectivity, and BayesAss estimates very low or zero recent dispersal (Table IV; Fig. 3). It is unlikely that the observed pattern is due to founder effects, where E. sezekorni only dispersed across the NPC during 1 short time window in the past. While islands within the Great Bahama Bank and within the Little Bahama Bank were intermittently connected by land or separated by shallow seas during Pleistocene glacial–interglacial cycles, the NPC is deep and did not fluctuate in width during the Pleistocene (for a review, see Speer et al., 2015). Despite the constant presence of the NPC since the formation of The Bahamas, there are sub-fossils of E. sezekorni found on Abaco in the Little Bahama Bank and Andros, the Exumas, New Providence, and San Salvador in the Great Bahama Bank (Morgan, 1989). This suggests that a single narrow window of dispersal of E. sezekorni across the NPC is unlikely, given the low nuclear microsatellite variation. Nuclear sequence data collected from E. sezekorni may further illuminate differences in recent and historical population structuring.

Finally, it is possible that E. sezekorni only rarely disperses across the NPC, but in high numbers. This pattern has been suggested in Cayman Islands bats as a result of hurricanes causing translocation of large numbers of bats to neighboring islands (Fleming and Murray, 2009). The pattern of high average dispersal, low recent dispersal, and greater population structuring in T. frequens matches what we might expect to see in cases of rare bouts of dispersal of many individual bats. While there is currently no direct evidence of storm-mediated dispersal in most islands of the Caribbean, including The Bahamas, we cannot eliminate this possible explanation.

Despite the high level of population differentiation detected in nuclear microsatellite loci of T. frequens, there were no mtDNA polymorphisms among bat flies sampled across the NPC or between islands. Similarly, in another co-roosting species of bat fly, Trichobius sp. 1, there were no polymorphisms detected between islands, but there was high interspecific polymorphism between T. frequens and Trichobius sp. 1. This is potentially due to the bacterial symbiont Wolbachia, which was detected in T. frequens (Table S5). Wolbachia, as an intracellular bacterium, relies on maternal inheritance for transmission (Cariou et al., 2017). To increase its fitness (i.e., transmission), Wolbachia can kill male embryos, change male embryos to female, and cause sterilization of uninfected females through infected males (Cariou et al., 2017). This leads to a decrease in mtDNA polymorphism and an increase in the prevalence of mtDNA lineage(s) associated with the originally infected female flies in a population (Cariou et al., 2017). It is unlikely that this pattern of high nuclear microsatellite differentiation and no mtDNA differentiation is caused by a slow mutation rate at the COII locus in bat flies. Previous research on mtDNA in bat flies and COII in other insects has indicated a rapid mutation rate and intraspecific variation (Maekawa et al., 1999; Witsenburg et al., 2015; van Schaik et al., 2018). The absence of substitutions within species, but high interspecific variation between T. frequens and Trichobius sp. 1, is consistent with observations of mtDNA sweeps induced by bacterial reproductive parasites pervasive in insects (Hurst and Jiggins, 2005). Our data contribute further evidence that mtDNA is problematic for examining phylogeographic patterns in insects due to the effects of symbiotic bacteria that are transmitted intracellularly (Jiggins, 2003; Hurst and Jiggins, 2005).

The NPC may not be a universal barrier to dispersal in bats. Instead, recent changes in roost and foraging habitat availability may impact gene flow of E. sezekorni between the Little and Great Bahama Banks more than factors associated with the NPC. Novel microsatellite markers developed for the bat fly T. frequens, and the first for the family Streblidae, provide evidence of low gene flow across the NPC, supporting findings from E. sezekorni that bat dispersal across the NPC occurs rarely or has recently decreased. Bat flies, and parasites in general, offer an opportunity to test hypotheses of host dispersal that cannot be measured or are difficult to measure using host genetic data alone (e.g., dispersal that does not result in gene flow, or recently or rarely occurring dispersal). Parasites provide an independent source of information about their hosts that may further support or offer novel hypotheses for examining the evolutionary history and ecology of their hosts.

Thanks go to Dr. Susan Perkins for generously providing funds to cover open-access page charges for this manuscript. This research was funded by the Sigma Xi Grants-in-Aid of Research, American Society of Mammalogists Grants-in-Aid of Research, American Museum of Natural History Theodore Roosevelt Memorial Fund, University of Florida Department of Biology John Paul Olowo Memorial Fund and Brian Riewald Memorial Fund, and University of Florida Women's Club Scholarship. We thank Dr. Stuart McDaniel, Adam Payton, Dr. Marina Ascunce, and Dr. J. Angel Soto-Centeno for valuable guidance on project design and data analysis. We thank Samantha Johnson, Kaylee Ludden, Nicole Norelli, Mitul Patel, and Claudia Zafra for assistance in data collection. We credit M. Brock Fenton and Sherri Fenton for photographs of Erophylla sezekorni used in Figures 2 and 3.