Vector host-feeding preferences drive transmission of multi-host pathogens: West Nile virus as a model system

Jennifer E. Simpson, Paul J. Hurtado, Jan Medlock, Goudarz Molaei, Theodore G. Andreadis, Alison P. Galvani, Maria A. Diuk-Wasser

Abstract

Seasonal epizootics of vector-borne pathogens infecting multiple species are ecologically complex and difficult to forecast. Pathogen transmission potential within the host community is determined by the relative abilities of host species to maintain and transmit the pathogen and by ecological factors influencing contact rates between hosts and vectors. Increasing evidence of strong feeding preferences by a number of vectors suggests that the host community experienced by the pathogen may be very different from the local host community. We developed an empirically informed transmission model for West Nile virus (WNV) in four sites using one vector species (Culex pipiens) and preferred and non-preferred avian hosts. We measured strong feeding preferences for American robins (Turdus migratorius) by Cx. pipiens, quantified as the proportion of Cx. pipiens blood meals from robins in relation to their abundance (feeding index). The model accurately predicted WNV prevalence in Cx. pipiens at three of four sites. Sensitivity analysis revealed feeding preference was the most influential parameter on intensity and timing of peak WNV infection in Cx. pipiens and a threshold feeding index for transmission was identified. Our findings indicate host preference-induced contact heterogeneity is a key mediator of vector-borne pathogen epizootics in multi-species host communities, and should be incorporated into multi-host transmission models.

1. Introduction

Zoonotic pathogens cause significant mortality, morbidity and economic loss to human, livestock and wildlife hosts throughout the world, constituting an estimated 75 per cent of emerging infectious diseases [13]. The complex dynamics characterizing these systems, combined with costly ecological data collection, make predicting epizootics and spill-over to humans (or livestock) difficult [4]. It is therefore critical to identify the major determinants of pathogen transmission, both to develop tractable models capturing natural transmission dynamics, and to inform surveillance and control efforts in the field.

For vector-borne zoonoses, pathogen transmission potential of a host community is determined by the relative abilities of host species to maintain the pathogen and transmit it to a feeding vector (‘reservoir competence’) and by ecological factors influencing contact rates between hosts and vectors. Previous research has focused on how heterogeneity in reservoir competence among species influenced pathogen persistence, transmission [5] and human disease risk [6]. However, recent developments in molecular methods allow estimation of species-specific vector–host contact rates by identifying the source species from vector blood meals. The potential for vector host preference to influence transmission has been postulated by several empirical studies evaluating mosquito feeding preferences [79]. In this study, we addressed the question of whether heterogeneity in contact rates between vectors and hosts influenced transmission dynamics, independently of the relative reservoir competence of the hosts.

The central aims of this study were to quantify vector host-feeding preferences at multiple field sites and determine the extent to which these preferences shape West Nile virus (WNV) transmission dynamics. To accomplish these aims, we developed an empirically informed model using an epidemiological susceptible–infected–recovered (SIR) framework to describe daily enzootic transmission of WNV in a community of multiple avian hosts and a single mosquito vector, Culex pipiens (figure 1). We modelled transmission over 1 year at four field sites in Connecticut (CT, USA; electronic supplementary material, figure S1). We measured and incorporated into our model, seasonal mosquito abundances using data collected with Centres for Disease Control and Prevention (CDC) light traps, bird abundances using data from field surveys, and Cx. pipiens feeding behaviour through DNA analysis to determine the blood meal source [10]. Using sensitivity analysis, we evaluated the influence of all model parameters on two model outcomes, intensity and timing of peak transmission. We found that host-feeding preferences by Cx. pipiens were the most influential determinant for WNV transmission. To our knowledge, our model is the first arbovirus transmission model to be both informed and validated by field data at multiple sites, and thoroughly analysed using robust sensitivity analysis. Our findings demonstrate the critical role vector host preferences play in multi-host vector-borne pathogen systems. From an applied perspective, this model could be extended to determine transmission and emergence of other infectious vector-borne diseases involving multiple host species.

Figure 1.

Flow diagram for the WNV model. Transmission was simulated among two host types (vector-preferred and alternative). See text for model equations and table 1 for list of parameter values and sources.

2. Methods

(a) Field surveys

