Introduction
Species associations within a metacommunity (i.e. a set of biological communities occupying multiple localities and potentially linked by dispersal; Wilson, Reference Wilson1992; Leibold and Mikkelson, Reference Leibold and Mikkelson2002) may be either random or non-random (e.g. Diamond, Reference Diamond, Cody and Diamond1975; Connor and Simberloff, Reference Connor and Simberloff1979; Gotelli, Reference Gotelli2000; Gotelli and McCabe, Reference Gotelli and McCabe2002; Ulrich, Reference Ulrich2004; Gotelli and Ulrich, Reference Gotelli and Ulrich2010). Non-random species associations indicate that communities within a metacommunity are structured, that is, variation in species composition across communities is, to a certain degree, predictable, and communities are thus organized by certain assembly rules (e.g. Diamond, Reference Diamond, Cody and Diamond1975; Patterson and Atmar, Reference Patterson and Atmar1986). The non-randomness of species associations (=co-occurrence pattern) is tested using null models (Gotelli, Reference Gotelli2000). In particular, if the frequency of co-occurrence is greater than expected by chance, then these species are positively associated, whereas if this frequency is lower than expected by chance, then these species are negatively associated (Gotelli, Reference Gotelli2000; Gotelli and McCabe, Reference Gotelli and McCabe2002). Alternatively, if the frequency of co-occurrence does not differ from that expected by chance, then the community is randomly assembled.
The pattern of species co-occurrence allows inferring the mechanisms of community structuring. For example, predominantly positive co-occurrences suggest a mechanism resembling environmental filtering in which abiotic and/or biotic factors prevent the establishment or persistence of species in a particular location (Maire et al. Reference Maire, Gross, Börger, Proulx, Wirth, Pontes, Soussana and Louault2012; Kraft et al. Reference Kraft, Adler, Godoy, James, Fuller and Levine2015). As a result, an environmental filter allows a community to contain only species possessing certain traits that are necessary for persistence in that environment (e.g. Ingram and Shurin, Reference Ingram and Shurin2009). Consequently, co-occurring species are expected to be similar in their traits. Predominantly negative co-occurrences may indicate interspecific competition (Diamond, Reference Diamond, Cody and Diamond1975) or the ‘ghost of competition past’, both resulting in the dissimilarity of co-occurring species in their traits (limiting similarity; MacArthur and Levins, Reference MacArthur and Levins1967). The latter concept states that species can coexist only if their overlap in resource use is limited due to morphological, ecological or behavioural dissimilarities. However, Ulrich and Gotelli (Reference Ulrich and Gotelli2013) argued that negative co-occurrences may result from a change in species composition (species turnover) determined by limited dispersal and/or a gradient in habitat quality, independently of trait similarity or dissimilarity. Nevertheless, the majority of the above-cited studies suggest that studies aimed at inferring the community assembly rules from the pattern of species co-occurrences should incorporate information on species traits.
Although patterns of species co-occurrence have been studied in a variety of biological communities (e.g. Gotelli and McCabe, Reference Gotelli and McCabe2002; Jenkins, Reference Jenkins2006; Ulrich and Zalewski, Reference Ulrich and Zalewski2006; Veech, Reference Veech2006; Ulrich and Gotelli, Reference Ulrich and Gotelli2007; Ulrich et al. Reference Ulrich, Jabot and Gotelli2017; Freilich et al. Reference Freilich, Wieters, Broitman, Marquet and Navarrete2018; Harikrishnan and Vasudevan, Reference Harikrishnan and Vasudevan2018), species trait information has rarely been incorporated in these studies (Korňan et al. Reference Korňan, Svitok and Krištín2019; Legras et al. Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019; Vinarski et al. Reference Vinarski, Korallo-Vinarskaya, Shenbrot, Warburton, Surkova, Khokhlova and Krasnov2020). Furthermore, the absolute majority of species co-occurrence studies considered pairwise co-occurrences. However, not every pair of species in a community is positively or negatively associated, but rather, a general non-random pattern arises due to associations of a subset of species only (Sfenthourakis et al. Reference Sfenthourakis, Tzanatos and Giokas2006; Veech, Reference Veech2006; Gotelli and Ulrich, Reference Gotelli and Ulrich2010). Therefore, consideration of groupwise, rather than pairwise, co-occurrences will likely improve pattern detection and will make inferring assembly rules more reliable (Morueta-Holme et al. Reference Morueta-Holme, Blonder, Sandel, McGill, Peet, Ott, Violle, Enquist, Jørgensen and Svenning2016; Legras et al. Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019). This can be achieved by applying a network approach to species co-occurrence studies. Recently, Legras et al. (Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019) proposed a new methodological approach based on the concepts of network theory and modularity and aimed to disentangle the deterministic and stochastics drivers of species co-occurrence patterns accounting for multiple functional traits. In brief, this approach involves calculating modularity in (a) a site × species network and (b) a trait × species network and then assessing the congruence between the 2 networks by testing whether the species from a given functional module belong to the same co-occurrence module. Legras et al. (Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019) demonstrated that this approach can accurately detect assembly rules acting at different spatial scales.
The scale-dependence (spatial or temporal) of co-occurrence patterns has been found in multiple taxa (Gotelli and Ellison, Reference Gotelli and Ellison2002; Sanders et al. Reference Sanders, Gotelli, Wittman, Ratchford, Ellison and Jules2007; Harms and Dinsmore, Reference Harms and Dinsmore2016; Legras et al. Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019). However, one of the methodological difficulties in community ecology studies is to identify community boundaries, which often cannot be unequivocally defined (Strayer et al. Reference Strayer, Power, Fagan, Pickett and Belnap2003). This is especially true for communities at different scales. From this perspective, parasite communities present convenient models for elucidating the patterns of species co-occurrence and underlying processes. This is because parasite community boundaries, both in space and time, can be identified more easily than those of free-living species, due to the fact that parasites inhabit a set of ‘islands’ represented by their hosts, whereas the environment between these ‘islands’ is mostly unfavourable. Furthermore, parasite communities are fragmented between host individuals (=infracommunities), populations of conspecific hosts (=component communities), and communities of multiple host species (=compound communities) (Holmes and Price, Reference Holmes, Price, Kittawa and Anderson1986; Poulin, Reference Poulin2007).
A pioneering study of species co-occurrence in parasite communities dealt with infracommunities of ectoparasites in marine fishes and reported an absence of any strong pattern, thus suggesting that these communities are mostly randomly assembled (Gotelli and Rohde, Reference Gotelli and Rohde2002), although the opposite appeared to be the case for ectoparasites of freshwater fishes (Bellay et al. Reference Bellay, Takemoto and Oliveira2012). Studies of co-occurrence patterns in ectoparasitic arthropods of terrestrial hosts reported predominantly positive species associations in different taxa, host species, and geographic regions (Krasnov et al. Reference Krasnov, Stanko and Morand2006a, Reference Krasnov, Matthee, Lareschi, Korallo-Vinarskaya and Vinarski2010, Reference Krasnov, Shenbrot and Khokhlova2011, Reference Krasnov, Shenbrot, Korallo-Vinarskaya, Vinarski, van der Mescht, Warburton and Khokhlova2019; Presley, Reference Presley2007; Tello et al. Reference Tello, Stevens and Dick2008; Veitch et al. Reference Veitch, Bowman and Schulte-Hostedde2020). By comparison, to the best of our knowledge, only one study has incorporated information on parasite traits in investigating their co-occurrence patterns (Vinarski et al. Reference Vinarski, Korallo-Vinarskaya, Shenbrot, Warburton, Surkova, Khokhlova and Krasnov2020). This study compared the probability of pairwise co-occurrences and pairwise trait dissimilarity in communities of fleas and gamasid mites, parasitic on small mammals. A significant, albeit weak, tendency of flea, but not mite, communities to be composed of functionally similar species has been found. However, using pairwise co-occurrences and trait similarity could bias the results due to the above-mentioned reason that a detection of non-randomness could be associated with only a subset of the entire community species pool (e.g. Sfenthourakis et al. Reference Sfenthourakis, Tzanatos and Giokas2006).
Here, we aimed to remedy this shortcoming and applied the framework proposed by Legras et al. (Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019) (see above) to investigate the relationships between species co-occurrence and trait similarity in flea communities at various scales (compound communities: across regions within a biogeographic realm or across sampling sites within a geographic region; component communities: across sampling sites within a geographic region; and infracommunities: within a sampling site). Fleas are holometabolous haematophagous insects, being the most diverse on small and medium-sized burrowing mammals. In the absolute majority of species, larvae are not parasitic and feed on all kinds of organic matter found in the hosts’ burrows/nests where pre-imaginal development takes place (Marshall, Reference Marshall1981; Krasnov, Reference Krasnov2008). Based on the predominantly aggregative structure of flea communities found in earlier studies (e.g. Krasnov et al. Reference Krasnov, Stanko and Morand2006a, Reference Krasnov, Shenbrot and Khokhlova2011) and the link between pairwise co-occurrence and trait similarity, we expected to find congruence between modules of co-occurrence and functional networks (see below definition of modules) at all hierarchical scales. For compound communities, we considered flea species co-occurrences not only among regions or sites but also among host species. In addition, we inferred the main processes (deterministic ones, such as limiting similarity/competition and environmental filtering, or stochastic ones) acting on community assembly.
Materials and methods
Data on flea species distribution
Data on flea species distribution across regions and host species within biogeographic realms (the Afrotropics, Australasia, Indo-Malay, Nearctic, Neotropics, and Palearctic) were obtained from various literature sources. At this scale (i.e. within a biogeographic realm), only host species that harboured fleas and belonged to tachyglossid Monotremata, Dasyuromorphia, Paramelemorphia, Diprotodontia, Macropodiformes, Didelphimorphia, Paucituberculata, Macroscelidea, Eulipotyphla, Rodentia and the ochotonid Lagomorpha were considered. This study included data on the distribution of 1313 flea species across 1153 mammal host species in 15 Afrotropical regions, 8 Australasian regions, 10 Indo-Malayan regions, 23 Nearctic regions, 17 Neotropical regions and 36 Palearctic regions (see lists of regions, maps and references in Krasnov et al. Reference Krasnov, Shenbrot and Khokhlova2022 and lists of flea and host species in the respective publications).
At the scale of regional compound communities, component communities and infracommunities, the sources of data were surveys of fleas on individual small mammals carried out in 1981–1984 across 83 sampling sites in Mongolia (Kiefer et al. Reference Kiefer, Klimaszewski and Krumpal1982, Reference Kiefer, Krumpal, Cendsuren, Lobachev and Khotolkhu1984, Reference Kiefer, Goncharov and Lobachev1990); in 1998–2015 across 38 sampling sites in Northwest Argentina (López Berrizbeitia, Reference López Berrizbeitia2018; López Berrizbeitia and Díaz, Reference López Berrizbeitia and Díaz2019); in 2006–2013 across 22 sites in Argentinian Patagonia (Sanchez, Reference Sanchez2012; Sanchez and Lareschi, Reference Sanchez and Lareschi2013, Reference Sanchez and Lareschi2019); in 1986–2000 across 13 sites in eastern Slovakia (Stanko et al. Reference Stanko, Miklisová, Gouy de Bellocq and Morand2002); in 1988 across 49 sites in Western Siberia [compiled from the database of the Omsk Research Institute of Natural Foci Infections (Omsk, Russian Federation) by NPK-V] (Surkova et al. Reference Surkova, Korallo-Vinarskaya, Vinarski, van der Mescht, Warburton, Khokhlova and Krasnov2018); in 2003–2018, across 36 sites in southwestern South Africa (Matthee et al. Reference Matthee, Horak, Beaucournu, Durden, Ueckermann and McGeoch2007; van der Mescht and Matthee, Reference van der Mescht and Matthee2017; Stevens et al. Reference Stevens, Stekolnikov, Ueckermann, Horak and Matthee2022); and in 2005–2008 across 18 sites in the Lushoto District, Usambara Mountains, Tanzania (Laudisoit et al. Reference Laudisoit, Leirs, Makundi and Krasnov2009a, Reference Laudisoit, Neerinckx, Makundi, Leirs and Krasnov2009b). Lists of flea and mammal species in each region can be found in Supplementary Material, Appendix 1, Table S1.
At the scale of component communities, 31 host species harbouring at least 8 flea species were selected from the regional data (see Supplementary Material, Appendix 1, Table S2). At the scale of infracommunities, 25 host species harbouring at least 6 flea species at least 1 sampling site were selected from the regional data (see Supplementary Material, Appendix 1, Table S3). The latter analyses were carried out for conspecific host individuals captured at the same sampling site at the same time.
Data on flea traits
Although a trait, by definition, is ‘any morphological, physiological or phenological feature measurable at the individual level’ (Violle et al. Reference Violle, Navas, Vile, Kazakou, Fortunel, Hummel and Garnier2007), information on parasite traits, measurable at the individual level, is largely unavailable. Consequently, we used species-specific attributes or characteristics that we further refer to as traits for simplicity. As such, species-specific attributes have been used in earlier studies of functional diversity (e.g. feeding habit in Pavoine et al. Reference Pavoine, Marcon and Ricotta2016). In the current study, each flea species was characterized by 4 ecological and 2 morphological traits. Ecological traits included the total number of host species exploited across a flea's geographic range; these hosts’ phylogenetic diversity; the latitudinal span of geographic range; and microhabitat preference as the relative time spent either in the hair of the host or in its/their nest/burrow (preference for hair or for nest, or else no clear preference). Morphological traits included the number and possession of sclerotized ctenidia and body size. The rationale for selecting these traits, information sources and calculation details can be found elsewhere (Krasnov et al. Reference Krasnov, Shenbrot, Khokhlova, Stanko, Morand and Mouillot2015, Reference Krasnov, Grabovsky, Khokhlova, López Berrizbeitia, Matthee, Roll, Sanchez, Shenbrot and van der Mescht2024). Data on the latitudinal span of geographic range were not available for the Australasian and the Indo-Malayan fleas, and data on microhabitat preference were not available for the Australasian fleas.
Data analyses
As mentioned above, the main idea of Legras et al. (Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019) is to compare the functional and the co-occurrence networks when the former is based on species traits and the latter is based on the number of species co-occurrences. In this case, the edges within each network represent similarity in either species traits or species co-occurrences, respectively. In a modularly structured network, species from some distinct groups ( = modules) are connected mainly with each other, while they have much fewer connections to other species belonging to other modules. Consequently, modularity measures the density of edges inside groups compared to edges between groups (Newman, Reference Newman2006). In the framework of Legras et al. (Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019), each network's modularity is calculated, and the modules’ species compositions are compared between the 2 networks, that is, testing whether species from a given functional module also belong to a co-occurrence module. If this is the case, then functionally similar species inhabit/are found in the same sites/regions/host species. This would strongly suggest that environmental or host-associated factors act as filters determining species composition. On the contrary, high dissimilarity between the 2 networks would suggest that functionally similar species are distributed among different sites/regions/host species, likely resulting from interspecific competition that, in turn, leads to reduced similarity between co-occurring species (i.e. limiting similarity sensu MacArthur and Levins, Reference MacArthur and Levins1967). However, a congruence between the module species composition of the functional and the co-occurrence networks that is not especially high and not especially low would indicate that species distribution among sites/regions/host species does not differ from random. This situation results either from stochastic mechanisms (e.g. dispersal limitations and/or population dynamics; Hubbell, Reference Hubbell2001) or from the fact that deterministic processes (environmental filtering and competition) act simultaneously, affecting species co-occurrences in the opposite directions. The latter may occur if the distribution of some species in a metacommunity is caused by environmental/host-associated filters, whereas the distribution of other species is determined by competition. In such a case, there is no clear dominance of one of the processes, so they cannot be disentangled.
In our study, the approach of Legras et al. (Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019) was slightly modified. Although the construction of co-occurrence networks and identification of co-occurrence modules followed Legras et al. (Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019), the construction of functional networks and identification of functional modules (=functional groups) was carried out using the best partitioning of the functional trait space and followed Rubio and Swenson (Reference Rubio and Swenson2022).
A functional matrix (flea species × traits with traits in columns) was compiled for all 1362 fleas in the dataset. Two presence/absence co-occurrence matrices for compound communities (either realm-specific or regional) and 1 presence/absence co-occurrence matrix for either component communities or infracommunities were constructed (see scheme in Supplementary Material, Appendix 2). For each compound community at the scale of a biogeographic realm (i.e. within a realm), these matrices were (a) a spatial co-occurrence matrix: flea species × regions (regions in columns) and (b) a host co-occurrence matrix: flea species × host species (host species in columns). For each compound community at the regional scale (i.e. within a region), these matrices were (a) a spatial co-occurrence matrix: flea species × sampling sites (sites in columns) and (b) a host co-occurrence matrix: flea species × host species (host species in columns). For each set of component communities (i.e. within a host species within a region across sampling sites), a spatial co-occurrence matrix was flea species × sampling sites (sites in columns). For each set of infracommunities (e.g. within a host species at the same sampling site across host individuals), a host co-occurrence matrix was flea species × host individuals. Although Legras et al. (Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019) recommended constructing the co-occurrence matrices using species abundances, presence/absence data appeared to be more reliable than abundance data for parasites (Gotelli and Rohde, Reference Gotelli and Rohde2002). This is because (a) parasites are aggregated across host individuals, being the most abundant only in a few host individuals, whereas the majority of host individuals harbour only a few or no parasites (Anderson and May, Reference Anderson and May1978); (c) ectoparasite counts are an unreliable measure of their abundance (see Krasnov et al. Reference Krasnov, Shenbrot, Khokhlova, Hawlena and Degen2006b), and (b) comparison of the results of analyses using abundance vs presence/absence data in studies of parasite community ecology, demonstrated similar performances of these analyses (Haukisalmi and Henttonen, Reference Haukisalmi and Henttonen1993; Krasnov et al. Reference Krasnov, Spickett, Junker, Bugmyrin, Ieshko, Bespyatova, Stanko, Khokhlova and Matthee2021; but see Brian and Aldridge, Reference Brian and Aldridge2021).
Then, a similarity matrix was calculated for each above-mentioned matrix, that is, a matrix of trait (=functional) similarity between species from a functional matrix and a matrix(ces) of co-occurrence similarity from co-occurrence matrix(ces). Functional similarity was calculated as (1-Gower distance index) for each species pair, using the ‘FD’ package (Laliberté, and Legendre, Reference Laliberté and Legendre2010; Laliberté et al. Reference Laliberté, Legendre and Shipley2014), implemented in the R statistical environment (R Core Team, 2023). Co-occurrence similarity was calculated as (1-binary version of the Bray–Curtis dissimilarity index), using the R package ‘vegan’ (Oksanen et al. Reference Oksanen J, Simpson, Blanchet, Kindt, Legendre, Minchin, O'Hara, Solymos, Stevens, Szoecs, Wagner, Barbour, Bedward, Bolker, Borcard, Carvalho, Chirico, De Caceres, Durand, Evangelista, FitzJohn, Friendly, Furneaux, Hannigan, Hill, Lahti, McGlinn, Ouelette, Ribeiro Cunha, Smith, Stier, Ter Braak and Weedon2022). The values of these similarity matrices were then used as weights of the edges between species (nodes) in the functional or co-occurrence networks.
For a functional or a co-occurrence network, modules were identified. Modules are distinct subsets of species that are mainly connected with each other, while they have far fewer connections to other species belonging to other modules (Girvan and Newman, Reference Girvan and Newman2002). For functional network, modules were considered to be equivalent to functional groups and identified from the functional similarity matrix using the ‘Nbclust’ function (with options distance = ‘euclidean’ and method = ‘ward.D’) in the R package ‘Nbclust’ (Charrad et al. Reference Charrad, Ghazzali, Boiteau and Niknafs2014). This function calculates the best partitioning of the functional trait space into functional groups (=modules) comparing results when using 30 different indices. For co-occurrence networks, identification of modules was done via the Louvain optimization algorithm, which is recognized as one of the best algorithms for module detection (see Legras et al. Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019), using the R package ‘modMax’ (Schelling and Hui, Reference Schelling and Hui2015).
Finally, the congruence between the functional and the co-occurrence networks was assessed by an index of module diversity for each functional group of Gauzens et al. (Reference Gauzens, Thébault, Lacroix and Legendre2015), modified by Legras et al. (Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019). This index was calculated for each group of species comprising a functional module, using the R code of Legras et al. (Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019). It varies from zero to unity, being zero if all species of a functional module tend to belong to the same co-occurrence module and unity if all species of a functional module tend to belong to different co-occurrence modules. The index values of all functional modules of a network were averaged, resulting in the DgM index. Subsequently, the value of the DgM index for each compound, component, or infracommunity was compared with distribution of the DgM values from 999 null models, in which the partition into functional modules and the number and size of co-occurrence modules were the same as in the original data, but species distribution among co-occurrence modules was random. The observed DgM index was compared with the distribution of DgM indices derived from null models. A significantly lower observed DgM value than expected by chance (<5% of the null distribution) indicates the predominance of environmental/host-associated filtering, whereas a significantly higher value than expected by chance (>95% of the null distributions) suggests the effects of competition/limiting similarity (Legras et al. Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019). An observed DgM value between 5% and 95% of the null DgM distribution suggests either the effect of stochastic mechanisms or the lack of clear dominance of one of the deterministic processes (Legras et al. Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019). The R code for null modelling was taken from Legras et al. (Reference Legras, Loiseau, Gaertner, Poggiale, Ienco, Mazouni and Mérigot2019). When the value of DgM was calculated for a given community, only modules containing species from this community were retained in a functional network.
Results
The best partitioning of the functional trait space resulted in the 21 functional groups (=functional modules) with average pairwise within-module similarity of 0.68 and between-module similarity of 0.59. Within-module functional similarity ranged from 0.70 to 0.93 (average 0.83 ± 0.04) for compound communities within biogeographic realms, from 0.70 to 0.82 (average 0.77 ± 0.01) for compound communities within geographic regions, from 0.62 to 0.95 (average 0.75 ± 0.01) for component communities, and from 0.62 to 0.98 (average 0.74 ± 0.02) for infracommunities (Tables 1–2, Tables S4–S5 in Supplementary Material, Appendix 1). The respective values of between-module similarity were 0.59–0.73 (average 0.65 ± 0.01), 0.60–0.65 (average 0.63 ± 0.01), 0.50–0.71 (average 0.62 ± 0.01) and 0.40–0.77 (average 0.60 ± 0.01) (Tables 1–2, Tables S4–S5 in Supplementary Material, Appendix 1).
Table 1. Congruence between functional (F) and co-occurrence (CoR: across regions; CoH: across host species) networks of compound flea communities in 6 biogeographic realms
M: modularity value (for co-occurrence networks), WS/BS: average within- and between-module similarity (for functional networks), Modules: number of detected modules, DgM: index of congruence (see text for explanation). p: proportion of DgM values from null models that are lower than the observed DgM; Process: the most likely process affecting community assembly inferred from comparison of the observed and null DgM values (EF, environmental filtering; HF, host-associated filtering; S, stochastic process(es) or no clear dominance of a deterministic process).
Table 2. Congruence between functional (F) and co-occurrence (CoL: across sampling sites; CoH: across host species) networks of compound flea communities in 7 regions
M: modularity value, WS/BS: average within- and between-module similarity (for functional networks), Modules: number of detected modules, DgM: index of congruence (see text for explanation). p: proportion of DgM values from null models that are lower than the observed DgM; Process: the most likely process affecting community assembly inferred from comparison of the observed and null DgM values (EF, environmental filtering; HF, host-associated filtering; S, stochastic process(es) or no clear dominance of a deterministic process).
Modularity in spatial co-occurrence networks ranged from 0.14 to 0.57 (average 0.34 ± 0.07) for compound communities within biogeographic realms, from 0.06 to 0.23 (average 0.16 ± 0.02) for compound communities within geographic regions, and from 0.01 to 0.28 (average 0.31 ± 0.01) for component communities (Tables 1–2, Tables S4–S5 in Supplementary Material, Appendix 1). Modularity values in networks of flea co-occurrences in host species (compound communities) or individuals (infracommunities) ranged from 0.14 to 0.35, from 0.06 to 0.15, and from 0.01 to 0.46 for compound communities at the scale of biogeographic realms or geographic regions and infracommunities, respectively (average 0.27 ± 0.03, 0.09 ± 0.01 and 0.20 ± 0.03, respectively) (Tables 1–2, Tables S4–S5 in Supplementary Material, Appendix 1).
The number of modules in functional networks ranged from 7 to 17 in realm-specific networks, from 5 to 13 in regional networks, from 2 to 8 in networks of component communities, and from 2 to 6 in infracommunities (Tables 1–2, Tables S5–S6 in Supplementary Material, Appendix 1). The Louvain algorithm detected 2–8 modules in spatial co-occurrence networks, and 2–12 modules in the networks of co-occurrence in host species/individuals (Tables 1–2, Tables S5–S6 in Supplementary Material, Appendix 1). In general, the number of detected modules in the co-occurrence networks was not associated with the number of species in a network, whereas the higher number of functional modules was detected in larger flea communities (see illustrative examples of module distribution in Figs 1–2 for compound communities and Figs. S1–S2 in Supplementary Material, Appendix 3 for component communities and infracommunities).
Figure 1. Modules based on trait similarity, spatial (across regions) co-occurrence similarity, and host species co-occurrence similarity for compound communities of fleas in the Palearctic. The number inside or near the circle is the number of species in the module. In trait-associated modules, the number in parentheses is average within-module similarity (above line) and between-module similarity (below line) between pairs of species. Edge width is proportional to average similarity between species belonging to the modules.
Figure 2. Modules based on trait similarity, spatial (across localities) co-occurrence similarity, and host species co-occurrence similarity for compound communities of fleas in Mongolia. The number inside or near the circle is the number of species in the module. In trait-associated modules, the number in parentheses is average within-module similarity (above line) and between-module similarity (below line) between pairs of species. Edge width is proportional to average similarity between species belonging to the modules.
Comparisons of the observed DgM values with the distributions of those obtained from null modelling demonstrated that functional and spatial or host-associated co-occurrence networks of fleas in compound communities within biogeographic realms were consistently congruent, except for spatial co-occurrence in Australasia and host-associated co-occurrence in Indo-Malay, and suggested that these communities (both across regions and across host species) were assembled via environmental or host-associated filtering, respectively (Table 1). In contrast, functional and spatial/host-associated co-occurrence networks, at the scale of regional compound communities (i.e. across localities), were mostly incongruent (except the spatial co-occurrence network in Slovakia and Tanzania and the host-associated co-occurrence network in Mongolia, both suggesting environmental/host-associated filtering) (Table 2). The incongruence between the networks indicated either stochastic processes or the lack of dominance of any deterministic process. Analyses of congruence between functional and either spatial (for component communities) or host-associated (for infracommunities) co-occurrence networks demonstrated that assembly rules in these communities varied among the host species harbouring them. In 18 of 31 component communities, stochastic processes in their assembly were indicated, whereas environmental filtering seemed to dominate in 4, and limiting similarity/competition in 9 communities (Fig. 3, see detailed results in Table S4 in Supplementary Material, Appendix 1). Stochastic assembly processes dominated in infracommunities harboured by 9 of 25 host species, limiting similarity/competition in infracommunities of 13 host species, and host-associated filtering in the infracommunities of 3 hosts (Microtus oeconomus, Sorex araneus and Lophuromys kilonzoi) (Fig. 4, see detailed results in Table S5 in Supplementary Material, Appendix 1).
Figure 3. Proportions of DgM values from null models that are lower than the observed DgM calculated for assessing congruence between functional and co-occurrence networks in component communities of fleas harboured by 31 host species. The most likely process affecting community assembly, as inferred from the comparison of the observed and null DgM values, is shown by bar colour (black: environmental filtering, grey: limiting similarity/competition, white: stochastic process(es) or no clear dominance of a deterministic process). M, Mongolia; P, Patagonia, WS, Western Siberia, Sl, Slovakia; SA, South Africa; T, Tanzania.
Figure 4. Proportions of DgM values from null models that are lower than the observed DgM calculated for assessing congruence between functional and co-occurrence networks in infracommunities of fleas harboured by 25 host species. The most likely process affecting community assembly as inferred from the comparison of the observed and null DgM values is shown by bar colour (black: host-associated filtering, grey: limiting similarity/competition, white: stochastic process(es) or no clear dominance of a deterministic process). M, Mongolia; P, Patagonia; WS, Western Siberia; Sl, Slovakia; SA, South Africa; T; Tanzania.
Discussion
We found substantial differences in the patterns of species co-occurrence and the processes dominating flea community assembly at different hierarchical scales. Moreover, these patterns and processes were consistent at some scales and contingent at other scales. In particular, regional compound communities within biogeographic realms were characterized by positive species associations across both regions and host species, suggesting the important role of environmental or host-associated filtering. On the contrary, local compound communities within a relatively homogeneous region demonstrated random species associations and, thus, indicating the role of stochastic factors. In component and infracommunities, patterns of species associations and, consequently, presumed assembly rules varied greatly among communities harboured by different host species. This supported the results of some, but not other, earlier studies on co-occurrence patterns in flea communities (Krasnov et al. Reference Krasnov, Matthee, Lareschi, Korallo-Vinarskaya and Vinarski2010, Reference Krasnov, Shenbrot and Khokhlova2011, Reference Krasnov, Pilosof, Stanko, Morand, Korallo-Vinarskaya and Vinarski2015 vs Krasnov et al. Reference Krasnov, Stanko and Morand2006a, Reference Krasnov, Pilosof, Stanko, Morand, Korallo-Vinarskaya and Vinarski2014; Vinarski et al. Reference Vinarski, Korallo-Vinarskaya, Shenbrot, Warburton, Surkova, Khokhlova and Krasnov2020; Gibert et al. Reference Gibert, Shenbrot, Stanko, Khokhlova and Krasnov2021, respectively).
The mechanism behind the aggregative pattern of flea co-occurrences in regional compound communities within a biogeographic realm, coupled with the trait similarity of co-occurring species, is associated with similar environmental and host-associated preferences of phylogenetically close flea species possessing similar traits (Krasnov et al. Reference Krasnov, Shenbrot, Khokhlova, Stanko, Morand and Mouillot2015). The environmental factors determining flea preferences are, first and foremost, air temperature and relative humidity (to which fleas are highly sensitive; Margalit and Shulov, Reference Margalit and Shulov1972), as well as vegetation structure (see review in Krasnov, Reference Krasnov2008). These factors affect the microclimate in hosts’ burrows where pre-imagoes develop, and adult fleas spend substantial time. An additional environmental factor is soil structure since soil particles are used by flea larvae for cocoon construction (Marshall, Reference Marshall1981). Particles of different soil types likely differ in their water-absorption abilities, while active water vapour uptake is highly important for cocooned fleas (Rudolph and Knulle, Reference Rudolph and Knulle1982). In Tanzania, flea abundance was found to be affected by soil properties such as exchangeable calcium and magnesium, phosphorus availability and base saturation (Meliyo et al. Reference Meliyo, Kimaro, Msanya, Mulungu, Hieronimo, Kihupi, Gulinck and Deckers2014). Obviously, environmental differences between regions of the same biogeographic realm are sharp, resulting in further classification into biomes and/or ecoregions (e.g. Olson et al. Reference Olson, Dinerstein, Wikramanayake, Burgess, Powell, Underwood, D'amico, Itoua, Strand, Morrison, Loucks, Allnutt, Ricketts, Kura, Lamoreux, Wettengel, Hedao and Kassem2001). Furthermore, regions within a realm differ not only in their abiotic environment but also in the species composition of small mammals that fleas can utilize as their hosts. The process of host selection by parasites is largely determined by the complementarity of co-evolving host and parasite traits (McQuaid and Britton, Reference McQuaid and Britton2013; dos Santos Cardoso et al. Reference dos Santos Cardoso, de Andreazzi, Maldonado and Gentile2021; see Krasnov et al. Reference Krasnov, Shenbrot, Khokhlova and Degen2016 for fleas). Host–parasite co-evolution consists of parasite adaptation to extract resources from their hosts, as well as to evade or overcome their defences, and reciprocal host counter-adaptation to avoid or tolerate infection (e.g. Buckingham and Ashby, Reference Buckingham and Ashby2022). Complementarity or matching between parasite and host traits may result in establishing novel parasite–host associations when a host or a parasite invades new areas (e.g. Schatz and Park, Reference Schatz and Park2021). Complementary parasite–host traits cause host species composition to act as a filter for assembly of a regional parasite community. The reason for not detecting environmental or host-associated filtering processes as assembly mechanisms, for compound flea communities in the Australasian and Indo-Malayan realms, respectively, could be the much poorer knowledge of either geographic distribution of fleas or host–flea associations in these than in other realms. For example, only 2 flea species parasitic on small mammals were recorded in Northern Territories of Australia (Dunnet and Mardon, Reference Dunnet and Mardon1974). Completely different types of analyses, carried out on almost the same data as in the present study, suggested the predominance of the dispersal-based mechanism of flea compound community assembly within biogeographic realms realized via host dispersal (Gibert et al. Reference Gibert, Shenbrot, Stanko, Khokhlova and Krasnov2021). However, these analyses did not consider trait difference/similarity between species and were based on a pairwise comparison of communities between continental sections/landscapes rather than between regions.
In contrast to regional compound communities within a biogeographic realm, functional and co-occurrence networks, at the scale of local communities within a region, were mostly incongruent. On the one hand, this could indicate stochastic birth/death processes of community assembly, resulting in ‘neutral’ communities that are primarily shaped by population dynamics and/or dispersal limitations (Hubbell, Reference Hubbell2001). On the other hand, counterbalancing of environmental filtering and limiting similarity/competition may lead to a pattern ostensibly resembling neutrality (Purves and Pacala, Reference Purves, Pacala, Burslem, Pinard and Hartley2005). The latter seems to be more likely because at least some, albeit only two, compound communities (Slovakian and Mongolian) suggested the acting of deterministic processes. The lack of a clear indication of any deterministic process at this scale may also be associated with relative environmental homogeneity within a region, so that environmental differences, as well as differences in host species composition, were not sharp enough for environmental/host-associated filtering to be detected (but see McNew et al. Reference McNew, Barrow, Williamson, Galen, Skeen, DuBay, Gaffney, Johnson, Bautista, Ordoñez, Schmitt, Smiley, Valqui, Bates, Hackett and Witt2021). Another reason for the lack of detection of deterministic processes at this scale is that compound communities across localities within a region are represented by the set of communities harboured by different host species (i.e. component communities). Processes affecting these communities’ assembly varied both between regions and between host species. Variation in the degree of congruence between flea functional and co-occurrence networks and flea community assembly processes could stem from ecological and/or physiological differences between hosts, such as mobility patterns, burrow/nest structure, population density and immunocompetence as well as from differences in relative numbers of flea species with a certain set of traits. For example, if only a small number of flea species possess certain traits allowing them to overcome host defences and exploit strongly immunocompetent hosts (Møller et al. Reference Møller, Christe and Garamszegi2005) (which the remaining fleas from the same community cannot), then the associations between these and the remaining fleas would be negative, leading to a pattern superficially resembling the effect of limiting similarity/competition. If, however, the number of flea species capable of overcoming host defences is high, then positive associations between these species would be found, and a process resembling environmental filtering would be detected. Therefore, the manifestation of a given assembly rule in parasite component communities might be the result of interplay between parasite traits, host traits and parasite community composition. It is well known that the pattern of species association within a metacommunity varies between pairs or groups of species, from positive to random to negative (Sfenthourakis et al. Reference Sfenthourakis, Tzanatos and Giokas2006; D'Amen et al. Reference D'Amen, Mod, Gotelli and Guisan2018). The conclusion regarding the predominant pattern of species associations and, consequently, the rule of community assembly would depend on the relative number of species pairs/groups with a given association pattern.
Functional and co-occurrence networks in flea infracommunities were congruent in some, but incongruent in other, host species. Among 16 hosts in the infracommunities in which between-network congruence was found, a limiting similarity/competition process was detected in 13 hosts and environmental filtering in 3 hosts. The occurrence of competition among fleas exploiting the same host individual contradicts studies on the patterns of flea co-occurrences (Krasnov et al. Reference Krasnov, Stanko and Morand2006a, Reference Krasnov, Matthee, Lareschi, Korallo-Vinarskaya and Vinarski2010) but supports experimental and observational studies that reported possible interspecific competition (Day and Benton, Reference Day and Benton1980; Krasnov et al. Reference Krasnov, Burdelova, Khokhlova, Shenbrot and Degen2005; Surkova et al. Reference Surkova, Korallo-Vinarskaya, Vinarski, van der Mescht, Warburton, Khokhlova and Krasnov2018). Interestingly, evidence for the interspecific flea competition found in experimental studies has suggested that it (a) is asymmetric and (b) occurs among larvae rather than imago fleas. Intraspecific larval competition may or may not be further manifested in the pattern of imago co-occurrence. Moreover, parasite infracommunities are ephemeral not only from the perspective of a host individual's lifespan (Poulin, Reference Poulin2007) but, in the case of fleas, from the fact that infracommunities’ presence and composition are temporally variable and may change daily (Krasnov et al. Reference Krasnov, Shenbrot, Khokhlova, Hawlena and Degen2006b). This explains the variation in the detection of deterministic and stochastic processes in infracommunity assembly.
In conclusion, assembly processes in parasite communities appeared to be scale-dependent, with different mechanisms acting at different scales (see also McNew et al. Reference McNew, Barrow, Williamson, Galen, Skeen, DuBay, Gaffney, Johnson, Bautista, Ordoñez, Schmitt, Smiley, Valqui, Bates, Hackett and Witt2021). The main reason for this is that community composition at the largest scale is predominantly affected by evolutionary, historical and biogeographical processes, whereas community composition at smaller scales is mainly determined by ecological and demographic processes. The approach combining species traits and patterns of groupwise, rather than pairwise, co-occurrence improves the reliability of detection and understanding of processes of community assembly. Regarding parasites, this approach could also be helpful from the perspective of the dynamics of vector-borne diseases (Makundi et al. Reference Makundi, Massawe, Borremans, Laudisoit and Katakweba2015).
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182024000969.
Data availability statement
Raw data on flea and host species at the scale of biogeographic realms are contained in the sources cited in Krasnov et al. (Reference Krasnov, Shenbrot and Khokhlova2022). The remaining data can be obtained from the corresponding author upon request.
Acknowledgements
We thank Geut Galai for help with the R code for producing figures. AL wishes to acknowledge the financial support from the Belgium funds for Scientific Research (FNRS and FRIA), the University of Antwerp (Belgium), the Sokoine University of Agriculture, Morogoro, Tanzania, and the logistical support from the Pest Management Centre that enabled the data collection in Tanzania. We thank 2 anonymous referees for helpful comments on the earlier version of the manuscript.
Author contributions
BRK conceived and designed the study. All authors collected the data. BRK performed statistical analyses and wrote the first draft of the article. All authors finalized the article.
Financial support
This study was partly supported by the Israel Science Foundation (grant 548/23 to BRK and ISK).
Competing interests
The authors declare there are no conflicts of interest.
Ethical standards
This study is based on published data, and therefore, ethical standards are not applicable.
Open access
