Introduction
It is commonly known that parasite evolution is tightly linked to the evolution of their hosts. More than 100 years ago, Heinrich Fahrenholz (Reference Fahrenholz1913) stated that the phylogenies of parasites and their hosts are expected to be congruent if speciation in parasite lineages is a response to speciation in the lineages of their hosts (= cospeciation = codiversification; see glossaries in the study by Althoff et al., Reference Althoff, Segraves and Johnson2014; Blasco-Costa et al., Reference Blasco-Costa, Hayward, Poulin and Balbuena2021). This idea was further formulated as Fahrenholz's Rule (Dogiel, Reference Dogiel1941; Stammer, Reference Stammer1957; Brooks, Reference Brooks1979; Poulin, Reference Poulin2007). However, the lack of appropriate analytical methods did not allow for explicitly testing various hypotheses of parasite–host coevolution until the late 1980s. Then, numerous studies of cophylogeny (i.e. the study of concordance between the phylogenies of 2 groups of interacting species; Blasco-Costa et al., Reference Blasco-Costa, Hayward, Poulin and Balbuena2021) in many host–parasite systems were carried out (e.g. Hafner and Nadler, Reference Hafner and Nadler1988, Reference Hafner and Nadler1990; Paterson et al., Reference Paterson, Gray and Wallis1993; Hafner and Page, Reference Hafner and Page1995; Beveridge and Chilton, Reference Beveridge and Chilton2001; Desdevises et al., Reference Desdevises, Morand, Jousson and Legendre2002; Krasnov and Shenbrot, Reference Krasnov and Shenbrot2002). The majority of these studies have demonstrated a general lack of full congruence of parasite and host phylogenies, thus suggesting that the history of parasites and their hosts is complicated by the predominance of various evolutionary events other than codiversification.
Fleas (Insecta: Siphonaptera) are a convenient model taxon for cophylogenetic studies. This is because (a) all members of this relatively small (about 2500 species) order are obligate haematophages, (b) they are monophyletic (Whiting et al., Reference Whiting, Whiting, Hastriter and Dittmar2008; Zhu et al., Reference Zhu, Hastriter, Whiting and Dittmar2015), and (c) the majority of their hosts belong to a few orders of a single vertebrate class (Mammalia) (only about 6% of species parasitize birds; Medvedev and Krasnov, Reference Medvedev, Krasnov, Morand, Krasnov and Poulin2006). Although earlier narrative descriptions suggested a tight coevolution between fleas and their hosts (Traub, Reference Traub1972, Reference Traub, Traub and Starcke1980, Reference Traub and Kim1985), cophylogenetic studies provided contradictory results. Krasnov and Shenbrot (Reference Krasnov and Shenbrot2002) demonstrated that cospeciation was generally absent in the evolutionary history of associations between fleas and their dipodid (jerboas) hosts, whereas frequent host switching and linear sorting (extinction of a parasite from a host lineage after a cospeciation event) events were widespread (see also Lu and Wu, Reference Lu and Wu2005 for leptopsyllid fleas and ochotonid hosts). However, mapping of host associations on flea phylogenetic trees showed predominant associations of certain flea phylogenetic lineages with certain host phylogenetic lineages (Whiting et al., Reference Whiting, Whiting, Hastriter and Dittmar2008). Moreover, at the large taxonomic scale, when a variety of flea and host species were considered in the analyses, the phylogenetic interaction effect appeared to be substantial and indicated that related flea species were more likely to be found on related host species (Hadfield et al., Reference Hadfield, Krasnov, Poulin and Nakagawa2014; Krasnov et al., Reference Krasnov, Shenbrot, Khokhlova and Degen2016). This was the case despite the fact that the degree of phylogenetic incongruence in hosts was concentrated in some, but not other, phylogenetic clades, whereas it was characteristic for the entire phylogenetic tree of fleas (Llaberia-Robledillo et al., Reference Llaberia-Robledillo, Lucas-Lledó, Pérez-Escobar, Krasnov and Balbuena2023). These contradictions could arise due to the substantial influence of ecological and geographic factors on flea-host common history (Krasnov and Shenbrot, Reference Krasnov and Shenbrot2002), with the effect of these factors varying between biogeographic realms (Krasnov et al., Reference Krasnov, Shenbrot and Khokhlova2022a). The results of the above-cited studies also suggest that phylogenetic congruence might arise from the matching of the deep, rather than the shallow, nodes of the interactors' phylogenetic trees (Blasco-Costa et al., Reference Blasco-Costa, Hayward, Poulin and Balbuena2021). Furthermore, congruent phylogenies could result from processes other than strict codiversification (Althoff et al., Reference Althoff, Segraves and Johnson2014; Russo et al., Reference Russo, Miller, Tooker, Bjornstad and Shea2018). For example, phylogenetic tracking is a pattern in which speciation events in one lineage (parasites) follow speciation events in another lineage (hosts), given that the former strongly depend on the latter but not vice versa, whereas vicariance results in congruent phylogenies if both interactors are subject to similar abiotic isolation events within the same geographic region so that the isolation of populations in both interactors causes parallel branching of their trees. Blasco-Costa et al. (Reference Blasco-Costa, Hayward, Poulin and Balbuena2021) argued that all these processes, including codiversification, are not mutually exclusive, and their combinations may produce patterns observed in nature.
The phylogenetic patterns of species associations have recently started to become an additional tool for understanding species assembly in communities of interacting partners (Corro et al., Reference Corro, Villalobos, Lira-Noriega, Guevara, Guimarães and Dáttilo2021; Krasnov et al., Reference Krasnov, Shenbrot and Khokhlova2022a; Blasco-Costa et al., Reference Blasco-Costa, Hayward, Poulin and Balbuena2021 for review). In addition, combining biogeography and cophylogeny will allow a better understanding of the common history of these interactions because the geographical location of present species and ancestral nodes in their phylogenies might indicate hotspots of various coevolutionary events (Berry et al., Reference Berry, Chevenet, Doyon and Jousselin2018; Blasco-Costa et al., Reference Blasco-Costa, Hayward, Poulin and Balbuena2021). Another way to infer the relationship between geography and cophylogeny is to detect geographic regions accounting for the independent diversifications of either parasite or host lineages and to test for the similarity of these regions' geographic positions for parasites and their hosts. In other words, this approach involves biogeographic regionalization incorporating phylogenetic information (Holt et al., Reference Holt, Lessard, Borregaard, Fritz, Araújo, Dimitrov, Fabre, Graham, Graves, Jønsson, Nogués-Bravo, Wang, Whittaker, Fjeldså and Rahbek2013; Falaschi et al., Reference Falaschi, Marta, Parrino, Roll, Meiri and Ficetola2023). Recently, Maestri and Duarte (Reference Maestri and Duarte2020) proposed the concept of evoregions, namely biogeographic regions where the majority of species evolved from one or a few ancestors inhabiting these regions. They also introduced the concept of evolutionary transition zones as regions characterized by high phylogenetic turnover. Furthermore, Maestri and Duarte (Reference Maestri and Duarte2020) proposed an approach for classifying evoregions and evolutionary transition zones based on fuzzy logic. An advantage of using a fuzzy approach is that it allows representing the boundaries between regions as gradients (Olivero et al., Reference Olivero, Márquez and Real2013).
Although fleas have been thought to originate in Gondwana (Zhu et al., Reference Zhu, Hastriter, Whiting and Dittmar2015), their diversification outburst likely happened in the Palearctic (Medvedev, Reference Medvedev1996, Reference Medvedev2005), resulting in the Palearctic having the highest flea species richness (about 900) as compared to other biogeographic realms (Krasnov, Reference Krasnov2008). This high diversity could stem from 2 not mutually exclusive processes, namely in situ diversification (Morrone and Gutiérrez, Reference Morrone and Gutiérrez2005) and dispersal from other realms (Krasnov and Shenbrot, Reference Krasnov and Shenbrot2002; Gibert et al., Reference Gibert, Shenbrot, Stanko, Khokhlova and Krasnov2021). Here, we applied the concept of evoregions and evolutionary transition zones to 206 species of Palearctic fleas and 265 species of their hosts for identification of evoregions and evolutionary transition zones for parasites (fleas) and their hosts (small mammals). We aimed to understand whether these evoregions and transition zones for flea and hosts match each other. Strong concordance between flea and host evoregions would indicate an important role played by cophylogeny in the common flea-host histories, whereas weak (if present at all) concordance would suggest that the predominance of host-switchings either occurred locally or followed dispersal or both.
Materials and methods
Geographic distributions
We selected 206 flea species for which data on geographic ranges were available (Krasnov et al., Reference Krasnov, Poulin, Shenbrot, Mouillot and Khokhlova2005, Reference Krasnov, Shenbrot, van der Mescht, Warburton and Khokhlova2018; Shenbrot et al., Reference Shenbrot, Krasnov and Lu2007; Maestri et al., Reference Maestri, Shenbrot and Krasnov2017). Then, we selected 266 host species (a) on which at least one of these fleas was recorded (information taken from various literature sources; see references in Krasnov et al., Reference Krasnov, Shenbrot and Khokhlova2022b) and (b) with geographic ranges mainly situated in the Palearctic. The lists of flea and host species considered in this study are presented in Appendices 1 and 2, Supplementary Materials.
Geographic ranges for all species with 5 or more occurrence records were used as polygons. For all flea species and for rodents of the superfamily Dipodoidea, family Cricetidae, subfamily Gerbillinae, and several murine species (Apodemus agrarius, A. flavicollis, A. mystacinus, A. ponticus, Nesokia indica, and Rattus pyctoris), polygons of geographic ranges were obtained from species distribution modelling (SDM), whereas polygons of the geographic ranges of other small mammal hosts were obtained from the International Union of Nature Conservation (IUCN, 2022), corrected following Burgin et al. (Reference Burgin, Wilson, Mittemeier, Rylands, Lacher and Sechrest2020). Further details on the calculation of geographic distributions are presented in Appendix 3, Supplementary Material.
Phylogenies
For fleas, we used the most comprehensive molecular phylogenetic tree (Zhu et al., Reference Zhu, Hastriter, Whiting and Dittmar2015) as a backbone. This tree contained data on all families, subfamilies, and genera (but not all species) of fleas used in this study. The positions of the species absent from the Zhu et al. (Reference Zhu, Hastriter, Whiting and Dittmar2015) tree were derived from their taxonomic positions based on morphological traits and dichotomous identification keys (Hadfield et al., Reference Hadfield, Krasnov, Poulin and Nakagawa2014; see references in Krasnov et al., Reference Krasnov, Shenbrot and Khokhlova2022a). We assigned all branch lengths to an equal length of 1 because no information on branch lengths was available; we then ultrametrized the tree using the function ‘force.ultrametric’ (with option method = ‘extend’) of the package ‘phytools’ (Revell, Reference Revell2012) implemented in the R statistical environment (R Core Team, 2023).
Host phylogeny was taken as a subset of 1000 random trees from the 10 000 species-level birth-death tip-dated completed trees for 5911 mammal species of Upham et al. (Reference Upham, Esselstyn and Jetz2019). Then, we constructed a consensus tree using the function ‘consensus.edge’. Upham et al.'s (Reference Upham, Esselstyn and Jetz2019) tree for mammals is, to a certain degree, outdated, and many recently distinguished species (some dipodids, gerbillines, and arvicolines) are absent from this tree. We manually added these species to the resultant tree, and their topological positions were taken from various sources (see details and list of sources in Appendix 4, Supplementary Material). We assigned branch lengths for these species to an equal length of 1. Then, we ultrametrized the tree as described previously and resolved polytomies using the function ‘fix.poly’ implemented in the R package ‘RRphylo’ (Castiglione et al., Reference Castiglione, Tesone, Piccolo, Melchionna, Mondanaro, Serio, Di Febbraro and Raia2018) and polytomous clades to non-zero length branches.
Data analyses
We calculated the presence/absence of each flea or host species over a grid map of the Palearctic with 2804 cells of 2° × 2° and assigned a species' presence in each cell if the species' range covered at least 12.5% of the cell (otherwise, too many cells were excluded from the analyses, especially for species with narrow geographic distributions). The empty cells were excluded from the analyses.
Identification of evoregions for either fleas or hosts is based on the phylogenetic turnover among grid cells, measured using the phylogenetic fuzzy-weighting method (Pillar and Duarte, Reference Pillar and Duarte2010) that takes into account both between-species phylogenetic distances and tree imbalance [i.e. the difference between the numbers of tips descending from internal nodes (Duarte et al., Reference Duarte, Debastiani, Freitas and Pillar2016)]. The procedure involves several steps (see details in the study by Maestri and Duarte, Reference Maestri and Duarte2020; Nakamura et al., Reference Nakamura, Rodrigues, Luza, Maestri, Debastiani and Duarte2023). In brief, the first step is constructing 2 matrices, namely a matrix of pairwise phylogenetic covariances between flea or host species standardized by marginal totals (matrix Q; see Maestri and Duarte, Reference Maestri and Duarte2020 for details) and a presence/absence grid cell × species matrix. Multiplication of these matrices produces a matrix of the phylogenetic composition of species assemblages in the grid cells (matrix P; Maestri and Duarte, Reference Maestri and Duarte2020) that reflects phylogenetic turnover as the difference in phylogenetic composition between cells. In other words, matrix Q illustrates the degrees of the phylogenetic membership of each species to every other species, whereas matrix P is a phylogenetic fuzzy matrix reflecting the phylogenetically weighted degree of each species' belonging (ranging from 0 to 1) to each cell assemblage. Then, matrix P is used as an input for computing the Principal Coordinates of Phylogenetic Structure (PCPS of Duarte, Reference Duarte2011), which represents the Principal Coordinate Analysis (PCoA) of this matrix, using square-rooted Bray–Curtis dissimilarities between cells to avoid negative eigenvalues (Borcard et al., Reference Borcard, Gillet and Legendre2018). Gradients of phylogenetic turnover across cells are reflected in the eigenvectors produced by PCPS (Duarte et al., Reference Duarte, Debastiani, Freitas and Pillar2016). Subsequently, PCPS eigenvectors (using the principal components with more than 5% of explained variance) are taken as input data for the Discriminant Analysis of Principal Components based on k-means non-hierarchical clustering (DAPC; Jombart et al., Reference Jombart, Devillard and Balloux2010) that is used to perform biogeographic regionalization. The above analyses are implemented in the function ‘calc_evoregion’ of the R package ‘Herodotools’ (Nakamura et al., Reference Nakamura, Rodrigues, Luza, Maestri, Debastiani and Duarte2023). The optimal number of clusters obtained from DAPC is automatically calculated by the ‘calc_evoregion’ function based on the ‘elbow’ method from the R package ‘phyloregion’ (Daru et al., Reference Daru, Karunarathne and Schliep2020).
After identifying evoregions, we calculated the degree of affiliation of each species assemblage in the grid cell with an evoregion that this assemblage was classified. The degree of affiliation (= degree of membership) with an evoregion varies between the assemblages of the grid cells within the evoregion because of varying phylogenetic similarity between these assemblages. This degree for a given cell was calculated as the mean phylogenetic dissimilarity between this cell and all other cells belonging to an evoregion, using the function ‘calc_affiliation_evoreg’ of the ‘Herodotools’ package. High affiliation values indicate assemblages that are highly similar to other assemblages of the same evoregion (i.e. assemblages characterized by low phylogenetic turnover), whereas low affiliation values indicate assemblages dissimilar from other assemblages of the same evoregion (i.e. assemblages characterized by high phylogenetic turnover). In other words, the latter assemblages constitute evolutionary transition zones (Maestri and Duarte, Reference Maestri and Duarte2020).
To understand whether a particular phylogenetic lineage mainly occurs in a particular evoregion, we calculated species affiliations with evoregions using the function ‘calc_spp_association_evoreg’ of ‘Herodotools’. Because many species occurred in more than one region, we applied the approach of Maestri and Duarte (Reference Maestri and Duarte2020) (albeit less conservatively) and considered a species to belong to an evoregion if 60% of the grid cells in which this species was recorded (70% in Maestri and Duarte, Reference Maestri and Duarte2020) were classified to this evoregion. Species that could not be affiliated with a single evoregion were denoted as widespread species (Maestri and Duarte, Reference Maestri and Duarte2020). Then, we estimated ancestral evoregions simulating stochastic character (evoregion-specific affiliation or being widespread) maps on a phylogenetic tree of either fleas or hosts using the R package ‘phytools’.
The congruence between evoregions identified for fleas and evoregions identified for their hosts was tested for the overall similarity of their spatial structure. This was done using the V-measure method developed by Nowosad and Stepinski (Reference Nowosad and Stepinski2018) to calculate the spatial association between regionalizations and derived from a measure used in computer science for comparing different clusterings of the same domain. The V-measure results from comparing 2 regionalization maps and ranges from 0 to 1 (from no congruence whatsoever to perfect congruence). This measure represents a harmonic mean of 2 metrics, namely homogeneity and completeness, that also range from 0 to 1 (see details in the study by Nowosad and Stepinski, Reference Nowosad and Stepinski2018). Homogeneity is the average homogeneity of the host evoregions with respect to the flea evoregions, whereas completeness is the average homogeneity of the flea evoregions with respect to the host evoregions. The V-measure, homogeneity, and completeness were calculated using the function ‘vmeasure_calc’ of the R package ‘sabre’ (Nowosad and Stepinski, Reference Nowosad and Stepinski2018).
Results
For fleas, we identified 5 evoregions (Fig. 1). Only one of these evoregions (evoregion E) was spatially continuous. It included the highly arid zones of North Africa and the Arabian Peninsula and was characterized by a single endemic flea (Parapulex chephrenis) (Fig. 2), whereas the other fleas inhabiting this evoregion also occurred in other evoregions. Evoregion A covered mainly temperate zones, including (a) almost all of Europe except the southern part of the Iberian Peninsula, the Caucasus, and the south of European Russia; and (b) Japan (Fig. 1). The affiliations of Japan's flea assemblages with this evoregion were weak (Fig. S1; Appendix 5, Supplementary Material). Flea lineages that were diversified in evoregion A were represented mainly by genus Palaeopsylla and subgenus Euctenophthalmus (Fig. 2). Evoregion B comprised the mountain zones of the Caucasus, western China, and the Chukotka Peninsula, as well northern Kazakhstan (Fig. 1). Evoregion C included the hot deserts of coastal North Africa, the Middle East, part of Mongolia, and the cold deserts of extreme northern Siberia (Fig. 1), although the degree of the affiliations of the Mongolian and northern Siberian flea assemblages was low (Fig. S1; Appendix 5, Supplementary Material). Flea assemblages in both evoregions B and C belonged to a variety of lineages from different families (Fig. 2). Evoregion D included most of Siberia and central Kazakhstan (Fig. 1). Lineages that predominantly diversified in evoregion D were the hystrichopsyllid subfamilies Rhadinopsyllinae and Neopsyllinae, leptopsyllid genera Ophthalmopsylla and Paradoxopsyllus (albeit only basal taxa of the latter), and fleas of the Coptopsyllidae family (Fig. 2). Reconstruction of ancestral states indicated that the majority of lineages at the deep phylogenetic level diversified in multiple regions, except some hystrichopsyllids such as Stenoponia, Rhadinopsylla, and Catallagia, as well as Coptopsyllidae, which likely diversified in the evoregion to which they belong (Fig. 2). Transition zones (i.e. zones of high phylogenetic turnover) between evoregions identified for fleas were represented by a boundary between (a) the highly arid and the coastal zones of North Africa (i.e. between evoregions B and E), (b) the arctic zones and the temperate zones of Europe and Siberia (i.e. between evoregion C and evoregions A and B, and (c) the Balkan Mountains and the Aegean Sea (i.e. between evoregions D and A to the west and B to the east; Fig. 3).
Evoregions (evolutionarily important biogeographic regions) for 206 species of Palearctic fleas. Different evoregions are denoted by different colours.
Phylogenetic tree of fleas with colours representing the predominant evoregions (at least 60% of a species' geographic range belongs to a given region). Widespread species are those in which 60% of their geographic range could not be attributed to a single evoregion.
Evolutionary transition zones between flea evoregions.
For hosts, 5 evoregions were also identified (Fig. 4). Three of these evoregions (B, C and E) were spatially continuous. Evoregion A mostly covered South Asia including the southern Far East and was associated with the diversification of Rattus/Niviventer rats and Mus mice, as well as Eothenomys, Craseomys, and Alexandromys arvicolines (Fig. 5). Evoregion B comprised arid and semi-arid Central Asia (including Anatolia, Turan, Iran, Kazakhstan, Mongolia, and Northwestern China). Lineages that diversified in this evoregion were mainly represented by true and mouse-like hamsters (Cricetinae and Calomyscidae, respectively), Alticola, Ellobius, and Neodon voles, the majority of jerboas (Dipodidae; except Jaculus), pikas (Ochotonidae), and gerbils of the Rhombomyina subtribe (Fig. 5). Evoregion C included North Africa and the Arabian Peninsula and was characterized as the diversification centre for Gerbillus gerbils, gundis (Ctenodactylidae), and Jaculus jerboas (Fig. 5). Evoregion D comprised Europe, southern and northernmost Siberia, the northern Far East, and part of eastern Siberia (except its easternmost portion). Lineages that diversified in evoregion D were lemmings (Lemmini), Chionomys and Microtus voles, dormice (Gliridae), some ground squirrels (Spermophilus), and water shrews (Neomys; Fig. 5). Finally, evoregion E comprised central Siberia and its easternmost part. It was difficult to discern the clades that radiated there; although this seemed to be the case for Sorex shrews, many of these shrews were identified as widespread (Fig. 5). Reconstruction of the ancestral states suggested that radiation of the clades classified to evoregions B, C and D took place in these evoregions (respectively), whereas species belonging to the remaining evoregions radiated elsewhere, and their communities in evoregions A and E were assembled via dispersal. This is partly supported by the low affiliation values of many assemblages of evoregion E (Fig. S2; Appendix 5, Supplementary Material), suggesting high phylogenetic heterogeneity within this evoregion. Transition zones between evoregions identified for small mammals are presented in Fig. 6. The transition zone between evoregions B and C covers the deserts of the Arabian Peninsula and the Zagros Mountains. The transition zone between evoregions A and B is represented by the Himalayas, whereas the transition zone between evoregions B and D occurs at the Atlas Mountains and the Strait of Gibraltar. In contrast, no transition zone can be envisaged between evoregions D and E.
Evoregions (evolutionarily important biogeographic regions) for 265 species of Palearctic small mammals harbouring at least one of 206 flea species for which evoregion regionalization was carried out. Different evoregions are denoted by different colours.
Phylogenetic tree of small mammals harbouring at least one of 206 flea species for which evoregion regionalization was carried out, with colours representing the predominant evoregions (at least 60% of a species' geographic range belongs to a given region). Widespread species are those in which 60% of their geographic range could not be attributed to a single evoregion.
Evolutionary transition zones between small mammal evoregions.
The V-measure (Nowosad and Stepinski, Reference Nowosad and Stepinski2018) was 0.33, with homogeneity being 0.33 and completeness being 0.34. This indicated relatively low similarity between regionalizations for fleas and hosts.
Discussion
We found that evoregions for fleas and their small mammalian hosts, in general, did not match. In other words, evolutionarily important bioregions for fleas and hosts did not coincide. Obviously, host diversification did not depend on flea diversification. However, the opposite is definitely true (Traub, Reference Traub, Traub and Starcke1980; Whiting et al., Reference Whiting, Whiting, Hastriter and Dittmar2008). This, however, does not mean that flea diversification strictly follows that of hosts because, otherwise, the phylogenies of fleas and hosts, across all of the trees, would be highly congruent, whereas this is generally not the case (Krasnov and Shenbrot, Reference Krasnov and Shenbrot2002; Lu and Wu, Reference Lu and Wu2005; Llaberia-Robledillo et al., Reference Llaberia-Robledillo, Lucas-Lledó, Pérez-Escobar, Krasnov and Balbuena2023). Moreover, Llaberia-Robledillo et al. (Reference Llaberia-Robledillo, Lucas-Lledó, Pérez-Escobar, Krasnov and Balbuena2023) analysed cophylogenetic congruence and incongruence between 130 small mammalian species and 202 flea species and found that the degree of cophylogenetic congruence was concentrated in some (e.g. the majority of rodents), but not other (e.g. pikas and shrews), host clades, whereas the degree of phylogenetic congruence and incongruence of fleas greatly varied within clades (i.e. families, subfamilies, and genera). This indicated that codiversification with hosts was characteristic for some, but not other, flea lineages even if they belong to the same higher taxon (e.g. family). As a result, some flea lineages exploit an extremely narrow range of closely related hosts, whereas other confamilial flea lineages exploit hosts of various phylogenetic positions (Traub, Reference Traub and Kim1985). For example, the ceratophyllid genus Oropsylla exclusively parasitizes marmotine rodents, whereas hosts of the ceratophyllid Thrassis include a variety of rodent families: Sciuridae, Heteromyidae and Muridae. Exploitation of hosts belonging to different phylogenetic lineages by fleas belonging to the same lineage suggests that host switching was one of the most frequent events in flea evolutionary history. In the framework of this study, frequent host switching seemed to be the main reason explaining weak (if present at all) matching between flea and host evoregions.
The main mechanism of host switching is the so-called ecological fitting (Brooks et al., Reference Brooks, León-Règagnon, McLennan and Zelmer2006; Hoberg and Brooks, Reference Hoberg and Brooks2008; Agosta et al., Reference Agosta, Janz and Brooks2010; Araujo et al., Reference Araujo, Braga, Brooks, Agosta, Hoberg, von Hartenthal and Boeger2015). Ecological fitting (Janzen, Reference Janzen1985) presents a scenario in which a species interacts with its environment in a way resembling a shared evolutionary history, whereas in reality, the traits allowing the species to exploit this environment evolved in a different place and as a response to different conditions. For parasites, ecological fitting may occur if its main requirement is the resource sensu stricto rather than its natural carrier (a given host species). Consequently, if this resource is shared by many host species, a parasite tracking the resource may (1) invade new areas where the resource is present, but where its original host is not; (2) switch to a co-occurring but novel host species if the original host is extinct or develops novel defence tools (Nuismer and Thompson, Reference Nuismer and Thompson2006); and (3) switch to an invading host (Brooks et al., Reference Brooks, León-Règagnon, McLennan and Zelmer2006). Obviously, the main resources required by fleas from their hosts are food (i.e. hosts' blood) and a place for pre-imaginal stages to develop (burrow/nest/den). Although these resources are shared by many small mammals, fleas' patterns of acquiring these resources often differ, even between closely related hosts. For example, a flea's acquisition of a host's blood depends not only on the host's morphological traits (e.g. skin thickness), which are similar in close relatives (e.g. Sokolov, Reference Sokolov1982), but also on the host's defensive abilities (anti-parasitic grooming or immune response), which may differ between close relatives (Khokhlova et al., Reference Khokhlova, Spinu, Krasnov and Degen2004; Goüy de Bellocq et al., Reference Goüy de Bellocq, Krasnov, Khokhlova and Pinshow2006). A flea's utilization of a host's shelter as a place for pre-imaginal development depends on the shelter's microclimate (because fleas are sensitive to temperature and relative humidity; Marshall, Reference Marshall1981; Krasnov, Reference Krasnov2008), which is determined by shelter architecture, as well as the soil and vegetation structure of the host's habitat (Krasnov et al., Reference Krasnov, Shenbrot, Medvedev, Khokhlova and Vatschenok1998). Therefore, the probability of host switching is driven not only by the co-occurrence of a potential host and a given parasite but also by the suitability of this host's traits for this parasite, which, in turn, depends on the traits of the parasite itself. Indeed, Krasnov et al. (Reference Krasnov, Shenbrot, Khokhlova and Degen2016) showed that fleas possessing a certain combination of traits, independent of their phylogenetic affinities, utilized hosts that were also characterized by certain trait combinations, although fleas exploiting the same host species were more phylogenetically related than expected by chance. In other words, the species compositions of fleas' host spectra, as well as the species composition of hosts' flea assemblages, are driven by complex interactions between the phylogenies and the traits of both fleas and hosts.
Affiliations of flea species with evoregions and our attempts to reconstruct ancestral geographic distributions (Fig. 2) indicated mostly ex-situ flea diversification, even at the family level. Assuming the monophyly of a majority of flea families (except Hystrichopsyllidae; Zhu et al., Reference Zhu, Hastriter, Whiting and Dittmar2015), this suggested that (a) initial diversification of these families took place out of the Palearctic and that (b) dispersal was the most likely assembly mechanism of flea communities in evoregions. In fact, the Gondwanan (Australia and South America, terrestrially connected via Antarctica until the upper Eocene) origin of the majority of flea families is now commonly accepted (Medvedev, Reference Medvedev1996, Reference Medvedev2005; Zhu et al., Reference Zhu, Hastriter, Whiting and Dittmar2015). However, many flea clades rapidly radiated only after their ancestors migrated to Eurasia and Africa (Zhu et al., Reference Zhu, Hastriter, Whiting and Dittmar2015). On the one hand, this suggests that if evoregions for fleas were identified at larger (e.g. global) scales, the geographic distribution of evoregions would differ from that of the smaller scale (e.g. the Palearctic as in this study). Indeed, the distribution of the evoregions identified for the majority of muroid rodents by Maestri and Duarte (Reference Maestri and Duarte2020) at a global scale substantially differed from those of the small mammalian hosts in our study (although muroids constitute the major proportion of them) and, to some extent, matched biogeographic realms. Moreover, the reconstruction of ancestral evoregions for muroids demonstrated that ex situ diversification was characteristic of many lineages at the deep phylogenetic level. On the other hand, consideration of evoregions at a smaller scale allows a better and more detailed understanding of geographic patterns of diversification for a given clade at the shallow phylogenetic level. Identification of evoregions for fleas at a global scale warrants further investigation. This will require a gargantuan effort to model the geographic ranges of the majority of flea species (more than 2500; Medvedev, Reference Medvedev2005) and to construct their molecular phylogeny. However, the data needed for this endeavour are largely unavailable.
No indication of in-situ diversification supports the results of Gibert et al. (Reference Gibert, Shenbrot, Stanko, Khokhlova and Krasnov2021), who applied a permutation-based algorithm to infer the relative roles of niche-based vs dispersal-based mechanisms in the assembly of regional flea communities in 4 biogeographic realms. They found that these communities' assembly was, to a great extent, governed by dispersal processes and, to a much lesser extent, by niche-based processes. The role of dispersal processes in the assembly of flea and host communities in evoregions is supported, albeit indirectly, by the occurrence of evolutionary transition zones. These zones were most likely determined by natural boundaries, such as deserts or mountain ranges, that could act as barriers to dispersal. This seems to be the case for hosts but not necessarily for fleas, because some of the transition zones detected for them did not correspond to landscape boundaries. For example, an occurrence of a transition zone between evoregions E and C in North Africa probably represented a transition between the climates of the Mediterranean coast (more humid) and the Sahara Desert (highly arid), so that flea species with pre-imagoes highly sensitive to desiccation inhabited evoregion C but not E. Similarly, the transition zone of extreme northern Eurasia could be associated with the sharp difference in soil temperature between evoregions C and A/D. The phylogenetic conservatism of either the degree of sensitivity to environmental factors or the preferable temperature/humidity regime in fleas has never been studied, but similar seasonality of reproduction and activity in many closely related fleas (Darskaya, Reference Darskaya1970; Krasnov, Reference Krasnov2008) suggests that this may be the case.
In conclusion, a comparison of biogeographic regionalization, coupled with phylogenetic information, between parasites and their hosts gives insight into the shared patterns and processes of 2 different histories. In particular, our study demonstrated that the application of the evoregion approach allowed a better understanding of the contribution of dispersal to cophylogenetic patterns.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182023000884.
Data availability statement
Raw data can be obtained from the corresponding author upon request.
Acknowledgements
We thank Renan Maestri and Gabriel Nakamura for their help with the R code of the ‘Herodotools’ package.
Author's contributions
BRK conceived and designed the study. Both authors collected the data and performed statistical analyses. BRK wrote the first draft of the article. Both authors finalized the article.
Financial support
This study did not receive any specific financial support.
Competing interests
The authors declare there are no competing interest.
Ethical standards
This study is based on published data, and therefore, ethical standards are not applicable.
Open access
