An eco-epidemiological modeling approach to investigate dilution effect in two different tick-borne pathosystems.

Disease (re)emergence appears to be driven by biodiversity decline and environmental change. As a result, it is increasingly important to study host-pathogen interactions within the context of their ecology and evolution. The dilution effect is the concept that higher biodiversity decreases pathogen transmission. It has been observed especially in zoonotic vector-borne pathosystems, yet evidence against it has been found. In particular, it is still debated how the community (dis)assembly assumptions and the degree of generalism of vectors and pathogens affect the direction of the biodiversity-pathogen transmission relationship. The aim of this study was to use empirical data and mechanistic models to investigate dilution mechanisms in two rodent-tick-pathogen systems differing in their vector degree of generalism. A community was assembled to include ecological interactions that expand from purely additive to purely substitutive. Such systems are excellent candidates to analyze the link between vector ecology, community (dis)assembly dynamics, and pathogen transmission. To base our mechanistic models on empirical data, rodent live-trapping, including tick sampling, was conducted in Wales across two seasons for three consecutive years. We have developed a deterministic single-vector, multi-host compartmental model that includes ecological relationships with non-host species, uniquely integrating theoretical and observational approaches. To describe pathogen transmission across a gradient of community diversity, the model was populated with parameters describing five different scenarios differing in ecological complexity; each based around one of the pathosystems: Ixodes ricinus (generalist tick)-Borrelia burgdorferi and I. trianguliceps (small mammals specialist tick)-Babesia microti. The results suggested that community composition and interspecific dynamics affected pathogen transmission with different dilution outcomes depending on the vector degree of generalism. The model provides evidence that dilution and amplification effects are not mutually exclusive in the same community but depend on vector ecology and the epidemiological output considered (i.e., the "risk" of interest). In our scenarios, more functionally diverse communities resulted in fewer infectious rodents, supporting the dilution effect. In the pathosystem with generalist vector we identified a hump shaped relationship between diversity and infections in hosts, while for that characterized by specialist tick, this relationship was more complex and more dependent upon specific parameter values.

INTRODUCTION

The dilution effect is the concept that higher biodiversity decreases pathogen transmission in multi‐host/multi‐vector pathosystems as different hosts/vectors present different degrees of transmission competence (Ostfeld & Keesing, 2000; Schmidt & Ostfeld, 2001). This phenomenon has been mainly investigated in the context of zoonotic vector‐borne pathogens, where higher biodiversity seems to reduce human disease risk (LoGiudice et al., 2003; Wood & Lafferty, 2013). Studies focusing on the interactions between pathogen transmission, disease risk and biodiversity have gained progressively more attention due to the current rapid environmental degradation, which seems to drive pathogen eco‐epidemiological dynamics and disease emergence/re‐emergence (Rogalski et al., 2017). Yet, evidence exists that higher diversity increases pathogen prevalence and disease emergence (Dunn et al., 2010; Loaiza et al., 2019), and complex, null, or idiosyncratic relationships between biodiversity and pathogen dynamics have also been observed (Faust et al., 2017; Salkeld et al., 2013). Critiques to the dilution effect regard its generality, the definition of disease risk, spatial scale, and the lack of a mechanistic framework specifically showing the interactions that reduce disease (Strauss et al., 2015).

Ostfeld and Keesing (2000) were among the first to describe the dilution effect, presenting a conceptual model of how Lyme disease risk for humans, in North America, was reduced in regions with higher host species richness and evenness. Lyme disease is caused by Borrelia burgdorferi s.l. and is transmitted by an ixodid tick bite to a susceptible host, such as the white‐footed mouse (Peromyscus leucopus), the main reservoir of the Ostfeld and Keesing (2000) system. The dilution effect has been empirically observed in more diverse communities, with the underlining assumption that the order in which species assemble is deterministic and positively correlates with their susceptibility, so that more diverse communities will support a relatively larger number of rare species representing low‐competence hosts (LoGiudice et al., 2003; Ostfeld & Keesing, 2000; Schmidt & Ostfeld, 2001). It has been observed in an assembly of nested communities that changes in infection were determined from shifts in host species richness (Lacroix et al., 2014; LoGiudice et al., 2003; Ostfeld & LoGiudice, 2003). Another aspect investigated was the relationship between community richness and the identities or functional traits of its species, with richer communities assumed to have higher functional complexity, and included species that reduced pathogen transmission rates to vectors, or that reduced the abundance of amplification hosts (Johnson & Thieltges, 2010; Pongsiri et al., 2009). This led to the hypothesis that pathogen transmission is likely to rise in disturbed and depleted, communities, especially subjected to anthropogenic stressors (Johnson et al., 2019). However, relationships between wildlife diseases and human‐modified landscapes are rarely consistent, and might depend on pathogen transmission types and landscape scale and features (Brearley et al., 2013).

Nevertheless, the mechanisms underlining dilution were observed to be not directly related to basic biodiversity metrics such as richness or diversity (Johnson, Ostfeld, et al., 2015). For example, in their seminal work in the Lyme disease context, LoGiudice et al. (2008) found that species richness was only weakly negatively correlated with nymphal infection prevalence (NIP), while no relationship existed between the Shannon diversity index and NIP. The relationship biodiversity–pathogen‐prevalence may realistically rely on the identities and frequencies of host species (Strauss et al., 2015). Namely, the dilution effect is connected to how richness scales with relative densities of focal host and non‐host species (Strauss et al., 2015), as well as how nested is the community assembly, and the ecology of pathogen and parasites involved (Wood, Lafferty, et al., 2014). In particular, the direction of the aforementioned relationship depends on the overall community competence (i.e., the sum of each host species’ competence multiplied by its abundance), and how this community competence changes with diversity. If the reservoir hosts abundance increases linearly with richness (additivity), and therefore increases community competence, then high‐diversity systems may support more infection (amplification effect) (Strauss et al., 2015). Alternatively, more diverse communities might include species that reduce reservoir host population sizes (substitution) due, for example, to competition or predation (Wojdak et al., 2014).

Experimental and modeling studies focusing on the dilution effect have traditionally assumed nested communities, which are often observed in nature (Johnson et al., 2012; Wright et al., 1997), and additivity (e.g., Johnson et al., 2013; Roche et al., 2013), although a realistic community might be assembled following a tradeoff between additive and substitutive extremes (Becker et al., 2014; Mihaljevic et al., 2014). It is not likely that total host density increases linearly (i.e., additively) with increasing species richness because the species added to the community might be competitors or predators of reservoir hosts (Levi, Keesing, et al., 2016). So, understanding how the diverse community species abundance differ from less diverse communities is critical in determining how parasite and pathogen transmission is altered when biodiversity is lost (Wojdak et al., 2014).

In vector‐borne pathosystems, such as Lyme disease, it has been thought that more richness has an additive effect on overall tick host density and a reduction in tick infection prevalence as transmission is frequency dependent (more wasted bites on non‐competent hosts; Dobson, 2004; Rudolf & Antonovics, 2005). Yet, adding vector host species could increase overall tick abundance (vector amplification) providing more feeding opportunities to ticks and so a larger number of infected ticks (Randolph & Dobson, 2012; Wood & Lafferty, 2013). In directly transmitted pathogens, such as Hantaviruses, infection probability in the rodent reservoir host decreased when density of sympatric non‐reservoir hosts increased (Khalil et al., 2016; Suzán et al., 2009). In Sweden, the decline of Tengmalm's owls, a main predator of bank voles, which are the PUUV hantavirus reservoir host, has been observed to have contributed to higher PUUV infection rates in that species (Khalil et al., 2016). Likewise, in Africa, the decline of large wildlife increased the landscape‐level prevalence of rodent‐borne diseases (Young et al., 2014). In a mesocosm experiment, diversity of predators was the best negative predictor of trematode infections in frogs (Rohr et al., 2015). In a similar pathosystem, Wojdak et al. (2014) manipulated the density of a focal host (green frog tadpoles), while manipulating the diversity of alternative species, so that alternative species either replaced focal host species thus maintaining the total number of individuals at a constant level (substitution) or added to total host density (addition). They found that parasite transmission remained nearly equal (or decreased slightly) in the first case, while it was higher in case of additivity.