Four CT field study sites (site A in East Haven, site B in New Haven, site C in West Haven and site D in Stratford) were selected among 91 fixed locations where the CT Agricultural Experiment Station maintains mosquito surveillance traps. These sites have consistently experienced WNV activity in years since its initial introduction and have typical Culex spp. habitat, including waterways, parks, undeveloped wood lots and temporary wetlands in densely populated residential areas [1114]. A 500 m buffer centred on the traps delineated the sites.

Bird point counts were conducted during the breeding season at 5–11 randomly distributed positions in each site, located at least 200 m apart. Each position was visited once between 22 June and 8 July 2006, and all birds seen or heard within a 50 m radius of a survey position during a 10 min period were recorded according to established methods [15]. Survey results are provided in the electronic supplementary material. The total bird population at each site was estimated as the number of observed birds divided by the area of all survey points (birds per square metre), multiplied by the total area of the site (πr2). No adjustments were made for variation in detection probabilities (see the electronic supplementary material).

CDC-type light traps at the centre of each field site collected host-seeking mosquitoes overnight at least once in every 10 days from 5 June to 26 October. Traps were gathered the following morning and collections were sorted and identified to species with the aid of a stereomicroscope and descriptive keys [16]. Blood-engorged mosquitoes were collected using backpack aspirators. Mosquitoes were collected in plastic mesh-covered containers and transported alive to the laboratory where they were frozen, sorted and identified to species.

3. Laboratory methods

(a) Blood meal analysis

