Comparative phylogeography between two generalist flea species reveal a complex interaction between parasite life history and host vicariance: parasite-host association matters

Background In parasitic taxa, life history traits such as microhabitat preference and host specificity can result in differential evolutionary responses to similar abiotic events. The present study investigates the influence of vicariance and host association on the genetic structure of two generalist flea species, Listropsylla agrippinae, and Chiastopsylla rossi. The taxa differ in the time spent on the host (predominantly fur vs. nest) and level of host specificity. Results A total of 1056 small mammals were brushed to collect 315 fleas originating from 20 geographically distinct localities in South Africa. Phylogeographic genetic structure of L. agrippinae and C. rossi were determined by making use of 315 mitochondrial COII and 174 nuclear EF1-α sequences. Both parasites show significant genetic differentiation among the majority of the sampling sites confirming limited dispersal ability for fleas. The generalist fur flea with a narrower host range, L. agrippinae, displayed geographic mtDNA spatial genetic structure at the regional scale and this pattern is congruent with host vicariance. The dating of the divergence between the L. agrippinae geographic clades co-insides with paleoclimatic changes in the region approximately 5.27 Ma and this provides some evidence for a co-evolutionary scenario. In contrast, the more host opportunistic nest flea, C. rossi, showed a higher level of mtDNA and nDNA spatial genetic structure at the inter-populational scale, most likely attributed to comparatively higher restrictions to dispersal. Conclusions In the present study, the evolutionary history of the flea species could best be explained by the association between parasite and host (time spent on the host). The phylogeographic pattern of the fur flea with a narrower host range correspond to host spatial genetic structures, while the pattern in the host opportunistic nest flea correspond to higher genetic divergences between sampling localities that may also be associated with higher effective population sizes. These findings suggest that genetic exchange among localities are most likely explained by differences in the dispersal abilities and life histories of the flea species. Electronic supplementary material The online version of this article (doi:10.1186/s12862-015-0389-y) contains supplementary material, which is available to authorized users.