From these examples, it seems that multiple factors, particularly the density of focal host species along the biodiversity gradient and their ecological interactions in the community, regulate the relationship between biodiversity and pathogen transmission, leading to either dilution or amplification effects. Ecological interactions among species might determine the dilution mechanisms originally described by Keesing et al. (2006). These include encounter reduction between reservoir hosts or between hosts and vectors; alteration of host species abundance, distribution, behavior, and overall population dynamics; and competition for resources or space (Civitello & Hartman, 2021; Strauss et al., 2015). Diverse host communities may thus provide two types of dilution hosts: species that are abundant and/or heavily parasitized deflecting contacts between hosts and/or vectors; species that can reduce the abundance of the most competent hosts through interactions such as competition and predation (Levi, Keesing, et al., 2016). However, in vector‐borne pathosystems, for instance those involving generalist ticks, even assuming no additivity, it is plausible that some hosts amplify pathogen transmission because the proportion of “wasted” bites do not compensate the increase in overall tick population due to the potential host species variability in sustaining tick population at different life stages (McCoy et al., 2013).

Hence, another aspect to be considered in the dilution effect debate is the degree of the pathogen/parasite/vector generalism (as opposed as host specialization). Theoretical studies have suggested that the dilution effect will be greatest for pathogens in which within‐species transmission is greater than between‐species transmission and for which transmission is frequency dependent rather than density dependent (e.g., Dobson, 2004; Keesing et al., 2006). Transmission dilution is not expected to occur in the case of highly specialized pathogens and vectors, or when a single host does not allow maintenance of the pathogen cycle, while highly specialized pathogens may facilitate dilution when vectors have a higher degree of generalism (Ostfeld & Keesing, 2012). A recent review suggested that vector‐borne generalist pathogens, especially if transmitted by generalist vectors, were those most likely to be affected by changes to biodiversity (Rohr et al., 2020). Specific studies on ixodid tick species host generalism suggested that generalist patterns certainly vary by tick species but might also be determined by environmental conditions and the scale under consideration, and thus require explicit consideration in epidemiological models of tick‐borne disease (McCoy et al., 2013). However, no mechanistic frameworks integrating realistic community (dis)assembly and the subject of vector degree of generalism exist that are able to investigate the dilution effect mechanisms in a pathosystem, not even in the Lyme disease case, which is among the most studied in the context of the dilution effect.

To consider less extreme and more realistic assembly patterns (Johnson, Ostfeld, et al., 2015), we therefore investigated through a mechanistic eco‐epidemiological approach two tick‐borne pathosystems with different degrees of generalism, and with a community assembly that expanded from purely additive to purely substitutive community structures. We chose two rodent–tick–pathogen systems as wild rodents are ideal wild host communities for studying multi‐host pathogen transmission in the eco‐epidemiological framework, and tick‐borne pathosystems are widely represented in literature allowing parameter estimation and comparison with previous studies. Rodents are abundant and resilient to anthropogenic disturbance (Han et al., 2015; Mendoza et al., 2020). They are reservoirs of numerous zoonotic pathogens of public health interest, such as Lyme disease, Hantavirus diseases, plague, leishmaniasis, various hemorrhagic fevers, human anaplasmosis, and babesiosis (Meerburg et al., 2009). They live in sympatry, sharing pathogens and parasites (e.g., Begon et al., 1999; Paziewska et al., 2012), and displaying a positive correlation between resilience and host reservoir competence (Previtali et al., 2012). Understanding eco‐epidemiological dynamics of rodent‐associated pathogens is crucial from an applied pathogen‐control perspective, as they are predicted to rise as a consequence of large mammal defaunation (Young et al., 2014).

In particular, the two tick‐borne pathosystems chosen were rodent–Ixodes ricinus–Borrelia burgdorferi s.l. and rodent–I. trianguliceps–Babesia microti. Both tick species have a three‐stage life cycle (larva, nymph, and adult). I. ricinus (sheep tick) is a generalist tick, of which larvae and nymphs feed on small mammals while adults mostly prefer larger‐sized animals (Mysterud et al., 2015), and it is considered a key vector for B. burgdorferi s.l. (the causative agent of Lyme disease in humans; Paziewska et al., 2010). I. trianguliceps is a specialist tick for small mammals (Mysterud et al., 2015), which, in the UK, has been found to be the key vector for B. microti, a potentially zoonotic protozoan, of which voles seem to be the main reservoir (Bown et al., 2008; Hussein, 1980; Siński et al., 2006). These two pathogens have also different infectious periods (Harrison et al., 2011; Hartemink et al., 2008; Randolph, 1995; Randolph et al., 1996). I. ricinus is sympatric with I. trianguliceps in some areas, but the role of the former in Babesia sp. transmission, and the zoonotic potential of the latter are still unclear (Bown et al., 2006; Kovalevskii et al., 2013).

We integrated unique empirical data on wild rodent communities and mechanistic eco‐epidemiological modeling to improve our understanding of the functional form of density dependence for the tick population, particularly with regard to the density of infectious nymphs (DIN), and how the presence of predators and competitors influences the density of amplification hosts. Ground‐dwelling rodent species were surveyed in Wales (UK) in order to collect ad hoc empirical data to parametrize a compartmental model including epidemiological and population dynamics of host and vector, as well as ecological relationships such as competition and predation, with species able to affect host and vector populations but incapable of maintaining pathogen transmission. Simulations of realistic, nonrandom community assemblages were performed to investigate how the interactions between host and non‐host species affected pathogen dynamics and the mechanisms behind the potential dilution or amplification effect. Our simulations aimed to identify (1) whether the dilution and amplification effects occurred under the assumption of partially additive host densities and which were the mechanisms involved; (2) whether these effects occurred in both pathosystems and whether there were differences due to the different degree of vector generalism; and (3) whether there were identifiable specific characteristics in each pathosystem determining dilution or amplification to provide practical disease risk management recommendations (i.e., if there were crucial information on the community or the vector in a specific area to be collected to make realistic predictions on potential disease risk). These were particularly relevant to the study area, where rodent tick‐borne zoonotic risk was not previously assessed despite the close interface between human‐livestock‐wildlife. The additional knowledge provided on tick‐borne pathosystem dynamics fits into a growing corpus of information required to improve wildlife management and zoonotic risk control policies in these unprecedented times of biodiversity decline.

METHODS

Ground‐dwelling rodent populations survey

Rodents were live trapped twice a year (spring and autumn) in seven different sites in Wales between June 2015 and June 2017, including a range of different habitats (semi‐deciduous woodland, conifer plantation, grassland, clear‐cut/scrubland). At each location, the square grid consisted of 36 trapping stations (6 × 6) with one Longworth and one Sherman trap at each station to maximize estimation of species composition (Anthony & Ribic, 2005). The traps were set with appropriate bedding material and food; each trapping occasion consisted of four consecutive days and nights, with the first day/night being pre‐baiting, and the traps were checked twice a day (early morning and sunset). The live‐trapping was performed during two different seasons to estimate individual densities of the pre‐breeding recruitment population (May–June) and the post‐breeding peak population (September–October). Each individual, at the time of first capture, was identified at the species level, sexed, assigned to an age class according to Telfer et al. (2002), weighed, individually marked by fur clipping, and finally released. On first capture only, ticks were collected from each individual, and stored for further molecular investigations (Occhibove, 2018).

Individual density was estimated by species, site, and trapping season using the POPAN algorithm (Schwarz & Arnason, 1996) within the software MARK (White & Burnham, 1999), assuming an open population, constant survival, and constant capture probability. Goodness‐of‐fit was tested with the RELEASE suite within the software. Growth rate, r, was estimated for each species according to the formula proposed by Lambin et al. (2000): r  = logN  − logN where N is the population density of species i at time t. This process generated distinct growth rates representing breeding and non‐breeding seasons, March to October and November to February respectively according to Hörnfeldt (1994) and Stenseth et al. (2002). Unless specified differently, all data analyses were performed in R (R Core Team, 2016).

