Introduction
It has long been known that geographic range varies immensely in size among different species (Brown et al. Reference Brown, Stevens and Kaufman1996). During the last several decades, multiple studies aimed at understanding the mechanisms behind this variation and the determinants of range size have been carried out (see reviews in Brown et al. Reference Brown, Stevens and Kaufman1996; Holt, Reference Holt2003; Gaston, Reference Gaston2003, Reference Gaston2008). In general, it is recognized that a species’ geographic range size (GRS) is the net result of the interplay between extrinsic and intrinsic factors. Among extrinsic factors, environmental limitations (Pigot et al. Reference Pigot, Owens and Orme2010), periodic global climatic changes (Jansson and Dynesius, Reference Jansson and Dynesius2002) and interactions with other species (including competition; Gaston, Reference Gaston2008) affect GRS. From the intrinsic perspective, the GRS of a species may depend on life-history or ecological traits such as body size (e.g. Gaston and Blackburn, Reference Gaston and Blackburn1996), density (Brown, Reference Brown1984) and the breadth of its ecological niche (Eeley and Foley, Reference Eeley and Foley1999). Another intrinsic factor that may play a role in determining a species’ GRS is associated with phylogeny. In other words, the GRS of a species could be, at least in part, a trait inherited from its ancestor(s) as a consequence of the phylogenetic heritability of traits that influence the environmental preferences or environmental tolerance of this species and thus determine its geographic distribution (e.g. Peterson et al. Reference Peterson, Soberón and Sánchez-Cordero1999, but see Mouillot and Gaston, Reference Mouillot and Gaston2007). As a result, closely related (e.g. sister) pairs of species would be more similar in their GRS than pairs of species taken randomly from a phylogenetic tree. This pattern is commonly defined as phylogenetic heritability or phylogenetic signal (PS) in a trait.
Starting from the pioneering work of Jablonski (Reference Jablonski1987), who reported significant positive correlations between the GRSs of closely related species pairs, multiple attempts to test for the effect of phylogeny on GRS have been undertaken (see Zacaï et al. Reference Zacaï2017 for review). These studies produced contradictory findings, from strong PS to none whatsoever (Zacaï et al. Reference Zacaï2017). One reason for these inconsistencies could be differences in methodological approaches (e.g. Webb and Gaston, Reference Webb and Gaston2003 vs. Hunt et al. Reference Hunt, Roy and Jablonski2005), while another reason could be that many studies ignored the role of environmental factors in shaping species ranges. Indeed, when the environment has been taken into account, it was found that phylogenetic heritability can be mediated by spatial and environmental processes (Mouillot and Gaston, Reference Mouillot and Gaston2009; Abellán and Ribera, Reference Abellán and Ribera2011; Zacaï et al. Reference Zacaï2017).
Phylogenetic heritability in GRS has been studied in a variety of extant and fossil taxa (plants and animals, invertebrates and vertebrates, in marine and terrestrial environments). To the best of our knowledge, it has never been studied in parasitic species. In part, this is because the spatial distribution of parasites is closely coupled with the spatial distribution of their hosts (Poulin et al. Reference Poulin2011), so the latter may mask the inherent patterns of the former. Nevertheless, studies of the phylogenetic heritability of parasites’ GRS may help to understand whether the unknown GRS of many species could be predicted from the known GRS of their close relatives (Mouillot and Gaston, Reference Mouillot and Gaston2009). This could be especially important for agents (e.g. helminths) and/or vectors (e.g. ectoparasitic arthropods) of infectious diseases.
PS in the GRS of parasites can be expected because PS in the traits potentially influencing their GRS has been found. Indeed, Krasnov et al. (Reference Krasnov2005a) reported a positive correlation between niche breadth in fleas (e.g. inverse host specificity in terms of the number of host species exploited) and their GRS, whereas the former was found to be phylogenetically heritable (Mouillot et al. Reference Mouillot2006). Further investigations of phylogenetic heritability in the niche breadth of fleas (Krasnov et al. Reference Krasnov, Poulin and Mouillot2011) demonstrated that host specificity, in terms of the number of host species, was phylogenetically heritable in the Palearctic but not in the Nearctic species. This suggests that the PS in these parasites’ GRS could vary geographically and be mediated by either the environment or the history of parasite divergence. On the one hand, phylogenetic heritability in GRS could be modulated by spatial and/or temporal environmental stability with GRS being conserved in stable, but not in variable, environments (Zacaï et al. Reference Zacaï2017). For example, PS in the range size of ammonites during the early Pliensbachian (Lower Jurassic) has been found mainly in physically stable (at that period) open-marine Mediterranean province but not in the less stable Northwest European province that represented a shallow epicontinental sea surrounded by continents (Zacaï et al. Reference Zacaï2017). Moreover, Zacaï et al. (Reference Zacaï2017) did not detect PS of ammonite range size during the last chronozone of the Pliensbachian known for numerous palaeoclimatic and palaeoenvironmental changes. On the other hand, GRS differences between species may accumulate proportionally to the time since divergence (GRS of close relatives are similar at the speciation event but become more and more different later), so that the PS in GRS would decrease with an increase in divergence time (Barraclough et al. Reference Barraclough, Vogler and Harvey1998; Cardillo, Reference Cardillo2015).
The main aim of our study was to infer whether the historical spatiotemporal stability of environment or time since divergence could be the main reason behind PS (if any) in GRS of fleas parasitic on small mammal. To answer this question, we investigated phylogenetic heritability in the GRS of fleas in five regions, each of a different biogeographic realm [Australasian (Australia), Afrotropics (Africa south of the Zambezi River, including South Africa, Namibia, Botswana, Lesotho, Swaziland, Zimbabwe and Southern Mozambique; further referred to as Southern Africa), Nearctic (Canada, Greenland and Alaska; further referred to as Canada), Neotropics (Venezuela) and Palearctic (Mongolia)]. In particular, we asked (a) whether GRS is phylogenetically heritable in fleas and, if yes, (b) whether the manifestation of phylogenetic heritability varies between these regions. We hypothesized that if phylogenetic heritability in the GRS of fleas is modulated by the spatiotemporal stability of the environment, then this heritability will be found in higher levels in the southern regions (Australia, Southern Africa and Venezuela) than in the northern regions (Canada and Mongolia) due to the consequences of Quaternary glaciation which was pronounced much stronger in the Holarctic realm than in the southern realms. This is because glacial/interglacial cyclicity could disrupt geographic ranges of a close relative by opening and closing refugia in an irregular fashion (Andersen and Borns, Reference Andersen and Borns1994) and thus erase any historic signatures. If, however, the manifestation of phylogenetic heritability is associated with time since divergence, then it is expected to be more pronounced in regions with flea faunas comprise a large proportion of species belonging to the most derived lineages (e.g. Ceratophyllidae in Canada and Leptopsyllidae in Mongolia) than in regions with a large proportion of species from the basal lineages (e.g. Stephanocircidae in Australia, Pulicidae and Doratopsyllinae in Southern Africa and Stephanocircidae, Doratopsyllinae and Rhopalopsyllinae in Venezuela) (Zhu et al. Reference Zhu2015).
Materials and methods
Data on GRS
Data on fleas parasitic on small mammals for the five regions and their geographic ranges were taken from Krasnov et al. (Reference Krasnov2005a) which, in turn, were taken from regional monographs (see references in Krasnov et al. Reference Krasnov2005a). Species with cosmopolitan distributions that were likely introduced to many regions by humans, livestock, pets and synanthropic rodents (Xenopsylla cheopis, Nosopsyllus fasciatus) were not considered in the analyses. The size of the geographic range for each species was calculated from a distribution map (see Krasnov et al. Reference Krasnov2005a for details). In brief, a distribution range map for each species was constructed as a polygon map, using ArcGIS Desktop 9.3 software based on published maps and points of occurrence taken from various sources. Then, both regional (across the considered region) and continental (across the respective continent) geographic ranges for each species were calculated. In total, we used data on the geographic ranges of 300 flea species (35 species in Australia, 106 species in Canada, 76 species in Mongolia, 58 species in Southern Africa and 25 species in Venezuela). Regional and continental ranges coincided or almost coincided (less than 10% difference) in all Australian fleas, 17 Canadian fleas, 13 Mongolian fleas, 49 Southern African fleas and three Venezuelan fleas.
Phylogenies
To construct phylogenetic trees for fleas of each region, we used the most recent molecular phylogeny of fleas (Zhu et al. Reference Zhu2015). This tree comprised approximately 10 and 46% (205 and 110, respectively) of the total numbers of species and genera, respectively, of the order Siphonaptera (Medvedev, Reference Medvedev2005). Most genera, but not species, in our dataset were represented by the tree of Zhu et al. (Reference Zhu2015). Consequently, the positions of taxa not represented in the original tree of Zhu et al. (Reference Zhu2015) were determined using their morphologically derived taxonomy (see taxonomic references in Poulin et al. Reference Poulin2006; Krasnov et al. Reference Krasnov, Poulin and Mouillot2011). This was done because the only available information on the majority of flea species is represented by brief descriptions and dichotomous identification keys, so we established most of the within-genus topology according to subgeneric or species group subdivisions based on morphological characters used for identification. All branch lengths were set equal to 1.0, and the trees were arbitrarily ultrametrized using the respective function of the Mesquite Modular System 3.31 (Maddison and Maddison, Reference Maddison and Maddison2017).
Data analysis
Values of GRS were log-transformed prior to the analyses. We estimated PS (the statistical non-independence among species trait values due to their phylogenetic relatedness) using Moran's index of spatial autocorrelation I (Moran, Reference Moran1950), Pagel's (Reference Pagel1999) λ, and Blomberg et al.’s (Reference Blomberg, Garland and Ives2003) K. The three indices are calculated using different approaches.
The Moran's I index (Moran, Reference Moran1948, Reference Moran1950) is the measure of autocorrelation originally used in spatial statistics and adopted for measuring the phylogenetic autocorrelation (Gittleman and Kot, Reference Gittleman and Kot1990). The advantages of using Moran's I include its insensitivity to tree size and tree balance (Lockwood et al. Reference Lockwood2002; Münkemüller et al. Reference Münkemüller2012) and the fact that it applies to the entire phylogenetic tree rather than to pairs of species (Abellán and Ribera, Reference Abellán and Ribera2011). Furthermore, Moran's I is not based on any evolutionary model (Pavoine et al. Reference Pavoine2008; Münkemüller et al. Reference Münkemüller2012). Its value varies from −1 to +1, with positive values indicating that a tested trait is more similar in close relatives than expected by chance and negative values indicating an absence of correlation between the trait values and interspecific phylogenetic distances.
Blomberg et al.’s (Reference Blomberg, Garland and Ives2003) K compares the observed PS in a trait to that under Brownian motion trait evolution. It estimates the strength of the PS as the quotient of the observed and expected mean square error (MSE) ratios, where the observed ratio is the MSE of the tip data measured from the phylogenetically corrected mean, divided by the MSE of the data calculated using the variance–covariance matrix derived from the tree, while the expected ratio is computed from the tree with Brownian motion as the evolutionary process. Higher values of K indicate stronger PS s. For any trait in which Brownian motion is the true evolutionary model (thus implying some degree of phylogenetic heritability), K equals 1, while zero K values (or those that do not differ significantly from zero) indicate no effect of phylogeny, and values >1 suggest strong conservatism of traits (Losos, Reference Losos2008).
Pagel's (Reference Pagel1999) λ tests whether the phylogeny correctly predicts the covariance patterns of a trait among species for a given trait according to a Brownian model of character evolution. Accordingly, λ is a multiplier of the off-diagonal elements of the variance/covariance matrix describing tree topology and branch lengths that transforms the phylogeny to ensure the best fit of trait data to a Brownian motion model. Values of λ range from zero to 1 (but see Münkemüller et al. Reference Münkemüller2012) and gradually eliminate the phylogenetic structure, so that zero values indicate the evolution of the trait to be independent of phylogeny, while values of 1 indicate a Brownian motion model of the trait's evolution on a given phylogenetic tree. Intermediate values of λ correspond to traits evolving according to a process in which the effect of phylogeny is weaker than in the Brownian model (Freckleton et al. Reference Freckleton, Harvey and Pagel2002).
An advantage of Blomberg et al.’s (Reference Blomberg, Garland and Ives2003) K and Pagel's (Reference Pagel1999) λ is that both indices are able to detect a PS when only partial phylogenetic information is available (e.g. when there is no information on branch length, so they are arbitrarily established; Freckleton et al. Reference Freckleton, Harvey and Pagel2002; Blomberg et al. Reference Blomberg, Garland and Ives2003).
All analyses were carried out using the packages ‘phylobase’ (Hackathon et al. Reference Hackathon2017) and ‘phylosignal’ (Keck, Reference Keck2015; Keck et al. Reference Keck2016) implemented in the R Statistical Environment (R Development Core Team, 2017). First, we tested for PS in GRS for each region using the three indices. Then, to test whether PS in GRS varies between clades within a region, we applied an adaptation of the local Moran's I (Local Indicators of Spatial Association; Anselin, Reference Anselin1995; Keck et al. Reference Keck2016) for phylogenetic data [Local Indicators of Phylogenetic Association (LIPA); Keck et al. Reference Keck2016] using function lipaMoran of the R package ‘phylosignal.’ These indicators allow the detection of species with phylogenetic neighbours similar in a given trait and species with phylogenetic neighbours different in a given trait. We calculated LIPAs for each tip in each phylogenetic tree using the proximity matrix based on the number of nodes in order to focus on clades.
The statistical significance of each index was tested for the null hypothesis of GRS values randomly distributed in the phylogeny (i.e. no PS). This was done by randomization for Blomberg et al.’s (Reference Blomberg, Garland and Ives2003) K and Moran's I (shuffling species across the tips of the phylogenetic tree using 999 permutations) and by a likelihood ratio test for Pagel's (Reference Pagel1999) λ as implemented in the function phyloSignal of the package ‘phylosignal’.
Results
The results of testing for PS in the GRS of fleas from the five regions are presented in Table 1. All three indices consistently indicated a significant PS in both regional and continental GRSs of fleas from Canada and in the continental GRS of fleas from Mongolia. Neither index showed significant PS in the GRS of fleas from Australia and Southern Africa. Venezuelan fleas demonstrated significant PS in regional but not continental GRS, but this was indicated by only two of the three indices.
Phylogenetic signal in regional and continental geographic range (GR) sizes of fleas from five geographic regions as measured by Moran's I, Blomberg et al.’s (Reference Blomberg, Garland and Ives2003) K and Pagel's (Reference Pagel1999) λ
The LIPA analyses revealed a significant local positive autocorrelation of GRS in some clades, even in those regions in which significant PS in GRS has not been detected across the entire phylogeny. In Australia, local PS in continental GRS was indicated for pulicid fleas and some Pygiopsylla (Fig. 1). In Canada, this was true for multiple ceratophyllid and hystrichopsyllid clades (Fig. 2). In Mongolia, a local positive phylogenetic autocorrelation in GRS was found for the leptopsyllid genus Paradoxopsyllus and ceratophyllid genera Megabothris, Ceratophyllus and Oropsylla and in two hystrichopsyllid species (Stenoponia singularis and Wagnerina antigua) (Fig. 3). In Southern Africa, the LIPAs of GRS were significant in some, but not other, clades of the large Xenopsylla genus (Fig. 4). Finally, among Venezuelan fleas, only two stephanocircid species (Sphinctopsylla tolmera and Cleopsylla monticola) demonstrated significant LIPAs for regional GRS (Fig. 5).
Local Moran's index values (LIPAs) for continental geographic range size (GRS) of each flea species in Australia. Black circles indicate significant index values. LIPAs allow to detect hotspots of positive and negative phylogenetic autocorrelation, i.e. species with similar neighbours and species with different neighbours.
Local Moran's index values (LIPAs) for continental GRS of each flea species in Canada. Black circles indicate significant index values. See text and legend to Fig. 1 for an explanation of LIPAs.
Local Moran's index values (LIPAs) for the continental GRS of each flea species in Mongolia. Black circles indicate significant index values. See text and legend to Fig. 1 for an explanation of LIPAs.
Local Moran's index values (LIPAs) for the continental GRS of each flea species in Southern Africa. Black circles indicate significant index values. See text and legend to Fig. 1 for an explanation of LIPAs.
Local Moran's index values (LIPAs) for the regional GRS of each flea species in Venezuela. Black circles indicate significant index values. See text and legend to Fig. 1 for an explanation of LIPAs.
Discussion
We found evidence for phylogenetic heritability in GRS in flea assemblages in northern (Canada and Mongolia), but not in southern (Australia, Southern Africa and Venezuela), regions. We also detected strong PS in GRS in some, but not in other, flea clades within a region. Furthermore, the values of the regional Moran's I were predominantly positive (except for the regional GRS in Mongolia), suggesting a tendency of related flea species to be more similar in their GRS than expected by chance. However, this was not the case for the LIPAs, which were both positive and negative, although all significant values were positive.
Phylogenetic heritability of GRS in fleas
Phylogenetic heritability of GRS in fleas suggests that some intrinsic characteristics of species inherited from their common ancestors may affect their geographic distribution. Brown et al. (Reference Brown, Stevens and Kaufman1996) argued that the evolution of GRS in any species can be explained by invoking two, not mutually exclusive, classes of mechanisms. The first class is associated with processes of colonizations and extinctions, whereas the second class includes limitations by environmental factors coupled with the niche requirements of a species, although environment undoubtedly plays a role in the first class of mechanisms as well. Furthermore, the effect of these mechanisms on GRS can be mediated by biotic interactions such as competition (Herrera-Alsina and Villegas-Patraca, Reference Herrera-Alsina and Villegas-Patraca2014). Another type of biotic interaction that may determine a species’ GRS, independently of phylogeny, is a parasitic way of life. Indeed, in contrast to free-living species, the most important section of a parasite's environment is represented by its host. If, for example, an area representing a range of abiotic environmental conditions that a given parasite can tolerate is not occupied by a suitable host species, then populations of this parasite cannot persist in this area. Consequently, the GRS of parasites is, to a great extent, determined by the hosts’ GRS (see Shenbrot et al. Reference Shenbrot, Krasnov and Lu2007 for fleas). As a result, the geographic range of a highly specialized parasite may cover only the geographic range of its specific host, whereas the geographic range of a host-opportunistic parasite may cover the geographic ranges of multiple host species. Krasnov et al. (Reference Krasnov2005a) found that host-specific fleas have more restricted geographic ranges than host-opportunistic fleas. Moreover, flea species with a broad geographic range are not only capable of exploiting more host species, but they also exploit host species from a wider range of taxa. Given the above-mentioned negative correlation between parasite host specificity and GRS (Krasnov et al. Reference Krasnov2005a), phylogenetic heritability in host specificity may thus be the reason behind phylogenetic heritability in GRS. Furthermore, the phylogenetic heritability of parasites’ GRS can be a consequence of the phylogenetic heritability of their hosts’ GRS. This, however, can only be true if (a) the hosts’ GRS is phylogenetically heritable, (b) the parasites are characterized by a high degree of host specificity and (c) the parasites co-speciated with their hosts. Although (a) may be the case (e.g. Carotenuto et al. Reference Carotenuto, Barbera and Raia2010; Cardillo, Reference Cardillo2015), neither (b) nor (c) is true for fleas. Indeed, fleas demonstrate rather loose host specificity with the distribution of the numbers of host species exploited being less right-skewed than in other parasite taxa (Krasnov, Reference Krasnov2008). In addition, the shared evolutionary history of fleas and their hosts was mainly characterized by events other than co-speciation (Krasnov and Shenbrot, Reference Krasnov and Shenbrot2002; Lu and Wu, Reference Lu and Wu2005; but see Hadfield et al. Reference Hadfield2014). Nevertheless, the tight phylogenetically constrained association between parasite and host geographic ranges could be the case for parasites in which host gene flow is an important determinant of a parasite's geographic range (Nuismer and Kirkpatrick, Reference Nuismer and Kirkpatrick2003). This, however, unlikely true for fleas (Lin et al. Reference Lin2014).
Although phylogenetic heritability in the degree of host specificity in fleas has been reported (Mouillot et al. Reference Mouillot2006), further investigation found that this was true mainly in the Palearctic but not in the Nearctic (Krasnov et al. Reference Krasnov, Poulin and Mouillot2011). In this study, we detected PS in the GRS of fleas from both the Palearctic and Nearctic regions. This suggests that a substantial component of GRS heritability may not necessarily be associated with niche heritability and that significant PS in the GRS of species can be found even without differences in biological traits (see also Mouillot and Gaston, Reference Mouillot and Gaston2007).
Whenever significant values of Blomberg et al.’s (Reference Blomberg, Garland and Ives2003) K were found, they were always <1. One of the reasons for K < 1 could be measurement errors in the tree topology or branch lengths (Blomberg et al. Reference Blomberg, Garland and Ives2003). For the flea phylogenies used in this study, branch lengths were unavailable, so they were arbitrarily set to 1, while errors in tree topology were also possible because not all the species used in this study were presented in the molecular phylogenetic tree (Zhu et al. Reference Zhu2015). However, Blomberg et al. (Reference Blomberg, Garland and Ives2003) demonstrated that errors in tree topology, branch length and even the tip data reduce, rather than inflate, the PS. On the one hand, K < 1 in our study indicates that the effects of phylogeny on GRS evolution in fleas are weaker than expected from a Brownian motion model (Freckleton et al. Reference Freckleton, Harvey and Pagel2002; Blomberg et al. Reference Blomberg, Garland and Ives2003), thus suggesting substantial evolutionary lability in flea GRS. The most likely reason for this is that the evolution of GRS in a flea species may be strongly affected by a variety of factors. In particular, the geographic distribution of fleas has been shown to be determined predominantly by abiotic environmental factors, whereas the effect of host distribution has been found to be relatively weak (Krasnov et al. Reference Krasnov2015). Adaptation to a particular environmental factor could occur only in some, rather than all, species in a clade. For example, two closely related Xenopsylla species demonstrated substantial differences in their preferences for ambient temperature and relative humidity (Krasnov et al. Reference Krasnov2001). On the other hand, uncertainties in phylogeny could be the reason for the lack of PS in the GRS of fleas from southern regions. Nevertheless, we believe that this is unlikely due to at least two reasons. First, the support for flea tree topology up to generic level is relatively high (Zhu et al. Reference Zhu2015). Second, it is reasonable to assume that uncertainties in tree topology are similar among trees constructed for different regions, but the significant PS was detected in Canada and Mongolia only. Therefore, biological processes are the most likely reason behind the above-mentioned lack of PS of GRS.
In addition, interspecific differences in colonization and dispersal abilities between sister species could be a factor erasing PS in GRS of fleas. It is generally unknown what is a predominant mode of speciation in the order Siphonaptera. However, sympatric speciation is thought to be more common in parasites than in free-living species (Théron and Combes, Reference Théron and Combes1995; McCoy, Reference McCoy2003). In fact, intrahost speciation has been suggested for fleas based on the ratio between the number of species and the number of genera especially in warmer regions (Krasnov et al. Reference Krasnov2005b). Consequently, flea sister species at the time of speciation could have approximately equal geographic ranges. If colonization and/or dispersal abilities differ between these species for some reason, then sizes of their future ranges would also substantially differ. The feasibility of this scenario is supported by a recent study by van der Mescht et al. (Reference van der Mescht, Matthee and Matthee2015). They found that two very closely related flea species (Dinopsyllus ellobius and Dinopsyllus abaris) that diverged as recently as during the Upper Miocene, have differential dispersal abilities and geographic ranges of different size.
Time since divergence, sampling artefacts and historical accidents
The size of a species’ geographic range is, to a great extent, determined by the breadth of its ecological niche (Brown, Reference Brown1984; see Slatyer et al. Reference Slatyer, Hirst and Sexton2013 for recent support). It has been shown that the ecological niches of sister species diverge little at or close to the speciation event, while niche differences appear to accumulate later (Peterson et al. Reference Peterson, Soberón and Sánchez-Cordero1999). Consequently, in older clades, the accumulation of GRS differences among closely related species over time can erase the PS. Our results partly support this statement. Indeed, we found PS in the GRS of fleas from Mongolia and Canada only. In these regions, about 60% of flea assemblages are represented by species belonging to the two youngest flea families Leptopsyllidae and Ceratophyllidae. These families likely diverged as late as during the Eocene, whereas the first divergence event in the order was that of Macropsyllidae, which already began in the Lower Cretaceous (Zhu et al. Reference Zhu2015). Species from the most basal flea taxa either were not present in Mongolia and Canada at all (Macropsyllidae, Stephanocircidae) or their proportion was low (Pulicidae, Doratopsyllinae; 8% in Mongolia and 4% in Canada).
In regions where no PS of GRS was detected, flea assemblages were, to a great extent, composed of basal flea taxa. Australian flea fauna was represented by species of Macropsyllidae, Stephanocircidae, Pygiopsyllomorpha, Doratopsyllinae and Pulicidae (of them, the youngest taxon is Pygiopsyllomorpha that diverged during the Upper Paleocene-Lower Eocene). The majority of Southern African fleas were pulicids, chimaeropsyllids and craneopsyllids, whereas younger taxa constituted only ca. 10% of the flea fauna. The youngest flea family (Ceratophillidae) was represented in Venezuela by only two species (Plusaetis dolens and P. smiti; ca. 8% of the entire flea assemblage).
The Local Indicators of Phylogenetic Association (LIPAs) revealed that in each region there were some clades within which the geographic ranges of species were more similar than expected by chance. Interestingly, this was found in both younger and older clades. For example, 38 species in Canada and 15 species in Mongolia had phylogenetic neighbours with geographic ranges of similar sizes. Of them, 20 Canadian and 13 Mongolian species belonged to the young families Leptopsyllidae and Ceratophyllidae. However, 18 Canadian species with significant LIPAs were hystrichopsyllids. Hystrichopsyllidae (sensu Hopkins and Rothschild, Reference Hopkins and Rothschild1962, Reference Hopkins and Rothschild1966 and Medvedev, Reference Medvedev1998) comprise several subfamilies, some of which have been proven to be paraphyletic and diverged at different times (Whiting et al. Reference Whiting2008; Zhu et al. Reference Zhu2015). Among 18 Canadian hystrichopsyllids with significant LIPAs of GRS, 13 belonged to the derived clade Rhadinopsyllinae 2 + Anomiopsyllinae 2 (sensu Zhu et al. Reference Zhu2015).
However, in other regions, significant LIPAs of GRS were unexpectedly found in basal taxa. For example, C. monticola and S. tolmera belonging to the ancient Stephanocircidae were the only species with significant LIPAs in Venezuela. This may be the result of a sampling artefact because these two species have a predominantly Andean distribution with Venezuela constituting only a small portion of their continental range (Tipton and Machado-Allison, Reference Tipton and Machado-Allison1972). As a consequence, the LIPAs of these species’ GRSs were significant for the regional (Fig. 5) but not the continental GRS (P > 0.1 for both). Another example is represented by significant LIPAs of GRS in pulicid fleas (Echidnophaga species and Xenopsylla vexabilis) in Australia. The occurrence of these fleas in Australia is considered to be a relatively recent phenomenon following historical accidents of their introduction by far-ranging hosts (birds and bats for Echidnophaga and commensal rats for X. vexabilis) and secondary adaptations to marsupials (Traub, Reference Traub, Traub and Starke1980). In other cases, significant LIPAs of GRS in ancient clades, such as those in two clades of African Xenopsylla (Fig. 4), are more difficult to explain. Nevertheless, it can be speculated that they are a consequence of a relatively recent burst of speciation in the genus. Indeed, Xenopsylla is one of the most speciose flea genera. It originated in Africa (Traub, Reference Traub, Traub and Starke1980; Whiting et al. Reference Whiting2008), although the basal beginning of the family Pulicidae is argued to be Neotropical (Zhu et al. Reference Zhu2015). Furthermore, it is possible that the radiation of African Xenopsylla (in particular, of the hirsuta group sensu Hopkins and Rothschild, Reference Hopkins and Rothschild1953) coincided with the radiation of their main hosts (the gerbilline clade Gerbilliscus + Gerbillurus) that started in the Upper Miocene (Alhajeri et al. Reference Alhajeri, Hunt and Steppan2015), whereas the family Pulicidae (to which Xenopsylla belongs) originated much earlier (during the Paleocene; Zhu et al. Reference Zhu2015).
In conclusion, GRSs in fleas may, to some extent, be determined by their phylogeny via, for example, the phylogenetic heritability of niche breadth (but not only). The manifestation of this heritability may be associated with the time of divergence. However, a pattern of PS in GRS can be distorted due to sampling artefacts and historical accidents. Nevertheless, patterns found in this study may not necessarily be true for other taxa of arthropod ectoparasites. Patterns and mechanisms behind the relationships between phylogeny and GRS in more host specific (e.g. lice) or more host opportunistic (e.g. ticks) ectoparasites could differ from those found in fleas. Further studies are needed to elucidate any general rule.
Acknowledgements
We thank four anonymous referees for helpful comments on an earlier version of the manuscript. This is publication no. 957 of the Mitrani Department of Desert Ecology.
Financial Support
This study was partly supported by the Israel Science Foundation (grant no. 149/17 to B.R.K. and I.S.K.). LVDM and EMW received financial support from the Blaustein Center for Scientific Cooperation (Jacob Blaustein Institutes for Desert Research, Ben-Gurion University of the Negev) (LVDM and EMW) and the French Associates Institute for Agriculture and Biotechnology of Drylands (Jacob Blaustein Institutes for Desert Research, Ben-Gurion University of the Negev) (LVDM).
Conflicts of Interest
None.
Ethical Standards
Not applicable.
