Specimen collection

We collected bat flies from three species of flying fox (Pteropus hypomelanus, P. vampyrus, and P. lylei) from localities across Southeast Asia (Figure 1, Table 1). In Malaysia, canopy mist nets were deployed near diurnal roosting trees or nearby feeding sites, and bats were immediately removed upon capture and held individually in cloth bags for ~1 hr before sampling. Bats were anesthetized using medetomidine/ketamine (0.025/2.5 mg/kg) administered intramuscularly or were restrained without anesthesia to collect ectoparasites and wing biopsies for host DNA [35]. Bat flies were collected with forceps by examining the pelage of anesthetized or restrained bats, and placed directly in 95% ethanol for preservation. All bats were handled and sampled in accordance with IACUC protocol (#AAAA3272) from Columbia University. GPS coordinates were recorded for all sampling localities at the point of capture, although in Vietnam and Cambodia local hunters captured bats alive in mist nets near feeding sites and the exact roost localities were not known but likely within 50 km of the sampling site. All bats were marked, photographed, and released at the site of capture except for individuals of P. lylei and P. vampyrus from Vietnam which were collected as voucher specimens currently accessioned at the Institute for Ecology and Biological Resources, Hanoi, Vietnam. From July 2004-May 2007, a total of 125 bat flies were collected from 66 individuals of the three host species (Table 1). Fly specimens are currently accessioned at Western Kentucky University, State University of New York at Buffalo, and the American Museum of Natural History. Sampling was not evenly distributed across host species; 104 flies were collected from three populations of P. hypomelanus, 17 flies from P. vampyrus, and four flies from P. lylei.

Laboratory methods

DNA was extracted from one randomly selected fly per bat using a technique where 1–3 legs per fly were removed (allowing retention of fly voucher). Qiagen DNeasy tissue extraction kits were used per manufacturer’s protocol with a 24 hr tissue digestion, and two combined elutions of 50ul of Buffer AE warmed to 55°C.

Bat flies were sequenced at three mtDNA genes, cytochrome B (cytB), cytochrome oxidase II (CoII), and cytochrome oxidase I (CoI), that were expected to show genetic differences at the population level based on previous intraspecific studies of Diptera [36–39]. Published forward and reverse primers used to amplify each gene were as follows: for cytB, A5 and B 1.1 [36]; for CoII, A-tLEU and B-tLYS [40]; and for CoI LepF1 and LepR1 [41]. For the cytB and CoII genes, we used PuReTaq Ready-To-Go™ PCR Beads (GE Healthcare) with 21 μl of molecular grade ddH₂0, 1 μl of each primer [10 mM], and 2 μl of template DNA. PCR conditions for cytB were identical to those in Dittmar and Whiting (2003), and for the CoII PCR conditions were an initial denaturation period of 3 min at 94°C, followed by 35 cycles of 94°C for 1min, 47°C for 1min, and 74 C for 1min, with a final extension at 74°C for 7min. PCR mix and cycling conditions for CoI followed Hebert et al.[41] although total reaction volume was reduced in half to 25 μl. Attempts to amplify the mtDNA control region were unsuccessful using primers and conditions from Oliveira et al.[42], multiple stutter bands were amplified and a single PCR product was never obtained for the control region. Negative controls were always included in PCR reactions to assess possible contamination.

PCR products were cleaned using Agencourt (Beverly, Massachusetts) AmPure magnetic beads, cycle sequenced with Big Dye v.3.0 terminator mix (Applied Biosystems, Inc., Foster City, California), and final DNA was precipitated using Agencourt CleanSeq magnetic beads. All sequencing was performed on an ABI 3730xl capillary sequencer (Applied Biosystems) at the Sackler Institute for Comparative Genomics at the American Museum of Natural History, New York. Sequences were edited in Sequencher v.4.6 and manually corrected for ambiguous base calls. Alignment of sequences was done using MAFFT v.6 using default parameters [43]. Aligned sequences were trimmed using MacClade 4.08 [44]. No insertions or deletions were present in sequences of any of the three genes.

Phylogenetic and population genetic analyses

Phylogenetic relationships and nucleotide/haplotype diversity were initially determined for each mtDNA gene fragment independently. As there was a paucity of variable and parsimony-informative sites for each marker alone, sequence data from the three mtDNA genes were concatenated for a total of 1744 bp for each specimen. Phylogenetic analyses on the pooled dataset were conducted using PAUP* for maximum parsimony (MP) analysis [45], RAxML v.7.0 for maximum likelihood (ML) analysis [46], and Mr. Bayes for Bayesian analysis [47]. Nodal support was evaluated with the nonparametric bootstrap method [48]. MP bootstrap analysis used 1000 replicates, TBR branch swapping, with a starting tree obtained by random stepwise addition and additional sequences added at random with 10 replicates. The optimum model for ML analysis was determined using Akaike’s Information Criterion to be GTR + G using Modeltest v.3.7 [49]. Optimal evolutionary models did not differ by locus, and data were not partitioned. RAxML analysis of 1000 bootstrap runs beginning with a random seed, and user-defined outgroup, was implemented on the CIPRES Portal webserver found at http://www.phylo.org/sub_sections/portal/. Bayesian analysis was conducted using MrBayes version 3.2 [47], with data partitioned by gene. The GTR + gamma + i model of evolution was used with a flat prior for topologies, uniform priors on gamma and alpha parameters, an unconstrained exponential prior on branch lengths, and Dirichlet priors on all other parameters. A total of ten million generations were run, and convergence was assessed using AWTY [50]. Sequences obtained from Eucampsipoda sundaica (#KB5) were used as an outgroup to Cyclopodia horsfieldi in all phylogenetic analyses because of its basal position within the Cyclopodiinae [14].

Aligned FASTA files were collapsed into variable sites and haplotypes for parsimony network reconstruction with the online tool FaBox 1.31 (http://users-birc.au.dk/biopv/php/fabox/). Fifty-one bat flies were sequenced, but outgroup and individual flies that did not have complete sequences for all three genes were excluded in subsequent analyses (PER08, PER02, PT15, PT12, and KB1), leaving 45 individuals. Statistical parsimony networks, useful for inferring relationships among sequences that have recently diverged, were created using TCS v.1.21 [51]. Parsimony networks were explored using the default connection limit of 95%, but due to one divergent haplotype in the network (RC24), final network was determined with a user-defined limit set to 19 steps. Haplotype (Hd) and nucleotide (π) diversity [52] were calculated in DnaSP v.4.5 [53]. Pairwise F_ST values were calculated in Arlequin using a method based on pairwise differences in sequence data, and statistical significance was assessed with 1000 permutations [54]. We tested for isolation by distance between C. horsfieldi populations from Malaysia only, using a Mantel test [55] with 1000 permutations implemented in the R package adegenet 1.3-6 [56]. All sequences are available on GenBank: CytB (KF273687 - KF273736); CoI (KF273737 - KF273783); and CoII (KF273784 - KF273833).

Analysis of molecular variance (AMOVA), was used to test hypotheses regarding the partitioning of genetic variation for C. horsfieldi among host species, among sampling localities within host species, and among individuals within sampling localities [57]. AMOVA was conducted in Arlequin v.3.1 [58] using the standard haplotype format with statistical significance assessed by 1000 permutations.

Bat fly morphology and demography

All bat fly specimens examined were all identified morphologically as Cyclopodia horsfieldi (Figure 2) with the exception of two individuals of Eucampsipoda sundaica[59] collected from a Pteropus vampyrus (#0709041) from Kuala Berang, Malaysia (Table 1). Some intraspecific morphological variation was noted, including minor differences in counts of ctenidial spines on male sternite and counts of dorsal abdominal setae on the females, but no systematic character variation related to host species or geography was observed. The sex ratio of C. horsfieldi specimens examined was male-biased, 1.85♂ to 1♀. We observed that most bats were parasitized by at least one fly, and that Pteropus hypomelanus individuals had higher numbers of parasites than P. vampyrus. Several P. vampyrus captured hosted no flies but this was very rarely observed for P. hypomelanus. Many P. hypomelanus examined harbored 4+ flies per individual (Figure 3); and some individuals had 10+ flies. Unfortunately, bats were not exhaustively sampled for flies and the number of fly specimens collected was generally limited to a few flies per bat regardless of parasite load, and quantitative data on intensity of infestation was not collected. The number of bat flies collected from each host was not significantly different (at the p<0.05 level) between P. vampyrus and P. hypomelanus using both a standard t-test (t=−1.85; df=13.8; p-value=0.085) and Mann–Whitney test (W=176.5; p-value=0.066). We also inadvertently collected Neolaelaps spinosa (Acari: Mesostigmata) mites as we removed individual flies with forceps. These mites appear to be phoretic with bat flies [60], a relationship that has been documented with mites and other dipteran species [61–63].

Mitochondrial sequence variation, Cyclopodia horsfieldi

Population genetic structure, Cyclopodia horsfieldi

Phylogenetic analyses of the combined mtDNA data set using ML (Figure 4) and MP (Figure 5) showed little resolution of population genetic structure for C. horsfieldi sampled across sites in Southeast Asia. No significant population structure that corresponded to either geography or host species was detected in any of the analyses (Figures 4, 6, 7). Despite an almost complete lack of resolution, some flies from geographically distant areas formed clades with greater than 50% support in the MP tree (e.g. PT13, KB4, and CM3, from Peninsular Malaysia and Vietnam), also evident in the Bayesian tree (Figure 6). Bat flies collected from Pualu Pangkor (PPH), off the west coast of Peninsular Malaysia, appear to be distinct in the haplotype network (Figure 7, yellow), but several individuals shared identical haplotype sequences with flies collected 300 km away on the east coast of Malaysia from Pulau Perhentian (PER) (Figure 7). The most genetically divergent bat fly individual, RC24, was a male morphologically indistinguishable from all the other Cyclopodia horsfieldi specimens.

Host specificity

Morphological and molecular species identification revealed that one fly species, Cyclopodia horsfieldi, parasitized all three Pteropus host species sampled from Malaysia, Vietnam and Cambodia. The only exceptions were two individuals of Eucampsipoda sundaica collected from an individual Pteropus vampyrus in Peninsular Malaysia. Minor morphological differences were observed among individuals of Cyclopodia horsfieldi, but these did not correlate with sampling locality or host species, and the lack of systematic morphological variation is supported by a similar lack of geographic or host-specific population genetic structure. The sex ratio of the bat flies we randomly sampled was male-biased. This observation agrees with previous reports of excess males in bat fly populations, and may be caused by grooming mortality patterns [70] or the diurnal activity of female flies during larviposition [71].

Experimental studies have shown that some bat flies exhibit behavioral preference towards their primary host over congeneric and confamilial species [72]. However, we observe a natural case of parasite oligoxeny among three bat species in Southeast Asia, which has also been documented in Cyclopodia spp. from flying foxes in Australia [16]. Oligoxeny is more likely in hosts that share geographic and habitat niches and/or co-roost at the same sites, and it appears to be an evolutionary advantageous strategy for an ectoparasite [2]. Our results suggest that co-mingling of Pteropus spp. in Southeast Asia is more common than previously assumed. On occasion, Pteropus vampyrus/P.lylei and P. vampyrus/P. hypomelanus co-roosting has been observed in the same trees in Ca Mau, Vietnam and Pulau Langkawi, Malaysia, respectively [69]; and similar observations have been made in Thailand (P. Duengkae, per. comm.). It is also possible that these species could be sharing bat fly parasites without simultaneous occupation of the same roost, i.e. sequential use of a roosting site within a 2–3 week window where flies may emerge from metamorphosis on the roost substrate.

The presence of Eucampsipoda sundaica flies found on a Pteropus vampyrus individual may represent a case of host-switching, which would suggest contact and ecological overlap between other fruit bat species in Peninsular Malaysia. Eucampsipoda sundaica is the sole ectoparasite of the Dawn bat, Eonycteris spelaea[24], although there are previous records of this species from Cynopterus sphinx in India, Pteropus in Myanmar, and Rousettus amplexicatus from the Philippines [59, 73]. Also, we found one genetically divergent bat fly, RC24, to be an outlier in haplotype network and phylogenetic analyses. This fly was morphologically indistinguishable from the other Cyclopodia horsfieldi examined, but could possibly represent a male of the sister species, C. sykesii[16]. C. sykesii is primarily associated with Pteropus giganteus from South Asia [16], and thus could possibly represent a case of host switching by secondary contact between these flying fox species. More fly specimens from P. vampyrus should be sampled and examined to confirm the rarity of these results.

Comparative host-parasite population structure

Overall, we found that populations of Cyclopodia horsfieldi lacked population genetic structure across geographically distant sites and host species in Southeast Asia. However, pairwise F_ST values between some fly populations, particularly from Pteropus hypomelanus, corroborate expectations of reduced gene flow based on the population genetics of their flying fox hosts [69]. Olival [69] found that island populations of P. hypomelanus were significantly differentiated at mtDNA markers from East to West in Malaysia (e.g. F_ST = 0.95, mtDNA control region), but had much higher levels of gene flow between the east coast islands (F_ST =0.0 to 0.4). Here we observe a similar pattern for the parasite, Cyclopodia horsfieldi, with significant pairwise F_ST values between flies from Pulau Pangkor off the west coast and Pulau Tioman (F_ST =0.560) and Pulau Perhentian (F_ST =0.435) off the east coast of Malaysia. In contrast, F_ST values among flies from the two east coast islands were low and not significant (F_ST =0.031). We also observed low and non-significant F_ST values between mainland populations of Pteropus vampyrus and island populations of P. hypomelanus, with the exception of Pulau Pangkor (PPH). These data suggest ongoing or recent gene flow among Cyclopodia horsfieldi parasites between mainland Pteropus vampyrus and eastern P. hypomelanus; and less contact between P. vampyrus and the small western island of Pulau Pangkor.

In contrast, Pteropus vampyrus was found to be essentially panmictic across a very large geographic range of thousands of kilometers in Southeast Asia at mtDNA markers [69], and these data were corroborated by satellite telemetry showing regular long-distance dispersal (100 s of kilometers), lack of roost fidelity, and large home range sizes (128,000 km²) across Malaysia, Thailand, and Indonesia [74]. The shallow branching pattern in our ML phylogeny and mixing of haplotypes observed in the statistical parsimony network suggest that Cyclopodia horsfieldi populations have recently diverged or are subject to ongoing gene flow. Frequent contact between flying fox host species and subsequent ectoparasite gene flow may best explain the lack of parasite population structure observed. In particular, we suggest that the highly volant Pteropus vampyrus may be acting as a “vector” spreading bat flies to other conspecific populations and Pteropus species in the region during long-distance dispersal events. Satellite telemetry studies have shown that P. vampyrus uses small islands as stopover sites when migrating to and from Peninsular Malaysia and Sumatra [74]. We also observed higher nucleotide and haplotype diversity values for flies from P. vampyrus relative to those from P. hypomelanus. This suggests that populations of Cyclopodia horsfieldi from Pteropus vampyrus may have larger effective population sizes and may be acting as a source for founding island populations of parasites [12, 75]. In summary, our observations of parasite population structure, combined with prior results from host population genetics and satellite telemetry, lend support to the idea that Pteropus vampyrus is actively dispersing parasites to the outlying island populations of P. hypomelanus and making contact, and not vice versa.

Three alternative explanations, beyond host-mediated gene flow of Cyclopodia horsfieldi, may explain our results. First, the observed lack of geographic or host species population genetic structure seen in C. horsfieldi could simply be due to invariability or insufficient variance of the molecular markers examined. We suggest this is not the case as a number of previous studies of winged, free-living dipterans using the same markers at similar geographic scales have found significant geographic population structure [36–39, 76]. All else being equal, free-living, volant flies should have lower levels of population genetic subdivision than parasites almost wholly dependent on their hosts for dispersal (i.e. bat flies). Second, endosymbionts, e.g. Wolbachia or Arsenophonus, may have influenced the observed lack of population structure and genetic variation via mitochondrial sweeps, as seen in other dipteran species [34, 39, 77, 78]. Differences in Wolbachia infection and immunity can create striking differences in mtDNA diversity and lead to speciation, e.g. in fig wasps [79]. A number of novel endosymbiont lineages have been identified in bat flies [80–82], and previous studies have suggested that selective sweeps caused by these endosymbionts may explain a lack of mtDNA diversity in bat flies [34]. While we cannot rule out this possibility, it seems unlikely that this selective pressure has influenced the demographic history of each bat fly population sampled. Also, even in the case of Wolbachia infection in insects, high levels of migration may still be a more prominent factor reducing genetic differentiation in species with potential to disperse long distances [83]. Third, demographic factors, i.e. populations bottlenecks with subsequent expansion, could also explain a lack of genetic variation and phylogeographic structure in Cyclopodia horsfieldi[84]. This scenario also seems unlikely, as it also would have to occur independently across multiple geographic localities for each sampled population. The association between Cyclopodia spp. and their flying fox hosts is likely not recent [16], and we believe that high gene flow among parasites is the most parsimonious explanation for the observed results.