Background
Paleoclimatic events have been put forward as one of the main mechanisms causing speciations and extinctions of evolutionary lineages [1][2][3][4]. The differences in ecology and life history of taxa, however, often result in differential species-specific responses to similar climatic events [5][6][7] and a simple explanation describing speciation processes hardly exists. For parasitic taxa the situation is further perplexed by complex life cycles, host biogeography, and the level of parasite-host associations [8][9][10][11][12][13][14][15][16]. Although some progress has been made regarding this topic [17][18][19][20] more quantifiable data are needed to make more accurate predictions on for example the factors affecting range expansions and connectivity among populations. For many parasitic taxa, this is desperately needed from a disease ecology perspective [21,22].
In southern Africa, several phylogeographic studies have been performed on small vertebrate taxa (for example see [23][24][25][26][27][28][29]). From these studies, paleoclimatic changes are once again suggested as one of the main drivers of evolution. Congruent phylogeographic patterns among lizards, rodents, shrews and elephant shrews also suggested the existence of vicariant biogeographic barriers in the region (for review see [7,30]). Of particular relevance to the present study is the support for regional vicariance found in Rhabdomys [29], Myosorex [27], Otomys [28] and Micaelamys [26]. These small mammals show distinct genetic clades that can be associated with the xeric western succulent biomes and the more mesic eastern predominantly grassland biome of southern Africa ( Fig. 1a; Additional file 1; also see [7]). In addition, Rhabdomys, Myosorex and Micaelamys show differentiation to a larger or lesser extent between the winter rainfall zone and the aseasonal rainfall zone in the Cape Floristic Region ( Fig. 1a; Additional file 1; also see [30]). From these studies it is evident that the diversification of small mammals in the region is driven by complex interactions between life history, climate and topographic barriers.
The effect of vicariant barriers on the evolution of parasites occurring in southern Africa is virtually unknown. An exception to this is provided by [13], who recently indicated partial congruence in spatial genetic structure between Rhabdomys and its species specific lice, Polyplax. Distinct genetic clades in the parasites and in the hosts supported the documented mesic eastern and xeric western divide in the subregion (see above). Co-evolution analyses, however, failed to support a strong signal of geographic co-differentiation between parasite and host and the authors suggested that the resultant pattern is due to the synergistic effects of parasite traits (host specificity), host-related factors (the vagility and social behavior of Rhabdomys) and the biogeography (vicariance) of the region [13]. To expand on these findings, we here present phylogeographic data on two generalist flea species occurring on Rhabdomys and other small mammal species in southern Africa.
Individual flea species are normally adapted to utilize a variety of hosts in the environment and are all obligatory blood feeders of endothermic vertebrates [31,32]. The life cycle consists of four stages (egg, larva, pupa, and adult) and fleas are dependent on the nest environments of the hosts to variable extent [31]. Immatures develop entirely in the off-host environment while adult stages have to spend some time on the host to obtain a blood meal. The length of time that is spent on the host, however, varies between flea taxa [31][32][33]. This observation has provided an opportunity to coarsely classify fleas according to their microhabitat preference as either a "fur" (adult stage spend more time on the host) or a "nest" flea (adults spend more time in the nest/burrow of the host; [31][32][33]). These differences can have profound effects on the ability of the flea species to disperse over the landscape and we therefore predict that fur fleas will show more phylogeographic congruence with the host/s (adult fleas spend more time on the host and can thus disperse over the landscape) and nest fleas will show more restricted gene flow across the landscape (since all juveniles and the adults spend the majority of their life cycle off the host and have limited dispersal capabilities).
To test the hypothesis that fleas with different levels of host association will differ in spatial genetic structure, the fur flea, Listropsylla agrippinae, and the nest flea, Chiastopsylla rossi [34,35] were selected. Both L. agrippinae and C. rossi occur widespread throughout southern Africa ( Fig. 1) [35,36] and span at least two well documented vicariant biogeographic barriers in the region ( Fig. 1; [7,30]). Apart from differences in the duration of time adults spend on the host, the two species also differ in the level of host specificity (niche breath). Listropsylla agrippinae has two principle host taxa (Rhabdomys spp. and Myotomys unisulcatus) whereas C. rossi has at least four main host taxa (Rhabdomys spp., Otomys irroratus, M. unisulcatus and Tatera brantsii) documented [34,35,37]. It has been suggested that host specificity may influence the level of intraspecific genetic divergences since more generalist parasite species will show a higher level of intraspecific genetic variation enabling them to infest a broader host range [38][39][40]. Furthermore, microhabitat preference (specifically referring to fur vs nest) is not always related to host specificity, meaning that fleas may differ in host specificity irrespective of microhabitat preference (see [32,34,36]). If niche breath and time spent on the host affect 1) genetic diversity of the parasite and 2) congruence in phylogeographic structure between parasites and hosts [38][39][40], we predict that L. agrippinae, when compared to C. rossi, will show lower levels of intraspecific genetic diversity and also more similarities to the vicariant patterns of the hosts they infest. Genetic distances among individuals from different sampling sites are expected to be higher for the nest flea, C. rossi, when compared to the fur flea, L. agrippinae, since the latter can utilize the host for dispersion. By considering the interplay between life history and geography, the present study should add valuable information needed to explain some of the evolutionary processes that shape ectoparasite distribution and diversity.

Parasite prevalence and distribution
Both R. pumilio and R. dilectus were trapped at a contact zone, Fort Beaufort (also see [29]), and represent the only  Table 5. The vicariant breaks in South Africa (1) winter and aseasonal rainfall break and (2) xeric and mesic biomes are indicated by white dashed lines, b the mapped distribution of L. agrippinae and c C. rossi are taken directly from [34] van der Mescht et al. BMC Evolutionary Biology (2015) 15:105 locality where more than one Rhabdomys species was trapped (see Additional file 2). Listropsylla agrippinae were recorded at lower prevalence at the majority of localities when compared to C. rossi ( Fig. 2; Additional file 2). Listropsylla agrippinae was found on six host species and was most prevalent on R. intermedius (23.02 %) (Fig. 2a). Chiastopsylla rossi was found on nine host species, but most prevalent on R. intermedius (38.89 %) and O. irroratus (37.31 %) compared to the other 4 host taxa (Fig. 2b). Individuals selected for sequencing were representative of the various host species trapped at each locality ( Fig. 2; Additional file 2).