Model structure and simulations

Similarly to Arino et al. (2004) and Venturino (1994, 2001, 2002), we developed a demographic model that accounted for interactions between different populations in which the pathogen spreads among host species but does not affect their survival rates. Indeed, B. burgdorferi s.l. and B. microti have not been observed to impact rodent mortality (Bown et al., 2008; Randolph et al., 1996; Telfer et al., 2008; Voordouw et al., 2015). Our modeling approach consisted of a deterministic multi‐host single‐vector SIR (Susceptible‐Infectious‐Recovered)‐SI (Susceptible‐Infectious) compartmental model made up of a set of differential equations (Appendix S1: Equations S1–S18 and model code in Data S1: Tick_borne_model.R). Figure 1 illustrates the structure of the model (Figure 1a) and the tick life cycle (Figure 1b). All the results were produced using the function rk4 in the R package deSolve (R Core Team, 2016), which is based on the classical Runge–Kutta fourth‐order integration. The model was not explicitly spatial. The area, 1 ha, was considered constant since it was constrained by the experimental sampling unit, therefore the populations were expressed in individuals/ha and the relevant parameters were scaled accordingly. Dilution was tested by modeling a progressively more complex community: species were added in turn following nonrandom assembly rules. The structure of the model was chosen to offer a compromise between complexity and parsimony, and it was parametrized, where possible, with data collected during live trapping.

FIGURE 1

(a) Tick‐borne disease compartmental model and (b) tick life cycle. Boxes represent epidemiological compartments in which each population is subdivided: L, tick larval stage; S, susceptible; I, infectious; R, recovered. Subscripts are n, tick nymphal stage; a, tick adult stage; i, rodent species (host). Vectors can feed also on non‐host competitor species (shrew) and specialist predators. Arrows indicate the direction of movement of individuals between classes (solid line, vector; dash‐dotted line, host). Arrows pointing outside the boxes represent mortality (dashed lines, host mortality through predation; dotted lines, vector natural mortality). Thicker arrows represent transmission. Each of the five scenarios in Table 1 were based on this model, with non‐host competitor occurring in scenarios 3 to 5, generalist and specialist predation occurring only in scenarios 4 and 5

TABLE 1

Species composition for each community assemblage scenario

Name	Community composition a	Species role in the community			Vector life stage preferentially hosted
					I. ricinus (generalist tick)	I. trianguliceps (specialist tick)
Scenario 1 (single host)	Bank vole	Pathogen host	Vector host	Rodent competitor	Larva, nymph	Larva, nymph, adult
Scenario 2 (two hosts)	Wood mouse	Pathogen host	Vector host	Rodent competitor	Larva, nymph	Larva, nymph, adult
Scenario 3 (rodents + competitor)	Sorex sp. shrew		Vector host	Rodent competitor	Larva, nymph	Larva, nymph, adult
Scenario 4 (Rodents + competitor + predation)	Specialist predator (mustelid)		Vector host	Predator of all species	Nymph, adult	…
	Generalist predator			Predator of all species	…	…
Scenario 5 (full community)	Ungulate (e.g., sheep or deer)		Vector host		Adult	…

Note: The communities are nested in that each sequential scenario includes the species of previous scenarios in addition to the additional species listed. Each of these five scenarios was simulated for each of the two pathosystems (Ixodes ricinus–Babesia burgdorferi and Ixodes trianguliceps–Babesia microti).

Each scenario incorporates each of the species in rows above in addition to those in the scenario row.

Five scenarios (Table 1) for the community assembly were simulated based around the structure outlined in Figure 1, with each scenario representing an increase in complexity from the previous one:Each of these five scenarios were simulated for both pathosystems (I. ricinus‐B. burgdorferi and I. trianguliceps‐B. microti) giving 10 different modeled scenarios in total. Although Scenario 5 in the second pathosystem was identical to Scenario 4 due to the tick specialization on small hosts and the null effect of ungulates on other species population dynamics.

Scenarios 1 and 2: as was found in field surveys, wood mice and bank voles were rodent host species for both the pathogen and the tick vectors. Scenario 1 involved one host species (bank vole), whereas Scenario 2 involved both wood mice and bank voles.

Scenario 3: Sorex ssp. shrews were added as a sympatric competitor species not hosting the pathogen. Evidence of their presence (shrew feces found in traps equipped with shrew escape holes) were found in the sampling sites, and they were suitable hosts for ticks.

Scenario 4: addition of two types of predators, specialists and generalists. Specialist predators may be represented by mustelids and are also suitable tick hosts. Generalist predators may be represented by birds and do not affect tick population.

Scenario 5: addition of a large grazing ungulate, i.e., sheep or deer, which is a suitable host for adult stages of I. ricinus but incapable of pathogen transmission. Otherwise ungulates are assumed to be a constant population not affecting, or affected by, any population in the system.