Blood-fed Cx. pipiens were processed individually to determine blood source and WNV infection status. The abdomens of blood-fed mosquitoes were separated from the body under a dissecting microscope, under sterile conditions to avoid cross-contamination. DNA was isolated from the abdominal contents of engorged mosquitoes using DNA-zol BD (Molecular Research Center, Cincinnati, OH, USA) according to established protocol [13,17,18]. Isolated DNA from the mosquito blood meal was screened using a polymerase chain reaction (PCR)-based technique taking advantage of the mitochondrial cytochrome b gene [13,17,18]. PCR products were sequenced at the W. M. Keck Foundation Biotechnology Resource Laboratory (Yale University, New Haven, CT, USA), a core DNA sequencing facility. The source of the blood meal was identified by comparison to the GenBank DNA sequence database (http://www.ncbi.nlm.nih.gov). The performance of the molecular-based assay was validated by isolating DNA from the blood of a number of known vertebrate species, subjecting them to PCR amplification, and sequencing [13].

(b) Feeding index

The ‘feeding index’ assesses the proportion of blood meals from a particular host species in relation to the proportional abundance of that species in the host community [10]. Therefore, a feeding index of 1 indicates opportunistic feeding habits, while a feeding index greater than 1 indicates preferential feeding. Confidence intervals for the feeding index were generated by bootstrapping with 1000 samples (STATA SE 10; Stata Corporation).

(c) Mosquito infection

Blood-fed mosquitoes were tested individually for infection with WNV. The head and thorax of individual mosquitoes were screened for viral infection by inoculating mosquito homogenates onto confluent Vero cell cultures as previously described [14]. RNA was extracted from virus isolates using viral RNA kit (Qiagen, Valencia, CA) and tested for WNV by real-time RT-PCR assays [19]. For estimating infection rate (IR), non blood-fed mosquitoes were pooled by species, collecting site, trap type, and date, and processed as above. The number of mosquitoes per pool ranged from 1 to 50. The IR and its skew-corrected 95% confidence intervals at each site were calculated using the Pooled IR add-in developed for Excel [20].

(d) Model

We constructed a model using a standard SIR epidemiological framework (e.g. [21]) to model enzootic transmission in a community of multiple avian hosts with a single mosquito vector for 1 year. The model incorporates WNV transmission between the vector, Cx. pipiens, and two categories of avian hosts: a primary avian species preferred by the vector, and all other bird species, referred to collectively as alternative avian hosts (figure 1). The resulting model equations (3.1)–(3.12) are given below.

Both primary and alternative avian hosts are further divided into classes of susceptible (S), infected (I) and recovered (R) individuals so the total population size is N = S + I + R. The variables for the preferred population of birds are indicated with subscript p, alternative hosts with subscript a, and vectors with subscript v. We assumed that mosquitoes do not recover from infection with WNV, so only susceptible (Sv) and infected (Iv) classes of vectors were modelled. The model is given by the system of differential equations:Embedded Image 3.1Embedded Image 3.2Embedded Image 3.3Embedded Image 3.4Embedded Image 3.5Embedded Image 3.6Embedded Image 3.7Embedded Image 3.8Embedded Image 3.9Embedded Image 3.10Embedded Image 3.11Embedded Image 3.12

We assumed that bird populations remain constant over the year, with balanced rates of recruitment and loss. To model seasonal dynamics of the mosquito populations, we forced the total number of mosquitoes to follow an empirically derived function for the population over time, r(t) in equation (3.7).

IRs for primary hosts, alternative hosts and vectors are written in terms of the force of infection (λp, λa and λv), respectively; equations (3.10)–(3.12). These equations describe WNV transmission between vectors and hosts, and incorporate the vector biting rate (v), reservoir and vector competence (β1 and β2, respectively), feeding preference (αv) and the number of infected individuals (I).

Model simulations for all four sites used parameter values and sampling ranges derived from published values [2233](table 1). Further information on model assumptions are described in the electronic supplementary material.

View this table:
Table 1.

Model parameters. (The rates have units per day.)

(e) Model validation

To validate the model, we compared the model output (daily abundance of infected Cx. pipiens) to the average monthly estimate of infected Cx. pipiens, or Vector Index (VI; described in [17]) which is calculated by the mean monthly Cx. pipiens abundance multiplied by the bias-corrected maximum-likelihood estimate (MLE) IR. A monthly estimate of the VI for each site for June through October 2006 was calculated (electronic supplementary material).

(f) Sensitivity analysis

We conducted sensitivity analyses to explore the model's behaviour and identify the key parameters driving transmission. We used a probabilistic method with Latin Hypercube Sampling (LHS), a type of stratified Monte Carlo sampling [3436]. Briefly, probability distribution functions were ascribed to each parameter (table 1), representing the degree of belief associated with each quantity. Sample values for each parameter are drawn from their distributions and LHS is used to better cover the whole space of possible combinations of parameter values with a moderate number of samples [34]. The model is run many times using the sampled parameter values to produce sets of model output. This sampling scheme generates input and output distributions useful in assessing model and parameter uncertainties in a ‘global' sense [34,35,37]. The partial rank correlation coefficient (PRCC) is a widely used and powerful test for parameter sensitivity [37] that calculates PRCC values for each input parameter (sampled by the LHS scheme) and model outcome variables. Since the PRCCs show which parameters have the largest influence on model outcomes, this analysis identifies which biological mechanisms are important determinants of intensity and timing of pathogen transmission.

4. Results

(a) Field results

To determine avian species abundances, a total of 489 birds representing 40 species were observed at point count stations in four field sites in CT (sites A–D). To parametrize the model, the abundance of American robins in each site was estimated for the total area of the site (site A = 88, site B = 133, site C = 25 and site D = 67 robins). The number of other bird species was calculated in the same manner (site A = 1613, site B = 1856, site C = 1100 and site D = 2483 birds).

We found that Cx. pipiens fed primarily on birds (139 of 147 blood meals, or 94.56%), and an average of 59.18 per cent (range 36.36–72.41%) of the avian blood meals were derived from American robins (T. migratorius) (figure 2). There were five species that were present in the blood meals but not observed during point counts: black-crowned night heron, Northern flicker, rose-breasted grosbeak, wild turkey and wood thrush. Since the blood meals and counts of non-robin species were pooled, we did not consider it necessary to make specific adjustments to the calculation of the feeding index, assuming that the proportion of blood meals on all other species is representative of the species abundances as estimated by the point counts. Using feeding indices, American robins were the preferred bird host at all sites (αv) >1 (site A = 6.70, CI: 2–18.90; site B = 8.53, CI: 5.49–16.05; site C = 31.89, CI: 14.52–47.14; site D = 24.83, CI: 11.92–98).

Figure 2.

Proportion of avian-derived Cx. pipiens blood meals at four study sites in CT, 2006 (n = 147). American robins comprised the greatest number of blood meals at all sites.

Mosquito abundance patterns as determined from light trap collections varied between sites. Abundance of Cx. pipiens peaked between late June and early August (electronic supplementary material, figure S2). The largest samples of Cx. pipiens were collected on 2 August (464) at site A, on 13 July (744) at site B, on 19 July (1276) at site C and on 26 June (419) at site D. The trap collections were used to parametrize the model. Further details on mosquito trapping results are provided in the electronic supplementary material.

(b) Model results

To address how strongly vector-feeding preferences influence WNV transmission relative to other biological mechanisms in the model, we varied the feeding index, αv, from baseline across a range of values (from 1 to 35). We observed its effect on two model outcomes; the intensity of transmission (measured as the peak number of WNV-infected Cx. pipiens at a site) and the timing of peak transmission (measured as the day of peak number of infected Cx. pipiens). Enzootic transmission of WNV was maintained at all sites in the model. Overall, sites with higher αv had more intense transmission occurring earlier in the season. Transmission peaked on day 295 (22 October) with 4.23 infected Cx. pipiens at site A, day 259 (16 September) with 52.72 infected Cx. pipiens at site B, day 193 (12 July) with 62.16 infected Cx. pipiens at site C, and day 200 (19 July) with 209.59 infected Cx. pipiens at site D.

We observed a consistent effect of αv on both transmission intensity and timing at all sites (figure 3). No transmission was maintained when αv was less than 6, and transmission intensity and timing was most sensitive to moderate αv values (7–15). For high αv (greater than 20), transmission intensity was approximately constant. This indicates that mosquito feeding preference has a large impact on transmission dynamics, and uncertainty in feeding preference estimates will most affect transmission intensity predictions when they are moderate.

Figure 3.

Model-predicted densities (log+1) of WNV-infected Cx. pipiens mosquitoes at four sites in CT in 2006. The influence of Cx. pipiens' preference for feeding on American robins (feeding index, αv) on peak abundance WNV-infected mosquitoes was compared. Site A (yellow diamonds) baseline αv was 6.70, site B (green squares) baseline αv was 8.53, site C (red triangles) baseline αv was 31.89 and site D (blue circles) baseline αv was 24.83. Transmission dynamics followed a consistent, nonlinear αv-dependent pattern across all sites.

These results are best illustrated when we simplify our model by assuming an equilibrium vector population of size (Nv). Under this assumption (see the electronic supplementary material) the basic reproductive number is given by:Embedded Image

R0 must be larger than 1 for disease to invade the disease-free steady state, and depends heavily on the value αv.

(c) Sensitivity analysis

We explored the contributions of each of the model parameters on the prediction precision of the two transmission outcomes (intensity and timing) using LHS and PRCC's. The sensitivity analysis showed that αv was the most important parameter for both transmission intensity (electronic supplementary material, table S1) and timing (electronic supplementary material, table S2) at all four sites. A strong positive correlation was observed between αv and number of WNV-infected Cx. pipiens at all sites, and a strong negative correlation between αv and day of peak transmission at three sites. The positive correlation between αv and timing of peak transmission at the remaining site, site A, suggests a lack of transmission when αv values are very low and a late date onset when αv values are just high enough to maintain transmission.

Mosquito biting rate (v) was ranked second in the PRCC analysis. The biting rate was positively correlated with transmission intensity at all four sites and negatively correlated with timing, except again at site A, where it was positively correlated for the same reason as for αv.

The PRCCs also show that the least influential parameter was the rate of new susceptible birds entering a site (b), either through births or immigration.

(d) Model validation

To determine whether model outcomes realistically represent enzootic transmission at the sites, we validated the model using a field-derived monthly mean Vector Index (VI) [38]. The VI is an estimate of the mean number of WNV-infected Cx. pipiens based on mosquito abundance and the MLE of the IR. A total of 16 676 adult female Cx. pipiens were collected and tested for WNV infection, of which 1267 individuals (53 pools) were from site A, 2823 (98 pools) from site B, 11 419 (273 pools) from site C, and 1167 (38 pools) from site D. The seasonal adjusted MLE (WNV-infected mosquitoes per 1000) for Cx. pipiens was 11.99 (13/53 pools tested positive for WNV), 4.93 (13/98), 10.91 (99/273) and 0.86 (1/38) at sites A–D, respectively. The first infected mosquitoes were identified on 29 June, the number of infected mosquitoes peaked in August or September and no positive mosquito pools were collected after 26 September at any of the sites.

We modelled the daily number of WNV-infected Cx. pipiens at each site when parameters were set to baseline. We show these results along with corresponding first and third quartiles from 1000 model runs (using LHS), and the monthly mean VI for validation in figure 4. The model predicts transmission outcomes falling within the confidence intervals of the monthly VI for three of the sites. Baseline model simulations matched the ranked order of the VI at three of the sites (low-level enzootic WNV transmission at site A, moderate at site B and higher at site C) but transmission was over-predicted at site D. The model was broadly able to corroborate outcomes, suggesting that our model realistically describes enzootic transmission at the study sites.

Figure 4.

Model results at baseline (black line), first (lower blue line) and third (upper blue line) quartiles (Q1 and Q3, respectively) from Latin hypercube sampling for the four sites in CT (A–D), along with 2006 Vector Indices (VI) represented as black triangles. (a) Site A, (b) site B, (c) site C, (d) site D. VI = monthly average maximum-likelihood estimate of infection prevalence × monthly average mosquito abundance. The lower (upper) error bar is derived from the lower (upper) limit of the MLE for infection prevalence multiplied by the monthly average mosquito abundance. Sites with low baseline αv had high variability, while sites with high baseline αv had less variability.

5. Discussion

Our model describes how vector–host contact heterogeneity owing to preferential mosquito feeding drives WNV enzootic transmission dynamics at four different locations in CT. Site-specific mosquito and bird abundances were represented in the model. From field data we identified host preferences for American robins by Cx. pipiens using the feeding index (αv) and with our model showed that this is the most influential parameter on both WNV transmission intensity and timing. This finding is consistent with host preferences by Cx. pipiens for American robins demonstrated by both field studies [7,8] and experimental trials [29]. As a consequence, it has been suggested that preference-induced contact heterogeneity [7,8,39] may influence pathogen transmission dynamics. Our model shows that strong host-feeding preference by Cx. pipiens for American robins is possibly a major driver for WNV transmission in CT.

Our model advances previous modelling efforts because it is both parametrized and validated using field data. We also used a parsimonious modelling approach: rather than modelling each site independently, we developed a single transmission model and subsequently incorporated site-specific field data from various locations of interest. This provides the advantage of a flexible model for exploring transmission dynamics both within and between multiple localities.

We identified a feeding preference threshold necessary for maintaining Cx. pipiens-mediated transmission of WNV (αv > 6). Increasing the degree of feeding on American robins resulted in more intense transmission occurring earlier in the season. The greatest influence was observed when preference was moderate (6 < αv <20); when αv > 20, transmission intensity and timing approached a plateau. This finding is consistent with how the threshold quantity Ro depends on αv (see the electronic supplementary material). Our model predicts maximum enzootic transmission potential at each site based on mosquito abundance and the relative abundance of preferred versus non-preferred hosts (figure 3). The level of host preference by Cx. pipiens, measured by αv, subsequently determines the timing and intensity of transmission and whether the plateau is reached. Inter-site differences in host preferences may be driven by landscape structure influencing the microhabitat availability of hosts to mosquitoes, availability of alternative hosts, etc; further research at a finer spatial scale may be warranted.

We used the VI to qualitatively validate that the model provides realistic estimates of infection in Cx. pipiens. The model agreed well with site-specific VI at sites A–C, predicting high densities of infected Cx. pipiens at site C, moderate densities at site B, and low densities at site A, but no agreement was observed between model and field data for site D. Lack of fit between model predictions and VI could be owing to simplified model assumptions, such as stable bird population densities, only two avian species represented, or similar infectivity and transmissibility for both avian hosts. On the other hand, the VI estimate is limited by field sampling biases and parameter estimation errors and low sample sizes, requiring using a monthly average of the VI. Unreliable IR estimates owing to small sample sizes for site D (38 pools tested, only 1 positive pool) may have led to the lack of fit with the model prediction.

Model results, along with IRs in Cx. pipiens in our study, further support the role of Cx. pipiens as the primary enzootic vector in this region. Culex pipiens has been established as an important enzootic vector by consistent isolations of WNV from mosquito trap collections [11,13,20], by its ornithophilic feeding behaviour [13,40,41], and associations between virus-infected mosquitoes and dead bird reports [4246]. This species has also been incriminated as a bridge vector in Illinois [47]. However, in our study, 94.65 per cent (139 out of 147) of blood meals examined were avian-derived and no human-derived samples were obtained, indicating that other mosquito species also serve as bridge vectors during epidemics in CT. Models to determine how feeding behaviours of bridge vectors influence epidemic transmission patterns would be useful to determine the risk of spill-over of WNV from birds to humans.

Our model builds on previously published epidemiological models that use generalized parameters to understand WNV transmission dynamics [27,48,49]. Previous models were rarely field-validated, since many were not empirically informed. Field-derived parameters are difficult to measure precisely and labour-intensive to obtain, but are vital to understanding real-world transmission dynamics scenarios. In contrast to previous models [27,4853], we used raw counts from mosquito traps to derive the rate of change in abundance. Other models have used differing assumptions to parametrize the vector population (reviewed in [54]). They may postulate some level of constant population growth [55], use a step-function [27], simulate growth over a season [51], or estimate production in a compartmental model consisting of egg, larvae and emerging susceptible adult stages [49,52]. By contrast, our method wholly captures the variation in the adult mosquito abundance data and does not require estimates for immature or other mosquito age classes.

WNV provides a suitable model for studying transmission dynamics of multi-host vector-borne pathogens because it is maintained in an enzootic cycle between mosquito vectors and multiple avian host species, with epizootic events in the United States every year. Given the strength of our findings, these results may be extended to other vector-borne diseases. A multitude of host-choice experiments provide convincing evidence that many malaria vectors, including Anopheles gambiae, prefer humans over cattle [5663], for example. Culex nigripalpus, a major vector for Saint Louis encephalitis in the southern United States, preferred chickens to bobwhite quail in a host-choice experiment in Florida, USA [64]. Host preferences are not exclusive to mosquitoes: Triatoma bugs that transmit Chagas' disease demonstrated strong feeding preferences for dogs over both chickens and cats in a recent study [65]. At least one species of tsetse fly (Glossina morsitans morsitans), vector of human (sleeping sickness) and animal trypanosomiasis in Africa, show feeding preferences towards buffalo and oxen but avoid other hosts such as waterbuck (reviewed in [66]). Extension of this model to other systems could determine how host preferences influence transmission dynamics.

In this study, we identified vector host preferences as the most important transmission parameter and quantified the contribution of preference-induced contact heterogeneity to enzootic transmission. The ‘dilution effect’ [5], whereby the presence of less pathogen-competent hosts in a more diverse community results in ‘wasted’ transmission events and therefore lower overall pathogen prevalence, relies on the assumption that vectors are generalists. Extensive application of the dilution effect to Lyme disease, for example, has been based on the assumption of no or low host preferences by ticks ([67], but see [68] for evidence of host preferences by ticks). Our study indicates that the host community composition experienced by the pathogen may be quite different from the host community we measured; incorporating host preferences and community reservoir competence into transmission models may be essential for modelling transmission dynamics and predicting epizootics.

Acknowledgements

We would like to thank the 2006 Yale field crew, especially Nick Bonomo for conducting bird surveys and to Elyse LeeVan, Anna Milkowski and Lindsay Rollend for their assistance in mosquito collection. We are also grateful to the CT Ag Station, John Shepard and Mike Thomas for providing equipment, mosquito trapping and identification, and Phil Armstrong for virus testing. We thank Angela Luis for helping with the R scripts to fit smoothed splines to the raw mosquito data. We thank the organizers of the EEID workshop (DEB—0722115, ‘Training Workshops on the Ecology and Evolution of Infectious Diseases’) for organizing the modelling training workshop that facilitated collaboration for several coauthors. We also thank Durland Fish, Jamie Childs, Corrine Folsom-O'Keefe and Paul Cislo at Yale University for their valuable scientific feedback and support. Funding was provided by the Ruth L. Kirschstein NIH Training Grant, CDC Laboratory Capacity for Infectious Diseases Cooperative Agreement no. U50/CCU6808-01-1, USDA-ARDS Cooperative Agreement 58-6615-1-218, USDA-ARDS Cooperative Agreement 58-0790-5-068, USDA Hatch Funds CONH00768 to the CAES, Yale Institute for Biospheric Studies Faculty Support Endowment Fund.

  • Received June 21, 2011.
  • Accepted July 28, 2011.

References

View Abstract