Introduction
One of the central aims of ecology is to understand patterns of biological diversity and processes behind these patterns. It is universally accepted that biological diversity is divided across scales into α diversity (diversity within a site), β-diversity (diversity among localities or timepoints) and γ diversity (total diversity within a landscape). The concepts of α and γ diversity can be easily grasped, whereas the concept of β-diversity is more complicated. One of the reasons for this is that spatial and/or temporal β-diversity of a biological community represents a combination of different components. For example, Pavoine and Dolédec (Reference Pavoine and Dolédec2005) and Pavoine et al. (Reference Pavoine, Marcon and Ricotta2016) argued that real biological units (e.g. communities) exist in a complex nested hierarchy of spaces such as patches within habitats, habitats within landscapes, landscapes within geographic regions, so that spatial β-diversity can be partitioned into (a) dissimilarity of communities among units of a finer scale within units of a broader scale and (b) dissimilarity of communities among units of a broader scale. Looking at β-diversity from another angle, Baselga (Reference Baselga2010) proposed to deconstruct β-diversity into two dissimilarity components, namely dissimilarity due to ‘pure’ spatial turnover and dissimilarity due to nestedness. Recently, Legendre and De Cáceres (Reference Legendre and De Cáceres2013) developed a method allowing disentangling between the effects of individual species and the effects of species assemblages ( = sites) on total β-diversity. In this approach, total β-diversity is estimated via total variance of the community data matrix and then it is partitioned into (a) the degree of relative importance of individual species for β-diversity across area [representing thus species contributions to β-diversity (SCBD)] and (b) indicators of the compositional uniqueness of the assemblages as compared with other sites [that is local contributions to β-diversity (LCBD)]. This partitioning of β-diversity appears to be helpful in understanding the reasons behind differential roles of species and assemblages in total β-diversity. On the one hand, the contribution of a species to β-diversity (SCBD) may be affected by various species-specific traits such as abundance (Vilmi et al., Reference Vilmi, Karjalainen and Heino2017; da Silva et al., Reference da Silva, Hernández and Heino2018), niche breadth (Heino and Grönroos, Reference Heino and Grönroos2017) or life history traits affecting the interactions of this species with other species or the environment. On the other hand, the contribution of an assemblage to β-diversity (LCBD) can be related to environmental factors (e.g. Tonkin et al., Reference Tonkin, Heino, Sundermann, Haase and Jähnig2016) as well as to the structure of the target biological community (e.g. da Silva et al., Reference da Silva, Hernández and Heino2018).
Development of the Legendre and De Cáceres’ (Reference Legendre and De Cáceres2013) approach for decomposing β-diversity triggered a burst of recent studies that related SCBD and/or LCBD to species-specific and site-specific characteristics, respectively, in a variety of plant and animal taxa (Tonkin et al., Reference Tonkin, Heino, Sundermann, Haase and Jähnig2016; Heino and Grönroos, Reference Heino and Grönroos2017; Kong et al., Reference Kong, Chevalier, Laffaille and Lek2017; Vilmi et al., Reference Vilmi, Karjalainen and Heino2017; da Silva et al., Reference da Silva, Hernández and Heino2018; Landeiro et al., Reference Landeiro, Franz, Heino, Siqueira and Bini2018; Tolonen et al., Reference Tolonen, Leinonen, Erkinaro and Heino2018). However, only two studies applied partitioning of β-diversity into SCBD and LCBD components to parasitic species. Biguezoton et al. (Reference Biguezoton, Adehan, Adakal, Zoungrana, Farougou and Chevillon2016) analysed β-diversity of ixodid ticks, whereas Poisot et al. (Reference Poisot, Guéeveneux-Julien, Fortin, Gravel and Legendre2017) studied β-diversity of fleas. The former study considered both SCBD and LCBD, but did not relate either SCBD to species traits or LCBD to site characteristics. The latter study considered LCBD only (both for parasites and hosts as well as their interactions) and analysed their relationships with climatic variables. Both studies thus missed an important component of parasite biology, namely their utmost dependence on host species. Indeed, contrary to free-living species, parasites have a dual environment including, on the one hand, traditional environmental factors such as climate and vegetation, and, on the other hand, an assemblage of hosts which has supreme importance for parasites because parasites cannot exist without their hosts. Consequently, SCBD and LCBD of parasites cannot be completely understood without testing for the effects of hosts on these metrics. In particular, SCBD of parasites are likely affected by their host specificity (i.e. the size and/or composition of host spectra), whereas LCBD of parasite assemblages are likely affected by host species composition of a site or region.
Here we asked what are the factors affecting SCBD and LCBD of fleas parasitic on small mammals in 45 regions of the Palearctic. We tested (a) whether variation in ecological, morphological, life history and geographic traits can predict SCBD of fleas and (b) if variation in flea and host community metrics, off-host environmental factors and host species composition can predict LCBD of flea assemblages. In addition, and following Tonkin et al. (Reference Tonkin, Heino, Sundermann, Haase and Jähnig2016) and da Silva et al. (Reference da Silva, Hernández and Heino2018), we tested whether spatial variables may describe the geographic distribution of regions with various degrees of uniqueness of flea species composition.
Materials and methods
Data on flea and host species composition
We used data on species composition of fleas and their small mammalian hosts (Soricomorpha, Rodentia and ochotonid Lagomorpha) in 45 regions of the northern and temperate Palearctic compiled from published surveys and used in our earlier studies [see map of the regions and references in Krasnov et al. (Reference Krasnov, Stanko, Khokhlova, Shenbrot, Morand, Korallo-Vinarskaya and Vinarski2011, Reference Krasnov, Shenbrot, Khokhlova, Stanko, Morand and Mouillot2015)]. In this study, we used only data from the surveys in which at least 500 host individuals were examined. We omitted from our data host species such as squirrels, flying squirrels, chipmunks and hedgehogs and their specific fleas (e.g. Tarsopsylla octodecimdentata and Archeopsylla erinacei) because these hosts require special trapping methods and they were either undersampled or not sampled at all. In total, the database included 202 flea species collected from 136 small mammal species.
Trait variables
Flea species were characterized by 10 traits. Three of these traits represented ecological variation (characteristic abundance, mean size of a host spectrum, phylogenetic distinctness of a global host spectrum). These quantitative traits have been shown to be true attributes of a flea species and vary significantly less among populations of the same species than among species (Krasnov et al., Reference Krasnov, Mouillot, Shenbrot, Khokhlova and Poulin2004a, Reference Krasnov, Shenbrot, Khokhlova and Poulin2006; Mouillot et al., Reference Mouillot, Krasnov, Shenbrot, Gaston and Poulin2006; Krasnov and Poulin, Reference Krasnov, Poulin, Morand and Krasnov2010). The size and phylogenetic distinctness of a host spectrum essentially represent niche breadth of a flea (Garnick, Reference Garnick1992, Poulin et al., Reference Poulin, Krasnov and Mouillot2011). Two traits represented bionomics ( = life history) of fleas: microhabitat preference and seasonal peak of reproduction. Microhabitat preference is determined by the proportion of time fleas spend on either the body of a host or in its burrow/nest. This proportion is a species-specific trait that distinguishes so-called ‘body’ fleas from ‘nest’ fleas [see details in Krasnov (Reference Krasnov2008)]. Seasonal peak of reproduction may occur during either the warm or cold period (‘summer’ or ‘winter’, respectively) or else a flea is able to reproduce all year without clear seasonality. Three morphological traits included armament (i.e. occurrence and/or the number of sclerotinized combs, an ordinal variable; see Krasnov et al., Reference Krasnov, Shenbrot, Khokhlova and Degen2016) as well as body size and the degree of sexual size dimorphism (both quantitative variables). The remaining two traits were biogeographic: size and latitude of the centre of geographic range. See Supplementary Material Appendix 1 for detailed explanations on the selection and measurements of flea traits as well as calculations of characteristic abundance, host specificity indices and sexual size dimorphism.
Following Heino and Grönroos (Reference Heino and Grönroos2017) and da Silva et al. (Reference da Silva, Hernández and Heino2018), we substituted the original values of bionomic and morphological traits with the ‘morphobionomic trait’ vectors. The latter were the first two principal coordinate axes extracted from the principal coordinate analysis of distance trait matrix constructed using Gower distance coefficient with the function ‘gowdis’ using package ‘FD’ (Laliberté and Legendre, Reference Laliberté and Legendre2010) of the R statistical environment (R Core Team, 2018). Principal coordinate analysis was carried out using function ‘cmdscale’ of the R base package ‘stats’. The two morphobionomic trait vectors (MBV1 and MBV2) explained 91.05% of variation and were continuous variables representing between-species differences in their morphological and life history traits. To guarantee that morphobionomic trait vectors adequately reflected original between-species distances in traits, we ran Mantel's correlation test between the matrix of Euclidean distances based on coordinates along MBV1 and MBV2 and the original matrix of Gower distances and found strong correlation between the two matrices (Mantel r = 0.84, P < 0.001). The first vector (MBV1) represented a gradient of body size of flea species from large to small, whereas the second vector (MBV2) showed a decrease in the degree of sexual size dimorphism (see details in Supplementary Material Appendix 1, Table S1). The relationships between ordinal (i.e. armament) and nominal (i.e. ‘nest’ vs ‘body’ vs ‘both’ preference) morphobionomic traits and the respective trait vectors are presented in Supplementary Material Appendix 1, Table S2 and Fig. S1. In brief, the negative area of MBV1 was mainly associated with species either preferring to stay on a host's body (i.e. ‘body’ fleas) or without clear microhabitat preference, demonstrating either summer or non-seasonal reproduction and having only a pronotal comb, whereas coordinates of ‘nest’ fleas with winter reproduction and either full (a pronotal and a genal) or no comb armament along this vector were mainly positive. ‘Nest’ species or those without clear microhabitat preferences having summer reproduction and either both combs or no combs at all were distributed mainly in the positive area of MBV2, whereas the opposite was true for ‘body’ fleas that reproduce in winter or all year round and have a pronotal comb only. There was no clear relationship between the pattern of flea armament and coordinates along MBV2 (Supplementary Material Appendix 1, Fig. S1).
Environmental (off-host and host-associated) variables
We characterized sampled regions by their off-host and host-associated environment. We calculated the sampling area in a sampled region based on data from an original source. The latitudinal and longitudinal positions of the region centres were determined using ArcGIS 10.6. Regional environmental variables included topography (mean altitude), amount of green vegetation (normalized difference vegetation indices), separately for autumn, winter, spring and summer], air temperature (mean, maximum and minimum as well as annual and monthly ranges) and precipitation (for autumn, winter, spring and summer). Environmental data were averaged over a region across 30 arc-second grids. Elevation data were obtained using ArcGIS 10.6. NDVI data were taken from the VEGETATION Programme (http://free.vgt.vito.be), whereas temperature and precipitation variables were extracted from WORLDCLIM (BIOCLIM) 2.0 package (Fick and Hijmans, Reference Fick and Hijmans2017). We ran principal component analyses of environmental variables and then substituted their original values with the scores of the three first principal components (PC1, PC2 and PC3). These three components explained 79.70% of environmental variation. PC1 represented mainly an increase in annual range of air temperatures and a decrease in the amount of green vegetation, PC2 correlated negatively with the mean, maximal and minimal air temperatures, whereas PC3 reflected an increase in spring and winter precipitation and mean altitude of a region (see details in Supplementary Material Appendix 2, Table S2).
To characterize a region from the perspective of host species composition, we calculated the numbers of host species belonging to one of 19 phylogenetic/ecological clades/groups (see Supplementary Material Appendix 3, Table S3) in each region [see justification in Krasnov et al. (Reference Krasnov, Shenbrot, Khokhlova, Stanko, Morand and Mouillot2015)]. This allowed us to take into account both numerical and phylogenetic facets of a regional host community composition. These host-associated variables were considered ordinal. To transform these variables into continuous predictors for further analysis, we applied principal coordinate analysis on a distance matrix based on Manhattan distance and then extracted the first three ‘host composition’ vectors (HCV1, HCV2 and HCV3). This was done using the R packages ‘vegan’ (Oksanen et al., Reference Oksanen, Blanchet, Friendly, Kindt, Legendre, McGlinn, Minchin, O'Hara, Simpson, Solymos, Stevens, Szoecs and Wagner2018) and ‘labdsv’ (Roberts, Reference Roberts2016). Similarly to morphobionomic trait vectors (see above), these host composition vectors were continuous variables describing differences between regions in host composition and explained 82.0% of among-region variation (see details in Supplementary Material Appendix 3, Table S4). HCV1 was responsible for more than half of the explained variation and reflected a compositional gradient from bank voles, moles and soricine shrews (negative area) to gerbils, jerboas and ground squirrels (positive area); HCV2 was associated with a decrease in the number of species of Microtus voles; whereas HCV3 correlated negatively with the number of Ochotona and Myodes species and positively with species richness of crocidurines. The matrix of Euclidean distances based on coordinates on HCV1–HCV3 and the original matrix of Manhattan distances were strongly correlated (Mantel r = 0.97, P < 0.001) demonstrating that the three vectors of host composition adequately reflected between-region differences in host composition.
Structure of flea and host communities
We used flea and host species richness as well as host compositional uniqueness as metrics of structure of their communities. Estimates of flea and host species richness across regions can be affected by unequal among-region sampling effort (i.e. number of host individuals examined; Morand and Poulin, Reference Morand and Poulin1998) and unequal sampled area (Rosenzweig, Reference Rosenzweig1995). Flea species richness was significantly, albeit weakly, affected by sampling effort but was not affected by the area of a region (r 2 = 0.10, F = 4.83, P = 0.03 and r 2 = 0.03, F = 0.52, P = 0.23, respectively), whereas host species richness was affected by the area (r 2 = 0.23, F = 13.03, P < 0.001). To control for the confounding effects of sampling effort and the area of a region, we regressed (a) the number of flea species recorded against the number of host individuals examined and (b) the number of host species against the sampled area in the log–log space. Then, we substituted the original values of flea and host species richness by their residual deviations from these regressions. Compositional uniqueness of host communities was assessed via their local contribution to total host β-diversity ( = LCBD for hosts; see below for details on calculation of this metric).
Spatial variables
Following da Silva et al. (Reference da Silva, Hernández and Heino2018), we applied the approach of Moran Eigenvectors Maps (MEM) [also known as the principal coordinates of neighbouring matrices (PCNMs)] (Borcard and Legendre, Reference Borcard and Legendre2002; Dray et al., Reference Dray, Legendre and Peres-Neto2006) in which spatial variables are calculated from geographic coordinates of each region. We built a truncated Euclidean distance matrix from coordinates and then extracted eigenvectors with positive eigenvalues that were used as predictors in further analysis (Borcard and Legendre, Reference Borcard and Legendre2002). In general, MEM allows evaluation of spatial structures over a range of scales with the first to the last eigenvectors representing ever-decreasing spatial scales from broad to fine (Borcard and Legendre, Reference Borcard and Legendre2002), although this is mainly true for a spatially regular sampling design (Dray et al., Reference Dray, Legendre and Peres-Neto2006; Borcard et al., Reference Borcard, Gillet and Legendre2018). We extracted 25 spatial variables ( = positive eigenvectors further referred to as PCNMs) using the R package ‘vegan’. Then, we tested for spatial autocorrelation in these variables using Moran's I test implemented in the R package ‘ape’ (Paradis et al., Reference Paradis, Claude and Strimmer2004). PCNMs that demonstrated significant spatial autocorrelation (13 of 25 PCNMs) were retained for the main analysis (see Borcard et al., Reference Borcard, Gillet and Legendre2011 for justification and explanations).
Data analyses
We calculated species and local ( = region) contributions to β-diversity of fleas (SCBD and LCBDf, respectively) using the approach of Legendre and De Cáceres (Reference Legendre and De Cáceres2013). We constructed a presence/absence region-by-species matrix which was subsequently Hellinger-transformed. This transformation takes the square root of the quotient of the value of each cell in a data matrix and its row sum and thus gives low weights to species with low counts and many zeros across regions. Then, we calculated SCBD values for each species and LCBDf values for each region using function ‘beta.div’ and Hellinger dissimilarity coefficient implemented in the R package ‘adespatial’ (Dray et al., Reference Dray, Bauman, Blanchet, Borcard, Clappe, Guenard, Jombart, Larocque, Legendre, Madi and Wagner2018). Local contributions to β-diversity of hosts ( = LCBDh) was calculated in a similar way using a Hellinger-transformed presence/absence region-by-host species matrix.
We modelled the responses of (a) SCBD to trait variables and (b) LCBDf to environmental (both off-host and host-associated) variables, spatial variables and variables describing flea and host community metrics. Because the response variables ranged from 0 to 1 (without attaining either 0 or 1 values), we applied β-regression with a logit link function. The β-regression approach is the most suitable for models with a response variable of this type because, in particular, it naturally incorporates heteroscedasticity and skewness inherent to these variables [Cribari-Neto and Zeileis, Reference Cribari-Neto and Zeileis2010; see also Heino and Grönroos (Reference Heino and Grönroos2017) and da Silva et al. (Reference da Silva, Hernández and Heino2018)]. The β-regressions were carried out using the R package ‘betareg’ (Cribari-Neto and Zeileis, Reference Cribari-Neto and Zeileis2010).
To better understand the effect of predictors of different nature on either SCBD or LCBDf, we ran three separate models for the former and four separate models for the latter. First, we analysed how SCBD of a flea species varies in response to (a) ecological flea traits (abundance, size of a host spectrum and its phylogenetic distinctness); (b) morphobionomic traits (seasonality, armament, body size, the degree of sexual size dimorphism represented by two morphobionomic trait vectors, see above); and (c) biogeographic traits (size and latitude of geographic range). Second, we tested for the dependence of LCBDf on (a) structure of flea and host communities (flea and host species richness and contribution of regional host communities to total host β-diversity); (b) off-host environment (represented by three principal components of environmental variables, see above); and (c) host-associated environment in terms of the number of host species belonging to the main phylogenetic/ecological clades (represented by three host composition vectors, see above). Finally, (d) we used PCNMs (see above) as explanatory variables to assess the spatial distribution of LCBDf across the Palearctic.
Results
More than one-third of the 202 flea species (71 species) demonstrated a higher than average contribution to total β-diversity (SCBD). These species belonged to different genera and families (see Supplementary Material Appendix 4). The β-regression analyses demonstrated that SCBD significantly increased with an increase in mean abundance of a flea species and a decrease in its phylogenetic host specificity (i.e. an increase in phylogenetic distinctness of a host spectrum) (Table 1, Fig. 1). Fleas with large and more northern geographic ranges contributed to total β-diversity significantly more than fleas with small and more southern geographic ranges (Table 1, Fig. 2). Furthermore, geographic traits of fleas explained more variation in their SCBD than ecological traits (compare pseudo-R 2 of the models, Table 1). In contrast to ecological and geographic traits, morphological and life history traits had no effect on SCBD (Table 1).
Fig. 1. Relationships between the contribution to total β-diversity of flea species (SCBD) and their characteristic abundances (A) and phylogenetic distinctness of host spectra (B).
Fig. 2. Relationships between the contribution to total β-diversity of flea species (SCBD) and the sizes (A) and latitudes of the centres (B) of their geographic ranges.
Table 1. Summary of β-regression analyses of the effect of ecological (model ET), morphobionomic (model MBT) and biogeographic traits (model GT) of flea species on their contributions to β-diversity (SCBD)
Morphobionomic traits were represented by two vectors produced by principal coordinate analysis (MBV1, MBV2; see text for explanations).
*–P < 0.05.
Almost half of 45 regional flea assemblages (21 assemblages) had above-average contributions to total flea β-diversity (see Supplementary Material Appendix 5, Fig. 2S). LCBDf did not depend on either flea or host species richness, but was significantly affected by compositional uniqueness of host assemblages (Table 2, Fig. 3). Variables associated with off-host environment had no effect whatsoever on LCBDf (Table 2), whereas the opposite was true for vectors describing composition of regional host assemblages (Table 2). For example, the contribution of flea assemblages to total β-diversity was higher in regions with a greater number of gerbil, jerboa, ground squirrel and marmot species than in regions with a high diversity of bank voles and Sorex shrews (Fig. 4). Among spatial variables, only PCNMs 2 and 4 were significantly associated with LCBDf (Table 2) with the former being positively and the latter being negatively related to LCBDf. Distribution of PCNMS and their values is presented in Fig. 5. Host compositional uniqueness and host composition explained almost equal proportions of variation in LCBDf, whereas the contribution of spatial variables was 1.5 times lower (Table 2).
Fig. 3. Relationship between the contribution of regional flea assemblages to total flea β-diversity (LCBDf) and the contribution of regional host assemblages to total host β-diversity (LCBDh).
Fig. 4. Relationship between the contribution of regional flea assemblages to total flea β-diversity (LCBDf) and a variable describing host composition of a region (host composition vector 1).
Fig. 5. Map of values of the PCNM variables selected to model LCBD of flea assemblages [(A) PCNM2, (B) PCNM4, see text for explanations]. Sizes of squares are proportional to positive (black) and negative (white) values.
Table 2. Summary of β-regression analyses of the effect of regional flea and host community metrics (model CM), off-host environmental variables (model OHE), variables describing host composition (model HCE) and spatial variables (model SP) on contributions of regional flea assemblages to β-diversity (LCBDf)
LCBDh is contributions of regional host assemblages to total host β diversity. Off-host environment was represented by three principal components (PC1, PC2, PC3); host composition was represented by three vectors produced by principal coordinate analysis (HCV1, HCV2, HCV3); spatial variables were represented by positive eigenvectors (PCNMs) calculated from distances between regions (see text for explanations). For model SP, only predictors with significant coefficients are shown.
*–P < 0.05.
Discussion
We found that SCBD was affected by ecological and biogeographic, but not morphological and life history traits of flea species. LCBD of flea assemblages was strongly related to LCBD of host assemblages and host species composition, but not to vegetation and climatic variables. In other words, dissimilarity of a given host assemblage from other host assemblages drives dissimilarity of its flea assemblage from other flea assemblages. In addition, LCBD of fleas demonstrated some spatial structure at a broad-to-moderate scale.
Positive relationships between species abundance and SCBD have also been found in stream insect (Heino and Grönroos, Reference Heino and Grönroos2017), fish (Kong et al., Reference Kong, Chevalier, Laffaille and Lek2017) and dung beetle (da Silva et al., Reference da Silva, Hernández and Heino2018) communities, although this was mainly true when SCBD was calculated using species abundance but not presence/absence data. In addition, in both Heino and Grönroos (Reference Heino and Grönroos2017) and da Silva et al. (Reference da Silva, Hernández and Heino2018), the quadratic term of the abundance variable was found to be significant, whereas no curvilinearity can be envisaged from the relationship between abundance and SCBD in our study (Fig. 1). As mentioned above, flea abundance has been shown to vary only within relatively narrow species-specific boundaries and thus can be considered as a true flea species trait (Krasnov et al., Reference Krasnov, Shenbrot, Khokhlova and Poulin2006). Our results thus demonstrate that species with inherently high abundance tended to contribute more to total β-diversity than species with low abundance. Heino and Grönroos (Reference Heino and Grönroos2017) suggested that the relationship between abundance and SCBD may result from a positive relationship between abundance and occupancy [see also Legendre and De Cáceres (Reference Legendre and De Cáceres2013)], so that species with high SCBD values are expected to show high local abundance and high occupancy as well as a large variation in abundance across locations. This is unlikely true for fleas because of low among-location variation in the abundance metric considered in this study (Krasnov et al., Reference Krasnov, Shenbrot, Khokhlova and Poulin2006). An alternative explanation for the effect of mean abundance on SCBD can be linked to high among-species variation in abundance (coefficient of variation is 0.97) and a small proportion of species with relatively high characteristic abundance (e.g. >1; 24 of 202 species). As a result, assemblages containing highly abundant species might be strongly dissimilar from the remaining assemblages.
The relationship between SCBD and niche breadth in free-living insects has been found to be either very weak (Heino and Grönroos, Reference Heino and Grönroos2017) or absent (da Silva et al., Reference da Silva, Hernández and Heino2018). On the contrary, this relationship appeared to be significant for fleas. This discrepancy might be a result of different approaches to the definition of niche breadth. In our case, we considered niche breadth as a diversity of host species exploited rather than defining it as a tolerance to environmental variables. Positive association between SCBD and flea niche breadth (inversely indicated by phylogenetic host specificity) can, at least in part, result from positive relationships between this trait and abundance. Indeed, host generalist fleas that tend to exploit phylogenetically/taxonomically unrelated species usually attain higher abundances (Krasnov et al., Reference Krasnov, Poulin, Shenbrot, Mouillot and Khokhlova2004b). Furthermore, our earlier studies demonstrated that phylogenetic/taxonomic distinctness of host spectra in fleas (i) is positively related to their geographic range (Krasnov et al., Reference Krasnov, Poulin, Shenbrot, Mouillot and Khokhlova2005) and (ii) increases towards higher latitudes (Krasnov et al., Reference Krasnov, Shenbrot, Mouillot, Khokhlova and Poulin2008). In addition, fleas with more northerly distributions have larger geographic ranges (Krasnov et al., Reference Krasnov, Shenbrot, Mouillot, Khokhlova and Poulin2008). These complex inter-relationships between characteristic abundance, phylogenetic host specificity and geographic distribution result thus in predictability of flea species contribution to total β-diversity. This is indirectly supported by the associations between SCBD and (a) geographic range position and (b) phylogenetic, but not numerical, host specificity. Our earlier study (Krasnov et al., Reference Krasnov, Shenbrot, Mouillot, Khokhlova and Poulin2008) demonstrated significant correlation between latitude of the geographic range and the former, but not the latter, measure of host specificity. In other words, highly abundant and broadly distributed host generalists with northern geographic ranges contributed more to total β-diversity than less abundant host-specific fleas with narrow and southern geographic ranges. The effect of geographic range position on SCBD seems to be further translated into spatial pattern of LCBD (see below). In addition, the distribution of data points on Fig. 1B is clearly triangular. This means that fleas exploiting closely related hosts invariantly demonstrated low SCBD, while SCBD of fleas recorded on unrelated hosts might be either high or low. This indicates an important role of the composition of a host spectrum.
LCBD of regional flea assemblages was found to be determined mainly by composition of regional host assemblages and their LCBD ( = LCBDh), but not species richness of either fleas or hosts and environmental variables. The higher LCBD of a host assemblage is, the higher LCND of the respective flea assemblage is as well. In other studies of free-living species (Tonkin et al., Reference Tonkin, Heino, Sundermann, Haase and Jähnig2016; Heino and Grönroos, Reference Heino and Grönroos2017; Vilmi et al., Reference Vilmi, Karjalainen and Heino2017; Landeiro et al., Reference Landeiro, Franz, Heino, Siqueira and Bini2018; Silva et al., Reference da Silva, Hernández and Heino2018; Tolonen et al., Reference Tolonen, Leinonen, Erkinaro and Heino2018) and parasites (Poisot et al., Reference Poisot, Guéeveneux-Julien, Fortin, Gravel and Legendre2017), the relationships between LCBD and environmental variables were found to be important. Nevertheless, the associations between LCBD and environmental predictors in some of these studies were weak (e.g. Heino and Grönroos, Reference Heino and Grönroos2017; da Silva et al., Reference da Silva, Hernández and Heino2018). The main reason behind these contradictory results can merely be the differences in life histories of taxa for which the effect of environment on LCBD was tested. Indeed, Landeiro et al. (Reference Landeiro, Franz, Heino, Siqueira and Bini2018) tested for the effects of environment on LCBD of 14 plant and animal taxa from the same region and found that LCBD of different taxa (or even different ecological groups within the same taxon) were, if at all, associated with different environmental variables. The most surprising discrepancy in the effect of environmental variation on LCBD is between our results and those of Poisot et al. (Reference Poisot, Guéeveneux-Julien, Fortin, Gravel and Legendre2017) that reported the effect of some climatic predictors on compositional uniqueness of flea assemblages. Although we and Poisot et al. (Reference Poisot, Guéeveneux-Julien, Fortin, Gravel and Legendre2017) used essentially the same data, the difference in methodological approach is the most likely reason behind this discrepancy between results. First, Poisot et al. (Reference Poisot, Guéeveneux-Julien, Fortin, Gravel and Legendre2017) related differences in LCBD among flea assemblages to differences in the environment using redundancy analysis, while we applied β-regressions for our analyses. Second and in contrast to Poisot et al. (Reference Poisot, Guéeveneux-Julien, Fortin, Gravel and Legendre2017), we did not use the data from surveys in which sampling effort was too low (see above). We also did not include in the input matrices those host species that required a special sampling technique, so that these hosts and their fleas were apparently undersampled. Including undersampled regions and species could introduce strong bias in the analysis. For example, it is highly likely that inclusion of hosts and fleas sampled in some regions but not in other regions where they undoubtedly occurred (e.g. squirrels and their fleas) would substantially and spuriously inflate LCBD values of the former. Finally, Poisot et al. (Reference Poisot, Guéeveneux-Julien, Fortin, Gravel and Legendre2017) measured environmental variables for each region around its centre, whereas we averaged these variables across the entire sampling area within each region. Environmental conditions in the centre of a region do not necessarily represent the entire range of conditions across larger areas.
The lack of the effect of vegetation and climatic variables on LCBD found in this study does not, however, indicate irrelevance of the off-host environment for diversity and/or composition of flea assemblages. In fact, the role of environmental conditions of various aspects of flea communities has been emphasized in our earlier studies (e.g. Krasnov et al., Reference Krasnov, Mouillot, Shenbrot, Khokhlova and Poulin2010, Reference Krasnov, Shenbrot, Khokhlova, Stanko, Morand and Mouillot2015), including the effect of environment on one of the components of flea β-diversity, namely species turnover along gradients (Maestri et al., Reference Maestri, Shenbrot and Krasnov2017). It is possible that the relative importance of the environmental effects varies among different facets of diversity. Moreover, the effect of environment on LCBD may be important or not in dependence of a taxon and/or a region of interest (e.g. Qiao et al., Reference Qiao, Li, Jiang, Lu, Franklin, Tang and Jiang2015 vs Tonkin et al., Reference Tonkin, Heino, Sundermann, Haase and Jähnig2016).
LCBD of flea assemblages were found to be mainly affected by composition of hosts that these assemblages exploit. As a result, regions with high compositional uniqueness of flea assemblages are also those with high compositional uniqueness of host assemblages. The role of among-region differences in host composition on differences in flea composition can be expected because the tight link between these metrics is well known for a variety of parasite taxa exploiting a variety of host taxa (e.g. Vinarski et al., Reference Vinarski, Korallo, Krasnov, Shenbrot and Poulin2007; Mihaljevic et al., Reference Mihaljevic, Hoye and Johnson2018) including fleas parasitic on small mammals (Krasnov et al., Reference Krasnov, Mouillot, Shenbrot, Khokhlova and Poulin2010). As a result, a particular host composition of a region and its representation across the entire study area determines the extent of contribution of a given parasite assemblage to total β-diversity. For example, higher LCBD values of flea assemblages were associated with host communities composed of gerbils, jerboas and ground squirrels, whereas the opposite was true for communities of bank voles and soricine shrews (see ‘Results’ section). The former host communities and their fleas are characteristic for arid habitats, whereas the latter host communities and their fleas are predominantly represented in forests. Arid habitats are less represented in the northern and temperate Palearctic than forests, so that communities of arid fleas harboured by arid hosts differ from the majority of other communities and thus contribute more to total β-diversity.
Obviously, host species are unevenly distributed across the continent and a particular host community is associated with a particular landscape. This may, at least partly, explain spatial patterns in LCBD. For example, comparison of model coefficients (Table 2) and geographic distribution of PCNMs (Fig. 5) suggests spatial structure of LCBD. At a relatively broader scale (PCNM2 in Table 2, Fig. 5A), higher LCBD values were mainly found in the predominantly western or predominantly eastern assemblages, whereas lower LCBD values were characteristic for the assemblages in the centre of the continent. At a somewhat finer scale (PCNM4 in Table 2, Fig. 5B), there was a trend of a south–north gradient in higher-to-lower LCBD values. This suggests that communities in western and eastern regions contained unique species not represented in the major part of the continent (e.g. western Atyphloceras nuperus and Ctenophthalmus bisoctodentatus and eastern Amphipsylla marikovskii and Ctenophthalmus congeneroides). The same was likely true for southern assemblages that contained species not distributed across most of the northern and temperate Palearctic (e.g. Xenopsylla conformis, Xenopsylla magdalinae, Mesopsylla hebes).
Surprisingly, flea species richness had no effect on LCBD values, although this effect has been found in studies of free-living organisms (Heino and Grönroos, Reference Heino and Grönroos2017; Kong et al., Reference Kong, Chevalier, Laffaille and Lek2017; Vilmi et al., Reference Vilmi, Karjalainen and Heino2017; da Silva et al., Reference da Silva, Hernández and Heino2018). Furthermore, the relationship between LCBD and species richness in different taxa and even in the same taxon in different habitats or sites could be either negative or positive or non-existent (Heino and Grönroos, Reference Heino and Grönroos2017; Kong et al., Reference Kong, Chevalier, Laffaille and Lek2017; Vilmi et al., Reference Vilmi, Karjalainen and Heino2017; da Silva et al., Reference da Silva, Hernández and Heino2018). Da Silva et al. (Reference da Silva, Hernández and Heino2018) suggested that the existence or the pronouncement of the LCBD–species richness relationship may depend on the proportions of rare and common species in the assemblages. This can be the case for fleas. In other words, LCBD of an assemblage can depend on a particular species composition rather than on a mere number of species.
In conclusion, SCBD of fleas was predictably associated with ecological/geographic but not morphological/life history traits. LCBD of flea assemblages was predictably associated with host composition, but not vegetation and climatic variables. Our results suggest that predictors of variation in SCBD and, especially, LCBD of parasites differ from those of free-living taxa due to substantial differences in the life history strategies. More studies on other parasites are required to support this conclusion and determine if parasitic taxa largely follow this general pattern.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182018001944.
Acknowledgements
This is publication number 988 of the Mitrani Department of Desert Ecology.
Author contributions
BRK and ISK conceive the study. GIS collected environmental information. BRK, ENS, SGM, NP, NE and BKK made measurements of body sizes. BRK, ISK, EMW and LVDM analysed data and drafted the manuscript. All authors participated in finalizing the manuscript.
Financial support
This study was supported by the Israel Science Foundation (grant no. 149/17 to BRK and ISK) and Section of General Biology, Department of Biological Sciences, Russian Academy of Sciences (project АААА-А17-117030310209-7 to SGM). EMW and LVDM received financial support from the Blaustein Center for Scientific Cooperation. ENS received the Scholarship of the President of the Russian Federation for short-term studies abroad.
Conflict of interest
None.
Ethical standards
Not applicable.
Author ORCIDs
Boris R. Krasnov 0000-0002-0382-3331