Rodent and shrew populations were modeled according to the Lotka‐Volterra system, namely a logistic growth tending to a species‐specific carrying capacity and limited by intra‐specific density dependent reduction and interspecific competition (Hanski et al., 1993; Lotka, 1925), as these are considered crucial population dynamics in rodents (Burthe et al., 2010; Occhibove, 2018; Smith et al., 2006). Interspecific competition among rodents, and between rodents and shrews was modeled with a density‐dependent competition term (Huitu et al., 2004; O'Regan et al., 2015), considering rodent species better competitors than shrews (Henttonen et al., 1989; Huitu et al., 2004). These species were all predated upon, but rodent species were considered preferential prey (Korpimäki, 1992; Korpimäki & Norrdahl, 1989). Generalist predation was modeled according to the alternative prey hypothesis (Holling type III functional response), while specialist predator population was modeled according to the Holling type II functional response based on the Rosenzweig‐MacArthur model (1963) with no preference among rodents (Erlinge, 1975; Hanski & Henttonen, 1996; Krebs & Myers, 1974). Specialist predation parameters were based on Mustela nivalis (least weasel; Table 2), which was widespread across sampled sites. Tick population was modeled with a density‐dependent fecundity reduction term (Norman et al., 1999; Ogden et al., 2007) and aggregation, i.e., non‐homogeneous distribution of vectors among the host population (Rosà et al., 2003). The aggregation term was parametrized according to the results of tick collection during the rodent survey (reported in Parameter estimation; Table 2).

TABLE 2

State variable initial values and parameter values (range of values included in sensitivity analysis are included in brackets)

Symbol	Description	Value	Source
N w	Wood mouse population (initial)	49	This study
N b	Bank vole population (initial)	75	This study
N j	Shrew population (initial)	20	This study
N p	Weasel population (initial)	3	This study
N d	Ungulate population (throughout)	5	This study
N v	Tick population (initial)	100	This study
cbw, cjw	Competition of wood mouse over bank vole, and shrew respectively	0.20 1.04	This study based on O'Regan et al. (2015)
cwb, cjb	Competition of bank vole over wood mouse, and shrew respectively	0.20 1.03	This study based on O'Regan et al. (2015)
cwj, cbj	Competition of shrew over wood mouse, and bank vole respectively	0.11 0.12	This study based on O'Regan et al. (2015)
d r d s d l	Molting‐feeding success on Rodents Shrews Larger hosts (weasel/ungulate)	0.14 (0.1–0.59) 0.16 (0.1–0.50) 0.21 (0.1–0.64)	LoGiudice et al. (2003)
G	Saturation rate of generalist predation	0.49	This study based on Turchin & Hanski (1997)
H	Prey density at which generalist predation rate is half of the maximum	9.9	This study based on Turchin & Hanski (1997)
K	Tick aggregation parameter	0.18	This study based on Rosà et al. (2003)
numegg	Maximum number of per capita adult female tick eggs production Ixodes ricinus I. trianguliceps	1500 1000	I. ricinus: Norman et al. (1999, I. trianguliceps: Krasnov et al. (2007)
q	Specialist predator–prey ratio constant	56	This study based on Turchin & Hanski (1997)
r b+ , r b−	Bank vole growth rate breeding season (+), and non‐breeding season (−)	0.007; −0.002	this study
r w+ , r w−	Wood mouse growth rate breeding season (+); non‐breeding season (−)	0.04; −0.006	this study
s v	Density‐dependent reduction of tick growth rate	‐	Formula in Ogden et al. (2007)
Α	Maximum rodent consumption rate of specialist predator	1	This study based on Turchin & Hanski (1997)
αs	Maximum shrew consumption rate of specialist predator	7.67	This study based on Turchin & Hanski (1997)
βsl βsn βsa	Encounter rate small host/larva I. ricinus; I. trianguliceps Encounter rate small host/nymph I. ricinus; I. trianguliceps Encounter rate small host/adult I. ricinus; I. trianguliceps	0.040; 0.040 0.040; 0.040 0.020; 0.040	Dobson et al. (2011), Hancock et al. (2011)
βll βln βla	Encounter rate large host/larva I. ricinus; I. trianguliceps Encounter rate large host/nymph I. ricinus; I. trianguliceps Encounter rate large host/adult I. ricinus; I. trianguliceps	0.025; 0.000 0.040; 0.000 0.060; 0.000	Dobson et al. (2011), Hancock et al. (2011)
Δ	Half‐saturation constant (rodent)	11.31	This study based on Turchin & Hanski (1997)
Δs	Half‐saturation constant (shrew)	22.62	This study based on Turchin & Hanski (1997)
νj	Shrew birth rate	2.60	This study based on de Leo & Dobson (1996), Bolzoni et al. (2008)
νp	Weasel birth rate	3.23	This study based on de Leo & Dobson (1996), Bolzoni et al. (2008)
ρj	Shrew death rate	1.04	This study based on de Leo & Dobson (1996), Bolzoni et al. (2008)
ρp	Weasel death rate	1.29	This study based on de Leo & Dobson, 1996, Bolzoni et al., 2008)
ρv	Tick death rate: larva nymphs adults	0.0014 0.0005 0.0004	Dobson et al. (2011)
σbb	Recovery rate Borrelia burgdorferi s.l.	0.0083	Harrison et al. (2011), Hartemink et al. (2008), Randolph (1995), Randolph et al. (1996)
σbm	Recovery rate Babesia microti	0.4	Harrison et al. (2011), Hartemink et al. (2008), Randolph (1995), Randolph et al. (1996)
τ	Reservoir competence of transmission: Host to vector Vector to host	0.50 (0.1–1) 0.80 (0.1–1)	Giardina et al. (2000), Harrison et al. (2011), Hartemink et al. (2008), LoGiudice et al. (2003)

Note: Rates are expressed per day in accord with the model time step.

A simulation time period of 20 years, with a daily time step, allowed equilibration. The model included susceptible, infectious, and recovered compartments (SIR) for each competent mammalian host; and susceptible and infectious (SI) for the vector (Porco, 1999). Non‐systemic transmission through tick co‐feeding was not considered because it was found to be a very minor, or inefficient route of transmission for the chosen pathogen (Jacquet et al., 2016). Both pathogens under consideration do not present vertical transmission, consequently larvae were not infectious, but could be infectious and molt into infectious nymphs (Gray, 2006; Randolph, 1995). Reservoir competence values summarized susceptibility, ability of the pathogen to magnify and persist in the host/vector and efficiency of transmission; all of the individuals of the same host species were equally competent, but different species might display different levels of competence. Initial density for each rodent species was set at the average value found in the field (bank vole, 75 ± 68 individuals/ha; wood mouse, 49 ± 53 individuals/ha). The inoculum was represented by a single infectious individual. Appendix S1: Equations S1–S7 represent the model of competent, non‐competent host species and vector dynamics in the absence of pathogen, while Appendix S1: Equations S8–S18 represent the SIR‐SI dynamics when the pathogen is introduced. Model variables and parameters are listed in Table 2 together with initial conditions for simulations and parameter values. Model code provided in Data S1.

The model outputs of interest were number of infectious hosts and number of infectious nymphs, named hereafter density of infectious nymphs (DIN; from Giardina et al., 2000). These were chosen because they are commonly used to estimate disease risk (LoGiudice et al., 2003; Piesman, 1989); in particular DIN is considered a better metric than nymphal prevalence to assess human zoonotic risk (LoGiudice et al., 2008). However, more relevance will be given to the actual numbers of infectious individuals as they provide the best interpretation of transmission dynamics without interference of effects from population dynamics across different communities (i.e., change in susceptible host population due to community assembly that altered the infectious proportion).

Model parameterization

Parameters were estimated, where possible, from rodent live‐trapping in this study, otherwise they were taken from relevant literature (Table 2 and Appendix S2). Rodent species growth rates (r) in the model were the averages for each species and for season (Appendix S3: Table S1) scaled by day. Seasonal growth rates (r) were used to simulate rodent intra‐annual fluctuations and the tick population diapause (Randolph, 2004).

Carrying capacity (K), and birth and death rates of non‐host species were estimated allometrically (Bolzoni et al., 2008; de Leo & Dobson, 1996; Appendix S2). The competition factor (c) was computed algebraically (Appendix S2). For the simulations, we used conservative values of competition (Table 2), according to Occhibove (2018). Predation parameters were estimated with methods in Turchin and Hanski (1997), separately for rodents and shrews, as the latter have a significantly lower body mass and are secondary prey items (Henttonen et al., 1989).

Parameters regarding ticks with relative source and value are listed in Table 2. Only the aggregation parameter (k) was calculated from empirical data (Appendix S2), while other parameters were taken from relevant literature allowing to model the differences in degree of generalism and ecology in general among the two tick species. Values for reservoir competence of hosts and vector, host–vector encounter rates, and molting success were taken from literature (Table 2). Recovery rates (σ) were pathogen specific and represented the reciprocal of the infectious period (Table 2).

Sensitivity analysis

By keeping constant the parameters estimated in this study and parameters characterizing the differences of the two pathosystems, a sensitivity analysis was performed on parameters that were not empirically estimated. The ranges for these parameters were selected from various studies on the most similar pathosystems, namely molting success on various hosts, and reservoir competence. These parameters were observed to greatly affect disease transmission and dilution output (Dunn et al., 2013; Roche et al., 2013), but uncertainty on the values existed for the pathosystems under consideration. We used Latin Hypercube Sampling (LHS) and partial rank correlation coefficients (PRCC) to assess the impact of uncertainty and the sensitivity of the outcomes of the simulations to variations in each parameter (Iman et al., 1981a, 1981b; Marino et al., 2008). To generate the LHS matrices, we assumed that all the selected parameters were uniformly distributed. Then a total of 1000 simulations of the models per scenario were carried out. The baseline conditions and the ranges used are given in Table 2; the ranges were chosen as minimum and maximum values found from empirical studies in literature for each parameter. Sensitivity analysis was performed separately for each host–tick–pathogen system, keeping the parameters not subject to sensitivity analysis constant. The response functions used were the density of infectious nymphs (DIN) and the number of infectious hosts at equilibrium.

RESULTS

During the 4968 trap‐nights of fieldwork, 258 bank voles (Myodes glareolus) and 230 wood mice (Apodemus sylvaticus) were captured, with signs of shrew visits in every trapping grid, i.e., shrew feces were recovered in trap tunnels. Sex ratios were balanced with a 1.08 males per female overall, while the ratio of wood mice to bank voles was 0.86. Inter‐seasonal fluctuations of population density were evident, and, as expected, higher density was found in the post‐breeding‐peak season (autumn). Non‐breeding season growth rates were mostly negative, or much lower than breeding season values in both species (Appendix S3: Table S1). General trends were similar between the two species. Density values across all field seasons were much lower than the estimated carrying capacity for each species, suggesting that in our sites the populations might be subjected to interspecific competition and predation. Hence, model structure reflected field observations; numerical simulations for rodent population sizes across different scenarios are shown in Appendix S3: Figure S1.

Ticks were found to be highly aggregated on hosts and the dispersion parameter of the negative binomial distribution fitting field data was used in the model to represent aggregation (k was 0.18, indicating a much more aggregated distribution that is frequently observed in ticks; see, e.g., Ferreri et al., 2014). In total, 225 ixodid ticks were collected from 120 rodents, 16.28% of total individuals sampled. Total infestation prevalence (i.e., including all occasions rodents were caught irrespective of previous trapping status) was 15.99%, with I. trianguliceps representing most of the sample set (196 specimens), followed by I. ricinus (10 specimens), while 19 specimens were not taxonomically identifiable.

Results of the sensitivity analysis showed that the parameters overall mostly influencing the outputs in both pathosystems were molting success on rodents, and transmission competence vector to rodent (Appendix S3: Figure S3). As these parameters made a substantial impact on the conclusions drawn from the main modeling results, results of model simulations were displayed with respect to these two influential parameters.

Model simulations

Firstly, host–vector interactions were investigated through simulations in the absence of pathogens to display the differences in dynamics of the two tick species. Literature indicated that fecundity of I. ricinus, is great er than I. trianguliceps. The inclusion of this difference in the no‐pathogen model resulted in a higher number of individuals across all scenarios for the generalist tick I. ricinus when compared to I. trianguliceps; more so when larger hosts were added to the community (Scenario 4/5; Appendix S3: Figure S2). Conversely, I. trianguliceps only feeds on small hosts (rodents and shrews) therefore there was no substantial change in tick numbers across different communities, apart from Scenario 3, when the shrew population was added, determining an overall host increase and consequent vector increase (Appendix S3: Figure S2). Tick populations were subjected to a density dependent fecundity reduction, therefore did not linearly increase with host population.

As our sensitivity analysis showed that molting success on rodents (henceforth referred to as simply molting success), and the competence of transmission from vector to rodent affected model output the most (see previous section), Figures 2 and 3 display model output variation across scenarios over the range selected for those parameters (i.e., the percentage change in the output of interest taken at equilibrium). In the I. trianguliceps–B. microti system, Scenario 5 was identical to Scenario 4 due to the specialization of the vector for small hosts and the absence of population interactions between ungulates and other hosts.

FIGURE 2

Percentage of change from Scenario 5 in (a, b) density of infectious nymphs (DIN) and (c, d) infectious hosts across the five scenarios (on the x‐axis) and the range of molting success on rodents (d ) (colored circles) in the two pathosystems: Ixodes ricinus–Babesia burgdorferi (top) and I. trianguliceps–B. microti (bottom). Black dashed line, no change. Box and whisker plots on top of colored circles represent the percentage of change distribution across the range of d values. Lower and upper box boundaries are 25th and 75th percentiles, respectively, line inside box is the median, lower and upper error lines are 10th and 90th percentiles, respectively, filled circles show data falling outside 10th and 90th percentiles. Decrease of percentage of change from less complex to more complex scenario represents dilution, while increase represent amplification. Other parameter values are listed in Table 2

FIGURE 3

Percentage of change from Scenario 5 in (a, b) DIN and (c, d) infectious hosts across the five scenarios (on the x‐axis) and the range of transmission competence of vector to rodents (τ, colored circles) in the two pathosystems: I. ricinus–B. burgdorferi (top) and I. trianguliceps–B. microti (bottom). Black dashed line, no change. Box and whisker plots on top of colored circles represent the percentage of change distribution across the range of τ values. Lower and upper box boundaries are 25th and 75th percentiles, respectively, line inside box is the median, lower and upper error lines are 10th and 90th percentiles, respectively, filled circles show data falling outside 10th and 90th percentiles. Decrease of percentage of change from less complex to more complex scenario represents dilution, while increase represent amplification. Other parameter values are listed in Table 2

As expected, in the I. ricinus–B. burgdorferi system, DIN increased together with the increase of the molting success, while the opposite trend could be observed for prevalence because of the rise of total tick population (not shown). Single‐host scenario (Scenario 1) and the most complex communities (Scenario 4 and 5) gave virtually identical results (at equilibrium), with the only difference being a slightly lower DIN in the first case at low molting success values (Figure 2a). On the other hand, Scenario 2 and 3 resulted in a lower DIN compared to the other communities across the entire range of parameter values and displayed very similar values. The DIN difference between the two groups of scenarios (Scenario 1, 4 and 5 vs. Scenario 2 and 3) decreased while increasing the parameter value; yet, at the highest end of the range, DIN was still higher ~25% in the simplest community (large population of competent host) and in the communities where tick population was amplified by a larger number of hosts (Figure 2a). This means that we could observe a dilution effect with the introduction of competitor species, but this was outweighed in Scenario 4 and 5 by the introduction of other host species for ticks. Changing the value of molting success only made very slight differences to the total number of infectious hosts. However, Scenario 2 and 3 (which gave almost identical results in this case) presented the highest number of infectious hosts; while the increase in community complexity in Scenario 4 and 5 resulted in a reduction of transmission at levels equal to or even lower than Scenario 1 (at low molting success values), producing dilution of transmission (Figure 2b).

Similarly, the increase of molting success determined an increase in DIN and infectious hosts in the second pathosystem, I. trianguliceps‐B. microti, although there were major differences between the two systems. DIN increase over the range of the parameter was significant, with all scenarios reaching similar numbers of infectious nymphs at the highest end of the molting success value (Figure 2c). In this case, scenarios are roughly ordered, in terms of DIN, from the least complex having lower values of DIN to the most complex. Namely, the increase in small‐sized hosts, for which this tick is a specialist, determined a DIN increase. Yet, the introduction of predation in Scenario 4, which determined a significant reduction of competent hosts but not shrews, resulted in slightly lower DIN. This is because the tick population was sustained by the shrew population, but it was slightly less likely for tick individuals to get infected on the less infectious hosts (Figure 2b). Therefore, only little dilution effect could be observed. Nonetheless, dilution could be identified in hosts, for any value of molting success, when predation was introduced (Scenario 4) with respect to Scenario 3. This occurred only above certain molting success values in the case of the first two scenarios (Figure 2d).

When varying the competence of transmission from vector to rodent parameter a very similar pattern compared to the previous case could be observed in the I. ricinus–B. burgdorferi system. The variation in the parameter had relatively little effect on the DIN in the same scenario. However, it decreased from Scenario 1 to 2, then did not significantly change with the introduction of shrews (Scenario 3), while it increased due to the introduction of predators and ungulates (Scenario 4 and 5; Figure 3a). This might have been the result of an interaction between the relaxation of the fecundity reduction with consequent tick release (less small mammals in Scenario 4 and 5), and the modeled aggregation on competent hosts (higher probability of finding an infectious host). Conversely, infectious hosts increased when the second host species was added to the community (Scenario 2), remained stable when shrew population was added (excluding a very low competence value; Scenario 3), while decreasing in the most complex communities simulated (Scenario 4 and 5), which presented the same results (Figure 3b). Nevertheless, at very low competence values, the single‐host community (Scenario 1) presented the least infectious hosts. DIN values and the number of infectious hosts did not increase linearly with competence as the tick population was regulated by a fecundity reduction function and was unevenly distributed in the host population (see Results: Ground‐dwelling rodent populations survey) and therefore no infinite increase of infections was possible. In addition, I. ricinus reproduced abundantly and being generalist had an abundance of hosts in every scenario, so low competence values were sufficient to sustain transmission.

The competence of transmission from vector to rodent appeared to be crucial in determining the final number of infectious hosts in the second pathosystem (I. trianguliceps–B. microti) when only rodents were in the community (Scenario 1 and 2; Figure 3d), while relatively no significant change could be observed across the parameter range for DIN in the same scenario (Figure 3c). Namely, a small level of competence was enough to saturate transmission in nymphs due to tick specialization, while an increase in competence of transmission from vector to rodent increased the number of infectious hosts in Scenarios 1–3. In this pathosystem, no dilution was observed when DIN was taken as metric of interest; on the contrary, the more hosts and the more complex the community, the more transmission seemed to be amplified (Figure 3c). However, shrews (Scenario 3) and predation (Scenario 4) determined a reduction of infectious hosts compared to the only rodents scenarios (Scenarios 1 and 2; Figure 3d).

The percentage of change at equilibria (Figures 2 and 3) suggested that, in both pathosystems, the dilution effect occurred in the competent host populations mostly due to the increase in functional complexity of the community (due to the addition of non‐rodent species from scenario 3 to 5), not by the increase of number of species per se. To appreciate differences over the course of epidemic, Figure 4 shows DIN and number of infectious hosts for each scenario and each pathosystem across the 20 years of simulations (parameters in Table 2).

FIGURE 4

Temporal differences of (a, b) DIN and (c, d) total number of infectious hosts in the two pathosystems over the entire simulation time (20 years). (e) Number of infectious hosts generated over the entire simulation time (20 years) by the specialist tick I. trianguliceps transmitting B. burgdorferi (modification of the second pathosystem). Black, scenario 1; purple, scenario 2; teal, scenario 3; light green, scenario 4; yellow, scenario 5. Parameter values are listed in Table 2

In the first pathosystem, no dilution could be observed regarding DIN (Figure 4a). Nonetheless, at the start of epidemic, Scenario 1 presents a much higher DIN compared to the other ones, while converging to Scenario 5 at the end. Although the addition of shrews determines a rapid DIN increase compared to Scenario 2, after a few years, despite the overall greater host population, the transmission is slowed down by “wasted bites” at the level of Scenario 2. On the contrary, in the I. trianguliceps–B. microti system, the addition of more hosts simply increases DIN (although the lower fecundity assumed determines overall smaller tick population and lower DIN; Figure 4b). Looking at the competent host population, no significant dilution occurred in the first pathosystem unless predators were added (Scenario 4), and the addition of ungulates (Scenario 5) further reduced the number of infectious hosts as they were dead ends in terms of transmission (Figure 4c). Nevertheless, comparing Scenario 2 and 3, it could be appreciated how the remarkable increase in total tick host population due to shrews did not correspond in such increase in infectious hosts, suggesting that the transmission reduction mechanism was an important source of infection dilution (Figure 4c). The I. trianguliceps–B. microti system presented a significantly lower number of infectious hosts over the entire course of epidemic no matter the scenario compared to the first system (Figure 4d). The epidemic is characterized by a lower peak compared to the first pathosystem and at equilibria only few individual rodents became infectious (Figure 4d). This was due to the lower tick fecundity, the higher tick specialization and most of all due to the much shorter infectious period of B. microti. In this case, the addition of a second rodent species (Scenario 2) markedly increased the number of infectious hosts compared to Scenario 1, while the addition of competitors (Scenario 3) and predators (Scenario 4) lowered the epidemic peak compared to Scenario 2 and decreased infectious hosts below the on rodent‐only scenarios after several years (Figure 4d).

To investigate the effect of pathogen infectious period on dilution, and disentangle its effect from tick ecology, simulations were performed assuming I. trianguliceps transmitting B. burgdorferi. The number of infectious hosts are shown in Figure 4e. Comparing these results with those obtained when we simulated I. ricinus transmitting the same pathogen (Figure 4c), we could observe a different course of epidemic, more similar to that in Figure 4c, i.e., a much higher number of infectious hosts compared to Figure 4d (but lower than Figure 4c), and overall amplification when increasing host population available to the tick. Only predation (Scenario 4) decreased infectious hosts compared to the previous scenario, but not at lower levels than Scenarios 1 and 2.

Thus, it seemed that a sort of hump shaped relationship between community complexity and transmission could be observed for the pathosystem including the generalist vector (I. ricinus). Namely, in our simulated scenarios, the single‐host one (Scenario 1) was comparable to the most complex ones (Scenarios 4 and 5), while Scenarios 2 and 3 displayed the peak of the hump (either negative in the case of DIN or positive in the case of infectious hosts). This occurred because ecological relationships (which altered host densities and consequently tick populations) and tick generalism both diluted infections in the most complex communities at the level of the single‐host scenario. This was performed through the mechanisms of transmission reduction (“wasted bites”) and susceptible host regulation (Figure 5). The addition of shrews in Scenario 3 did not impact results as one would expect, considering the significant increase in host population compared to Scenario 2. This is because the high amount of transmission reduction and the density dependent fecundity reduction in ticks did not allow an infinite tick population growth. In Scenarios 4 and 5, amplification and dilution effects occurred simultaneously in vector and host populations respectively.

FIGURE 5

Conceptual model of the dilution mechanisms described by the model. (a) Transmission reduction: addition of non‐competent hosts to the community (scenario 3–5) leads to a reduction in numbers of infectious hosts because of wasted bites. Infectious ticks, in black, feed on non‐competent hosts halting pathogen transmission to hosts (and molting in non‐infectious adults), especially when ticks do not display preferences between competent and non‐competent hosts (e.g., pathosystem with generalist vector). Dashed line separates competent hosts, on the left, from non‐competent hosts, on the right. (b) Susceptible host regulation: increase in community diversity (scenario 1–5) leads to competent host population reduction and so to a reduction of infectious hosts; however, as the number of available hosts for the ticks also rises, there might be a concurrent increase of infectious ticks, especially if ticks feed preferably on competent hosts (e.g., pathosystem with specialist vector)

In the pathosystem with the specialist tick, this relationship was more complex and more dependent upon parameter values. However, amplification occurred when considering DIN, with a progressive increase in infectious nymphs from Scenarios 1 to 4; while dilution could be observed when considering infectious hosts, and was more significant when predation was introduced (Scenario 4). Susceptible host regulation seemed the dilution mechanism more consistent for the pathosystem with the specialist vector. The effect of “wasted bites” observed in the pathosystem with a generalist vector was not crucial here, and determined dilution only when the pathogen infectious period was short.

In summary, there were significant differences between the pathosystems:

The host specialization of I. trianguliceps (as it fed only on small‐sized hosts) was a significant difference that led to a number of results: (a) a lack of dilution effect with regard to DIN in every scenario; (b) a potential amplification of infectious hosts in Scenario 3 when the shrews, able to sustain tick population, were added to the community.

The larger number of “wasted bites” of the generalist I. ricinus led to dilution in hosts in each scenario comprising alternative non‐competent hosts (Scenarios 3–5).

The generalism of I. ricinus (likely combined with the higher fecundity) determined that low molting success or low competence values were sufficient to reach the peak of infections.

DISCUSSION

The rodent species growth rates estimated were comparable with other live‐trapping studies (Careau et al., 2013; Huitu et al., 2004; Lambin et al., 2000; Merritt et al., 2001). Population densities suggested that field sites hosted rather diverse communities. For most of the sites, a species list including avian and terrestrial predators of our focal species was available, and the mid‐Wales sites were included in a pine marten reintroduction project; shrew feces were recovered in trap tunnels at every site (Occhibove, 2018). Interspecific competition among rodents might also have played a role in keeping population sizes under control, with bank voles seemed to be overall slightly better competitors than wood mice (Occhibove, 2018). This competitive interaction could be hypothesized to be due to the variability in intra and interspecific contacts among rodent species at different sites. This might be explained by different resource and space use, and by the inverse relationship between proportion of interspecific contacts and individual density (Occhibove, 2018).

The numbers of ticks recovered from rodents were in agreement with other studies in which collection was made on living, anesthetized individuals (e.g., Paziewska et al., 2010). The proportion of the population parasitized by ticks and fleas was very small, supporting the “20/80 Rule,” which means ticks were highly clustered on a relatively small proportion of the host population (Occhibove, 2018). In general, bank voles were more parasitized than wood mice, which have been observed to be more heavily parasitized by I. trianguliceps than murine species in England (Hussein, 1980). A low number of I. ricinus was found at the sampled sites, but this does not exclude the potential presence at nearby sites, as a very fine‐scale distribution of this species in various habitats was observed in the UK (Dobson et al., 2011). Indeed, tick distribution seems to be extremely random and patchy, being determined by local host–vector interaction dynamics (Lorusso et al., 2011); while the ability to recover them it is highly dependent on time and type of sampling, which might hide the real tick species composition (Dobson et al., 2011). Another explanation might be that at our field sites the most frequent ungulate species were sheep and cattle, while deer were absent, and the latter were found more positively associated with I. ricinus presence in the UK (Medlock et al., 2008). So, the higher prevalence in bank voles may be determined by the overall dominance of I. trianguliceps, since, when I. ricinus has been found to be the dominant tick species, wood mice were the most parasitized hosts (e.g., Gray et al., 1999; Kurtenbach et al., 1995).

This study integrated empirical data about local host communities with mechanistic models to develop a realistic eco‐epidemiological approach for the investigation of tick‐borne pathogens in pathosystems differing in vector degree of generalism. Our main result was the identification of a hump‐shaped relationship between community complexity and transmission in the pathosystem including the generalist vector (I. ricinus). This was due to the interaction of ecological relationships, which altered host densities and consequently tick populations, and tick generalism that determined dilution of infection through the mechanisms of transmission reduction (“wasted” bites) and susceptible host regulation (Figure 5). Conversely, in the pathosystem with the specialist tick, this relationship was more complex and seemed to be variable depending on parameter values. In both cases, dilution occurred when considering the host population but not the vectors.

Ixodes ricinus, which displays a higher degree of host generalism (Mysterud et al., 2015), seemed to realize a higher number of “wasted” bites, in terms of pathogen transmission, while I. trianguliceps largely fed on competent hosts. Nonetheless, a lower value of molting success on rodents (the competent hosts), compared to I. trianguliceps, was sufficient to saturate transmission in hosts due to the greater population. On the contrary, the I. trianguliceps smaller population determined that a higher value of molting success was necessary to saturate transmission in hosts. In this pathosystem, molting success on rodents appeared to be the predominant regulator of pathogen transmission. Whereas, competence of transmission was crucial in determining the epidemic size only in the scenarios where fewer competent hosts were available. In general, in both pathosystems, variability in competence of transmission vector to rodents resulted in less variability in DIN and infectious hosts compared to molting success hosts, suggesting that competence of focal hosts might not be a limiting factor to transmission.

Our model highlighted the differences in the two host–vector–pathogen systems, justifying the independent modeling of the two pathosystems and supporting the empirical observation of tick and pathogen ecological differences (Mysterud et al., 2015). Our simulations revealed vector population density and proportion between life stages to be comparable with other studies (Dobson et al., 2011; Harrison et al., 2011). The numbers of infectious individuals predicted by our model were plausible in both pathosystems, with B. burgdorferi s.l. being significantly more prevalent than B. microti (Hussein, 1980; Obiegala et al., 2017; Welc‐Falęciak et al., 2008).

The first dilution mechanism detected in our model was transmission reduction (Figure 5a). This was observed after the introduction of non‐competent hosts into the community, e.g., shrews, which provided a source of “wasted” tick bites (Halsey et al., 2018). This was more marked in the I. ricinus–B. burgdorferi system, as the generalist tick fed also on the larger hosts included in Scenarios 4 and 5. In this pathosystem, despite the remarkable rise of overall hosts available for ticks, there was no such increase in infectious hosts, although the actual number of infectious nymphs was amplified. In more diverse community assemblages, the greater host availability for the vector led to a rise in the overall vector population, including infectious vectors (amplification effect), but infectious hosts decreased again, suggesting that alternative non‐competent hosts are an effective source of “wasted” bites. This phenomenon could be appreciated because we took into account the actual number of infectious vectors instead of prevalence. Prevalence is more frequently reported in the literature, but generally leads to misleading interpretation of “spurious” dilution results (Dobson & Auld, 2016; Roche et al., 2013). Nevertheless, the dilution potential of a non‐competent host population depends on its relative abundance, degree of tick infestation, interspecific competition and competence; and this requires it to be evaluated in more detail at local scale, as shown in other pathosystems (Mysterud et al., 2019). Most likely this mechanism was detectable due to the assumption of partial species additivity, and the simulation of tick aggregation, which realistically modeled the uneven distribution of ticks across the entire host population (Levi, Keesing, et al., 2016).

In the pathosystem involving the generalist tick, this mechanism acted in synergy with the second dilution mechanism observed, susceptible host regulation (Figure 5b), which in our model was mainly caused by predation (Levi et al., 2012; Ostfeld & Holt, 2004). Yet, dilution due to susceptible host regulation was observed in the specialist tick pathosystem as well. Predators are able to maintain rodent prey at low abundance across a wide range of resource availability (Korpimäki et al., 2004), and have been found to have additionally indirect negative effects on transmission (Hofmeester et al., 2017). In particular, mustelids do not markedly contribute to increases in the tick population (Lorusso et al., 2011; Meyer‐Kayser et al., 2012), suggesting that preserving top‐down forces as predators might be beneficial in decreasing disease risk (Johnson et al., 2019). In different communities, top predators were also observed to regulate deer populations (Korpimäki et al., 2004).

The role of ungulates in transmission dynamics might be more complex than in our model (Kilpatrick et al., 2017), but the effect of large grazing ungulates on ground‐dwelling small mammals’ behavior was beyond our scope, and it is still under investigation (Navarro‐Castilla et al., 2017). Nonetheless, we tested the effect on pathogen transmission of a large‐sized host population, such as sheep or deer, as this is preferred by I. ricinus adults and considered to be a source of Lyme disease amplification (Eisen et al., 2012). Our results were in agreement with Huang et al. (2019): i.e., the nymphal infection prevalence reduction outweighed tick population amplification, and this was observed in particular in highly biodiverse contexts (Sprong et al., 2020). Considering the generalist vector system and assuming low density, despite the contribution to a minor DIN increase, the ungulate population diluted transmission in hosts. Further analyses are needed to investigate the relationship between deer abundance and human disease risk at different relative densities (Huang et al., 2019).

A less significant reduction of susceptible hosts (i.e., rodents) resulted from shrew competition, but we chose a conservative estimate, so this remains an important area of uncertainty that requires further investigation (likewise competition among rodent species). In a directly transmitted disease, shrew density has been observed to be negatively correlated to infection probability in the focal rodent host, although the density of the two species were positively related. This suggests that the role of shrews as diluters might include being a source of wasted bites for ticks and a reservoir host competitor, and they may also alter the reservoir host behavior and habitat utilization (Khalil et al., 2016). Information on this aspect might be crucial to better evaluate the role of such competitor species in cases such as our small mammal specialist vector pathosystem. Our results showed that, in the pathosystem characterized by the small mammal specialist tick, the role of shrews as diluter or amplifier was complex and fluctuated as parameters related to the vector or the pathogen varied.

Our results support the findings of Ruyts et al. (2016, 2018), who also identified a sort of hump‐shaped DIN response to the increase in community diversity, together with an overall increase of nymphal population, although they did not test the relationship mechanistically, and assumed that community diversity was correlated to habitat complexity. From their results, it might be expected that the public health risk associated with Lyme disease transmitted by I. ricinus in Europe is not decreased by higher diversity, in disagreement with what was found in the United States. This is likely due to the differences between the ecology of tick species involved in those pathosystems (characterized mainly by I. scapularis), and the differences in host communities. Nevertheless, in the United States, the relative density of the reservoir host was observed to modulate the effects of species richness on DIN and B. burgdorferi host prevalence. So much so that high biodiversity did not always reduce those effects, while the best predictors of tick abundance were deer abundance and temperature at ground level (Werden et al., 2014). This suggests that despite pathosystem and community differences, patterns of community (dis)assembly and the way in which relative densities of competent and non‐competent hosts respond to these patterns might be a common trait to investigate in the context of the diversity–disease relationship.

This relationship is likely to be nonlinear and may be unimodal with a peak at some intermediate level of diversity (Strauss et al., 2016), therefore the European and North American communities might be at different sides of the peak. Very degraded habitats or very low‐diversity communities might have such few host species remaining that the likelihood of some pathogens persisting is extremely low (e.g., Gray et al., 1995; Richter & Matuschka, 2006). Therefore, to reduce disease risk by conserving or restoring biodiversity, it is crucial to determine on which side of the peak a given location falls before attempting an intervention (Strauss et al., 2016). For example, in Britain and Ireland, low‐biodiversity tick habitats are common where grass‐fed livestock are predominant and can alone maintain I. ricinus populations; here, low biodiversity can decrease rather than increase the abundance of infected ticks, contrasting with what is observed in the United States (LoGiudice et al., 2003), and emphasizing the complexity of the diversity–human‐disease risk of B. burgdorferi s.l. in Europe (Gray et al., 2021). Whereas, in habitats where the overall diversity might be increased by nonnative species, these might represent tick population ecological boosters, potentially increasing disease risk (e.g., Doi et al., 2021).

Considering the differences we have observed modeling the two pathosystems, the variety of aspects of tick ecology might also be crucial to determine disease risk in such pathosystems. For example, I. ricinus is more prolific than I. trianguliceps (Bown et al., 2006). It determines an overall larger infectious tick population despite the density‐dependent fecundity reduction function modeled and the greater effect of dilution mechanisms. With regard to tick specialization, this has been observed to follow a pattern of global generalism and local specialism, which needs to be better investigated to better understand the circulation of tick‐borne pathogens and exposure risks for humans and livestock (McCoy et al., 2013). Anthropogenic disturbance, and consequent habitat degradation, seems to favor generalist over specialist parasites (just like it has been observed in free living species; Dharmarajan et al., 2021), and this might be the case for vectors and particularly for ticks as well. However, if biotic homogenization of parasite communities usually leads to an increase in disease emergence risk, then pathosystems characterized by generalist ticks could be considered to be the keystone for the dilution effect hypothesis. Another aspect to be taken into account is vector feeding preferences in multi‐host communities and how it aggregates on hosts with different transmission competence levels. It has been demonstrated that this changes with host diversity, decreasing overall pathogen spread due to nonrandom sorting of viruliferous vectors between preferred and non‐preferred host species (Shoemaker et al., 2019).

Our model presented some limitations; first, it did not account for co‐infection, i.e., multiple pathogens influencing relative prevalence (e.g., Diuk‐Wasser et al., 2016; Hersh et al., 2014). Second, it was not spatially explicit, i.e., it did not include links between tick population dynamics and habitat (Maaz et al., 2018), the effect of patch connectivity, individual dispersion, and meta‐populations on pathogen transmission (Cohen et al., 2016), nor the influence of spatial scale and latitude (Magnusson et al., 2020), but these were beyond the scope of the study. Yet, the importance of considering spatial scale in evaluating zoonotic risk has been extensively demonstrated, especially in the case of Lyme disease in the context of forest fragmentation (e.g., Allan et al., 2003; LoGiudice et al., 2008). Third, despite collecting empirical data from wild rodents, crucial epidemiological parameters remained uncertain in our models. Uncertain parameters make unreliable quantitative predictions. Nonetheless, their inclusion in models does highlight important knowledge gaps and this indicates which experimental studies have the most potential to improve modeling predictions (Johnson, de Roode, et al., 2015). However, despite the uncertainty we were still able to identify two distinct dilution mechanisms, the species responsible for each mechanism to occur, and the differences in terms of community‐complexity–transmission relationship in the two pathosystems characterized by generalist and specialist vectors. Focusing on the qualitative patterns of our results, we demonstrated that dilution and amplification effects are not mutually exclusive in the same community, and that they depend on the epidemiological metric under consideration. Our modeling approach was innovative as it expanded from purely additive versus purely substitutive community structures (Johnson, Ostfeld, et al., 2015), creating more ecologically sound communities and a framework that embraced both theoretical and observational aspects (Mihaljevic et al., 2014).

The consistency of our main results with other European studies, i.e., that the density of infected ticks is hardly affected by dilution, supports the idea that in Europe the best strategy to lower tick‐borne disease risk might be the reduction of the tick‐human contact rate. This may be achieved through actions such as guided visitor flows in forests or frequently mowed forest trails instead of trying intervention to generically increase biodiversity (Ruyts et al., 2018). Managing disease risk through dilution might be very difficult for several reasons: when hosts compete for resources reintroduction of diluters might cause undesirable effects on the density of focal hosts; or diluters might fail to control disease in focal hosts because they are constrained by competition (e.g., Becker et al., 2014; Lacroix et al., 2014); or in extreme cases spillover events might occur (Strauss et al., 2015). Additionally, local species interactions, which may potentially interfere with disease transmission, may also exact other ecological consequences. Strauss et al. (2016) found that the observed disease risk–diversity correlation was spurious, and derived from differences in habitat structure, thus they suggest that another option to manage disease risk is to restore habitats. Indeed, in the context of rodent‐borne diseases, if the goal is to minimize disease emergence/abundance without species‐specific knowledge, preserving habitat may be a useful strategy, as land‐use change has been connected to a significant rise in rodent zoonotic reservoir species abundance (Mendoza et al., 2020). This approach might be the most effective option, not only relative to human disease risk, but also in case the pathogen under consideration was to threaten host populations. According to our results, preserving species functional diversity, and in particular top‐down predation, was crucial for the dilution effect to occur in host populations, therefore conserving habitat structure might be effective in avoiding local predator extinctions, and to increase functional diversity in general. Future research might focus on establishing the reservoir role of host/vector species in different habitat types, the behavior of ticks and hosts in different forest types and vegetation features, and the influence of habitat specific properties to identify potential generalities (Kilpatrick et al., 2017).

In conclusion, this investigation into the influence of multiple and diverse ecological components on pathogen transmission identified distinct dilution mechanisms and the species responsible for each mechanism to occur. In both pathosystems simulated, we demonstrated that more complex communities led to fewer infectious hosts providing evidence for dilution. Nonetheless, especially in the system with the more generalist vector, I. ricinus, the increase of host availability led to an amplification of human disease risk via the increased density of infectious nymphs (DIN). This was important in the system with the specialist vector only when considering the addition of shrews (potentially alternative hosts); yet, our conservative approach in terms of interspecific competition might have hindered the dilution potential of the competitor population in that system. Hence, our results supported the observations that dilution and amplification effects might not be mutually exclusive in the same community or pathosystem (Halsey, 2019), and that the relationship between diversity and disease risk is nonlinear (Ruyts et al., 2016; Strauss et al., 2016). The interpretation of this relationship depends on the epidemiological metric selected according to the focus of the study (Johnson & Thieltges, 2010; Lou et al., 2014), and patterns of community (dis)assembly in terms of both species identity (i.e., functional role in the community) and relative variation of reservoir species density (Johnson et al., 2019). Instead of trying to find a one‐size‐fits‐all approach to understanding disease–diversity relationships, our results emphasize the value of using eco‐epidemiological modeling supported by empirical data collection to investigate specific pathosystems. In particular, to improve the assessment of wildlife and human disease risk, there is an urgent need for more refined mechanistic approaches and local studies to avoid global ecological fallacy (Kilpatrick et al., 2017). This is crucial to plan ecological interventions to improve wildlife management and zoonotic risk control policies, and to achieve current public health and conservation goals (Sokolow et al., 2019). Currently, without specific knowledge on any given pathosystem the best strategy to achieve these goals seems to be habitat structure conservation.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

Appendix S1

Click here for additional data file.

Appendix S2

Click here for additional data file.

Appendix S3

Click here for additional data file.

Data S1

Click here for additional data file.