Characteristics of molecular markers
Mitochondrial COII data for 126 L. agrippinae and 189 C. rossi specimens were generated and attempts were made to generate a similar nuclear EF1-α data set. Nuclear amplification was less successful but we nonetheless generated 94 L. agrippinae and 254 C. rossi nuclear alleles (Table 1) [GenBank: KR 263182-263844]. Haplotype and nucleotide diversity values were lower for L. agrippinae when compared to C. rossi for mitochondrial COII and nuclear EF1-α data (Table 1). Unexpectedly, the nuclear EF1-α data for both species resulted in higher diversity estimates in terms of haplotypic and nucleotide diversity when compared to the mitochondrial COII data (Table 1).

Phylogeographic analyses Listropsylla agrippinae
TCS statistical parsimony analysis based on the mitochondrial COII data showed two distinct clusters (L1 and L2; Fig. 3a and b) that could not be connected within the 95 % confidence interval. These two clusters are separated by an average corrected mtDNA sequence divergence of 3.02 % (±0.36) and this large differentiation is further supported by the Neighbor-Net tree (Fig. 3c). The clusters correspond to the xeric western (L1) and eastern mesic (L2) zones of the country. The BAPS analysis corroborated the two clusters detected by TCS and the Neighbor-Net tree, and also indicate some additional substructure within cluster L1 (P = 1.00; log likelihood of optimal partition = − 1179.80) ( Fig. 3a and c). These two subclusters loosely correspond to the winter and aseasonal rainfall divide (Fig. 1). The spatial genetic patterns were less pronounced when the nuclear EF1-α data is compared to the mitochondrial DNA network (Figs. 3 and 4). There is some evidence for regional clustering as indicated by the grouping of similarly shaded colors in the parsimony network and Neighbor-Net tree ( Fig. 4) but more importantly, however, the most common DNA allele is shared between the far eastern and far western side of the sampling distribution. Mantel tests for mitochondrial COII data indicated that there were weak isolation by distance when calculated for all localities (r = 0.11; P = 0.00) and when calculated for each cluster L1 (r = 0.11; P = 0.01) and L2 (r = 0.16; P = 0.01) separately. The time of divergence between the two main clusters within L. agrippinae date back to the late Miocene which was estimated to be 5.27 Ma (95 % HPD interval: 2.31, 9.58 Ma).
Genetic landscape interpolation surface plots for the COII data indicated that there were marked differences between some sampling localities when the distribution of genetic diversity found at each is compared across the landscape (Fig. 5a). The graphical representation suggests that the area of greatest genetic diversity was found in the contact zone close to the 'Bedford gap' (Fig. 5a). Fixation index values of mitochondrial COII for L. agrippinae were significant at all levels and the highest level of differentiation was recovered among subclusters (60.00 % of variation) ( Table 2). Fixation index values of nuclear EF1-α for L. agrippinae were significant at all levels but in this instance the highest level of variation was recovered within localities (63.26 % of variation) ( Table 2). Pairwise Φst values among sampling localities showed that the majority of sampling sites are significantly  Summary statistics for the mitochondrial DNA (COII) and nuclear intron (EF1-α) markers sequenced from L. agrippinae and C. rossi. The number of sequences/alleles, total number of haplotypes, total number of singleton haplotypes, fragment length (bp), number of polymorphic sites (P), nucleotide diversity (π) and haplotype diversity (h) is indicated   differentiated from each other for both the nuclear and mtDNA data (Table 3).

Chiastopsylla rossi
TCS analysis based on mitochondrial COII data showed a remarkably different genetic pattern for C. rossi when compared to L. agrippinae. The majority of the localities sampled are characterized by unique divergent haplotypes ( Fig. 6a and b), and although two distinct genetic clusters separated by an average sequence divergence of 3.01 % (±0.32) were also obtained for this species (C1 and C2; Fig. 6b), it was not congruent with vicariant host patterns (Fig. 6). Instead, limited haplotype sharing was evident among distant sampling sites and the majority of populations were characterized by unique/closely related locality specific haplotypes ( Fig. 6a and b). The lack of clear geographic structure is further supported by the Neighbor-Net tree (Fig. 6c) and moreover by the BAPS analyses suggesting at least six distinct subclusters present within C. rossi (P = 0.99; log likelihood of optimal partition =−1862.85). The spatial genetic pattern of the nuclear EF1-α network and the Neighbor-Net tree suggest more lineage sharing across the landscape when compared to the mitochondrial DNA COII data (Figs. 6 and 7). Similarly to L. agrippinae, the nuclear data also indicate some level of population differentiation (closely related haplotypes confined to single localities). The nuclear TCS analyses suggest a distinct clade comprising individuals from KP, MN and AF and this assemblage is supported as one of the six mtDNA BAPS subclusters. This subcluster is mainly found in the north eastern part of South Africa, and although it may have some biological meaning, the Neighbor-Net analysis (Fig. 7c) show that comparatively the level of differentiation is low. The majority of the subclusters identified by the mtDNA and nDNA data do not corresponded to any of the known biogeographic breaks previously documented for the majority of the host species and in fact are scattered across the landscape (Fig. 6a). Mantel tests indicated that there were virtually no isolation by distance  Pairwise Φst values among L. agrippinae sampled localities for mitochondrial DNA (COII) below and nuclear intron (EF1-α) above the diagonal. Significant values (p < 0.05) are highlighted in bold. (Locality codes as in Table 5) within C. rossi when calculated for all localities (r = 0.06; P = 0.00). The scattered distribution coupled to high interpopulational genetic variation in C. rossi is supported when the genetic landscape surface plot is considered. In contrast to L. agrippinae, C. rossi show several genetic discontinuities (peaks) throughout the total sampling distribution (Fig. 5b). Fixation index values of the mitochondrial COII for C. rossi were significant at all levels and the highest level of differentiation was recovered among subclusters (51.29 % of variation; Table 2). The highest level of differentiation for the EF1-α was recovered within localities (48.79 % of variation; Table 2). Pairwise Φst values for the mitochondrial and nuclear data also showed significant differentiation among the majority of sampling localities and the highest level of differentiation was again generally detected between KP, MN and AF and the remainder of the populations (Table 3).

Discussion
Parasites are generally considered to have a high mutation rate, small effective population size and a limited dispersal ability. From this, it is predicted that pronounced spatial genetic structuring will be evident among sampling sites due to reduced gene flow and increased genetic drift [8,[41][42][43]. Overall, the results presented for the two generalist parasite species in this study conform to these suggestions. Most of the geographically distinct populations sampled show significant differentiation among sampling sites at both the mtDNA and nDNA level (Tables 3 and 4) a scenario attributed to reduced gene flow across the landscape (also see [38,44]).
Fleas are considered generalist parasites and are thus predicted to show very little phylogeographic congruence with host genetic structure [38,39,44]. The present study provides the first evidence for the converse and, by making use of a comparative approach between two species, Fig. 6 Chiastopsylla rossi mitochondrial COII. a A map of geographical distribution of sampling localities in South Africa for C. rossi. Bayesian analysis of population structure (BAPS) revealed six subclusters indicated by different colours. b Statistical parsimony mitochondrial COII haplotype network colour coded according to sampling locality. Circle size depicts frequency, branches depict single mutational steps, small black circles display intersections and cross hatching indicate missing haplotypes/mutational steps. Clusters are bound by black (C1) and deep yellow (C2) line boxes. c Neighbour-Net phylogenetic network for C. rossi labelled according to haplotype groupings from statistical parsimony results and are bound by black (C1) and deep yellow (C2) line boxes highlights the importance of the level of association between parasite and host in shaping genetic diversity across the landscape. The fur flea (L. agrippinae) which seems to have a narrower niche breath (i.e. narrower host range) [34,35,37], and spend longer time periods on the host, show three distinct mtDNA phylogeographic clades which are markedly congruent to the previously published regional vicariant biogeographic regions based on the patterns obtained in Rhabdomys [29], Micaelamys [26], Otomys [28] and Myosorex [27] (see Additional file 1). Further support for the notion that the large scale dispersal of L. agrippinae is closely dependent on host movement can be obtained from the observation that the 95 % confidence interval of the time of divergence between the two major clades found in L. agrippinae (5.27 Ma; 95 % HPD interval: 2.31, 9.58 Ma) corresponds reasonably well with the timing of the divergence of the majority of rodent host lineages in the region [7]. In contrast, the more host opportunistic nest flea (C. rossi) show, in agreement with other phylogeographic studies on fleas [38,39], virtually no congruence in phylogeographic patterns between host and parasite. The distinct clustering of C. rossi fleas sampled in the north east of the country (KP, MN, AF: Fig. 6a and Fig. 7) is not clearly depicted by longer branches in the Neighbor-Net analyses (Fig. 7c) and also not congruent with well documented biogeographic provinces. We argue that although this finding is interesting, it may simply be an artefact of the sampling distances between the various sampling sites.
The documented contact zone between R. pumilio and R. dilectus at Fort Beaufort ( [29]; Fig. 1; Additional file 1) provides an interesting scenario to further explore host specificity among the two flea species (no contact zones have been described for the other hosts species sampled in this study). In the zone of contact, the two Rhabdomys species harbored distinct L. agrippinae mtDNA lineages represented by L1 (found exclusively on R. pumilio) and L2 (found exclusively on R. dilectus;  (Figs. 6 and 7). Rhabdomys species is considered the principle hosts of L. agrippinae and M. unisulcatus is considered to be an auxiliary host, whereas C. rossi seems to have a larger host range [34,35,37]. Albeit based on a very small sample size, the tighter host association displayed by L. agrippinae might explain why we found distinct mtDNA lineages for L. agrippinae on Rhabdomys spp., but no differentiation in C. rossi in the contact zone.
When the intraspecific phylogeographic structures of L. agrippinae and C. rossi are compared to each other some more discrepancies are evident. First, our data provide additional support for the hypothesis that more generalist parasites will show higher levels of genetic diversity when compared to more specific parasites ( [16]; Table 1). The connectivity among sampling sites is also markedly different between the two flea species. The more host specific fur flea L. agrippinae, show lower inter-populational divergence within clades (Fig. 3) when compared to the more host opportunistic nest flea, C. rossi (Fig. 6). The same trend is evident when the nDNA data are considered (Figs. 4 and 7). The higher level of inter-populational divergences of C. rossi are best reflected by the large number of site changes among locality specific haplotypes (Fig. 6) and the numerous peaks on the landscape interpolation plots (Fig. 5b). These findings could be interpreted to be the direct result of differences in the dispersal abilities of the two fleas and can best be ascribed to differences in host association. Both parasites show isolation by distance but in both cases the correlation is extremely weak (indicated by the r values). In the case of the host specific fur flea, L. agrippinae, host movement provide hitchhiking possibilities within clades resulting in a higher level of connectivity among sampling sites. In the host opportunistic nest flea, C. rossi, the parasite has less opportunity to spread via host movement, resulting in a pattern of distantly related haplotypes at most sites. At specific geographic sites, it is possible that populations experience high levels of genetic drift that is homogenizing the locality specific signals [8,[41][42][43]. Long distance dispersals of the host opportunistic nest flea studied herein is a rare event and happen in a random way utilizing a wide niche breath (in the absence of strong isolation by distance), culminating in a pattern where some shared haplotypes are found on opposite ends of the geographic scale (for example see haplotype sharing between the group of localities with different shades of blue (FB, AL, TC, HB, DE) and geographically distant localities (BR and also DF; Fig. 6).
From our study it was also evident that L. agrippinae occurred at a lower prevalence than C. rossi overall and the pattern was also apparent at most localities where  the distribution of the two species overlap. The same trend in prevalence and abundance for L. agrippinae and C. rossi was recorded in previous studies in South Africa [35,37,45]. The two flea species differ in terms of body length with L. agrippinae being larger (on average 3650 μm) compared to C. rossi (on average 1750 μm) ( [46]; also supported by body size index in [35]). Studies on free-living taxa and more recently on ectoparasitic mites of small mammals recorded a negative relationship between body size and abundance [47][48][49]. In the latter study it was suggested that the pattern may be due to higher host-induced mortality (grooming) associated with larger bodied ectoparasites [49]. Listropsylla agrippinae is almost twice the size of C. rossi and it is possible that host grooming (allo-and autogrooming) resulted in lower L. agrippinae prevalence. Interestingly, in the present study the difference in body size between the two flea species seem to support their level of host specificity. Studies on microbes and diatoms and more recently on ectoparasitic mites recorded a negative relationship between body size and niche breadth [49][50][51]. It appears that smaller bodied species are more adaptable to environmental fluctuations which results in a larger niche breath or in the case of parasites a larger host range. This is in contrast to larger species that have a smaller niche breath/host range due to narrower tolerance levels [49][50][51].
If prevalence is positively correlated to abundance in fleas [52,53], then it could be argued that the higher abundance (and niche breath) of the more host opportunistic nest flea, C. rossi, on the hosts will facilitate dispersal among sampling sites. This pattern is however contrary to what we found in the present study where C. rossi is more structured between sampling sites when compared to L. agrippinae. One possible explanation for this may relate to differences in the effective population sizes between the two flea species. A higher genetic diversity, as found in C. rossi, is expected for populations with larger effective population sizes [41,44] and it is furthermore reasonable to predict that nest fleas will predominantly have higher abundances in the nests of their hosts [32,39]. Given the nest bound nature of C. rossi, the number of dispersing individuals on the hosts may thus not be large enough to overcome the effects of drift when introduced into the nests of the hosts elsewhere (containing a large population of local genotypes in the new environment). In contrast, for L. agrippinae, the comparatively lower effective population sizes can facilitate the signature of more haplotype sharing among localities within clades.

Conclusions
In the light of re-emerging flea borne diseases worldwide, it is important to have a thorough understanding of the mechanisms that are involved in shaping flea distribution and movement [22]. From our study it is evident that host association (microhabitat preference and host specificity) plays an important role in flea dispersal and subsequent gene flow within and between geographic locations. This study also provides the first evidence of congruent phylogeographic vicariant patterns between a generalist parasitic flea and its hosts. Unfortunately our conclusions are based on a single study comprising two distinct life histories. More parasite taxa and gene fragments needs to be evaluated in order to formulate stronger hypotheses in this regard.

Sample collection
Small mammal trapping was performed during 2010-2013 at 20 localities (natural areas and low density grazing farms) in South Africa ( Fig. 1a; Table 5). Baited Sherman-type live traps were set in a line transect and sampling varied between 4 to 7 days per locality. All adult specimens of species that are listed as potential hosts for the two flea species were selected and juveniles were released at the trap site. Trapped animals were placed in a plastic bag before they were euthanized with sodium pentobarbital (200 mg/kg; ethical approval reference number SU-ACUM11-00004). The bodies were brushed over a white plastic tray and all fleas were collected. The brush was inspected and cleaned (using 96 % ethanol) after each animal was processed and new brushes were used for each host species at each locality. Individual fleas were placed in separate tubes filled with 96 % ethanol. Before DNA extraction, L. agrippinae and C. rossi individuals were preliminary identified using a Leica stereoscopic microscope (Leica Microsystems, Wetzlar, Germany) and the taxonomic key of [34]. After DNA extraction the exoskeletons of all extracted fleas were mounted (see [35]) and a thorough morphological identification was done under a Leica DM 3000 light microscope (Leica Microsystems, Wetzlar, Germany) using the key of [34]. Most of the host species that were trapped during the study are quite common and widely distributed throughout South Africa. As a result voucher specimens of the species are readily available in several museums. In the case of the fleas, voucher specimens will be deposited in the Museum of the Department of Conservation Ecology and Entomology (Stellenbosch University) and the National Flea Collection in Johannesburg (both in South Africa).

DNA extraction, amplification and sequencing
Total genomic DNA was extracted with a Qaigen, DNeasy® Blood and Tissue kit (Qaigen, Valencia, CA, USA) following the protocol of the manufacturer. Whole flea specimens were placed in the extraction buffer containing Proteinase K (600 mAU/ml solution or 40 mAU/mg protein), and digested at 56°C overnight. After digestion, flea exoskeletons were removed for identification purposes (see above).

Alignment and phylogenetic analyses
Sequences for each gene fragment were edited and aligned using BioEdit Sequence Alignment editor 7.2.5 [56]. Mitochondrial sequences were translated into amino acids using EMBOSStranseq (www.ebi.ac.uk/Tools/st/ emboss_transeq) to confirm functionality. Heterozygous positions in the nuclear fragments were resolved in DNASP 5 [57] using PHASE 2.1.1 [58,59]. The algorithm was run for 1000 generations with a thinning interval of 1 and burn-in of 100 generations. Phases with a 0.9 probability or higher were considered resolved and the analysis was performed three times to see if there was any significant difference in results between runs [58].

Diversity indices and population level analyses
Nucleotide diversity (π) and haplotype diversity (h) values were obtained using DNASP 5 [57]. The evolutionary relationships between haplotypes were investigated by constructing statistical parsimony networks with 95 % confidence intervals in TCS 1.21 [60]. The best-fit models of sequence evolution were determined for each fragment under the AICc [61,62] in jModelTest 2.1.4 [63,64]. To further gain an evolutionary perspective on the association between clusters, individual HKY-corrected networks of each flea species were drawn for mitochondrial COII and nuclear EF1-α using the Neighbor-Net method [65] implemented in SplitsTree 4.5 [66]. HKY-corrected sequence distances among mitochondrial COII and nuclear EF1-α were calculated using PAUP* v4.0b10 [67]. To investigate population structure without a-priori assumptions, a Bayesian Analyses of Population Structure (BAPS) was performed in BAPS 6.0 [68] on the mtDNA data. Spatial genetic mixture analyses of individuals and of groups were performed independently using a vector of maximum K values (each replicated 5 times; [68,69]). Locality information from where specimens were obtained for each of the two flea species, with codes corresponding to Fig. 1 Analysis of molecular variance (AMOVA; [70] and pairwise Φst statistics between sampled populations were performed in ARLEQUIN 3.5.1.2 [71]. Only sampling localities with more than 5 individuals were included. The AMOVA higher level group differentiations were defined based on the subclusters obtained in the BAPS analysis (only sampling localities with more than 5 individuals were used). Mantel tests [72] were performed to test for isolation by distance in ALLELES IN SPACE [73]. Spatial genetic structure was further explored by making use of genetic landscape shape interpolation surface plots constructed in ALLELES IN SPACE [73].

Dating main phylogenetic events
We used a relaxed exponential Bayesian molecular clock as implemented in BEAST 2.1.3 [74] to estimate divergence time between clusters. Siphonaptera lacks a useable fossil record and as a calibration point we employed the 2.3 % per million years estimated for various arthropod taxa [75,76]. The HKY + G model was used and the birth-death process of speciation with exponential priors specified. MCMC simulation ran for 20 million generations, sampling every 10 000 generations for each of the two runs performed. Convergence and mixing were assessed in Tracer 1.6 [74] and the first 25 % of the trees were discarded as burn-in. A maximum clade credibility tree was produced in TreeAnnotator 2.1.3 [74].

Availability of supporting data
DNA sequence data are available in GenBank [GenBank: KR 263182-263844].
Additional file 2: Host identity, abundance and parasite prevalence. Host identity and abundance, parasite prevalence and the number of specimens sequenced for each flea species per locality. *Could not determine which host the samples were sequenced from.
Additional file 3: Primer and PCR specifications.docx. Primers used for PCR amplification of mitochondrial and nuclear genes for the two flea species. COII amplification consisted of a denaturation cycle of 1 min at 95°C followed by a 10 cycle loop of 1 min at 95°C, 45°C, and 72°C, respectively. A 30 cycle loop was then performed with denaturation for 1 min at 93°C followed by annealing for 1 min at the primer specific temperature, and 1 min extension at 72°C, followed by a final extension period of 5 min at 72°C. General PCR cycling conditions for the EF1-α region included an initial denaturation of 5 min at 94°C followed by 40 cycles of 30 s denaturation at 94°C, 45 s annealing at primer specific temperature, and 1 min extension at 72°C, followed by a final extension period of 7 min at 72°C. *Primers from [54].

Competing interests
The authors declare that they have no competing interests.

Authors' contributions
This study forms part of the PhD research of L.V.D.M who conducted laboratory and field work, as well as the bulk of the data analyses. L.V.D.M. conceptualized some of the ideas and produced the first draft of the manuscript. C.A.M and S.M. conceptualized and initiated the study, provided samples and supervision during the field and laboratory work and significantly contributed towards finalizing the manuscript.

Authors' information
Luther van der Mescht studied functional and evolutionary ecology of fleas in South Africa as part of his PhD research as a member of the Department of Conservation Ecology and Entomology and Evolutionary Genomics Group at Stellenbosch University. Conrad A. Matthee is a professor and member of the Evolutionary Genomics Group at Stellenbosch University. His research interests include systematics and phylogeography of southern African taxa. Sonja Matthee is a senior lecturer and researcher in the Department of Conservation Ecology and Entomology at Stellenbosch University. Her research interests include parasite taxonomy and parasite-and disease ecology.