Co-feeding transmission facilitates strain coexistence in Borrelia burgdorferi, the Lyme disease agent.

Coexistence of multiple tick-borne pathogens or strains is common in natural hosts and can be facilitated by resource partitioning of the host species, within-host localization, or by different transmission pathways. Most vector-borne pathogens are transmitted horizontally via systemic host infection, but transmission may occur in the absence of systemic infection between two vectors feeding in close proximity, enabling pathogens to minimize competition and escape the host immune response. In a laboratory study, we demonstrated that co-feeding transmission can occur for a rapidly-cleared strain of Borrelia burgdorferi, the Lyme disease agent, between two stages of the tick vector Ixodes scapularis while feeding on their dominant host, Peromyscus leucopus. In contrast, infections rapidly became systemic for the persistently infecting strain. In a field study, we assessed opportunities for co-feeding transmission by measuring co-occurrence of two tick stages on ears of small mammals over two years at multiple sites. Finally, in a modeling study, we assessed the importance of co-feeding on R0, the basic reproductive number. The model indicated that co-feeding increases the fitness of rapidly-cleared strains in regions with synchronous immature tick feeding. Our results are consistent with increased diversity of B. burgdorferi in areas of higher synchrony in immature feeding - such as the midwestern United States. A higher relative proportion of rapidly-cleared strains, which are less human pathogenic, would also explain lower Lyme disease incidence in this region. Finally, if co-feeding transmission also occurs on refractory hosts, it may facilitate the emergence and persistence of new pathogens with a more limited host range.

1. Introduction

Wildlife and humans are frequently infected with multiple pathogens and their genetic variants, with important implications for human health (Walter et al., 2016; Hanincova et al., 2006; Andersson et al., 2013; Pedersen and Fenton, 2007). Tick-borne pathogens constitute a useful model to study mechanisms of pathogen coexistence because of high levels of co-infection in ticks and hosts (Walter et al., 2016; Andersson et al., 2013; Diuk-Wasser et al., 2016; Moutailler et al., 2016). In eastern North America, at least seven human pathogens are transmitted horizontally between Ixodes scapularis ticks and a wide range of vertebrate hosts (Diuk-Wasser et al., 2016; Donahue et al., 1987; Vuong et al., 2014). Furthermore, local populations of Borrelia burgdorferi sensu stricto, the most common I. scapularis-borne pathogen and the agent of Lyme disease in the United States, often contain genetically diverse strains (Wang et al., 1999). The ospC gene, which codes for the outer surface protein C (ospC), is a commonly used genetic marker to differentiate among strains (hereafter ‘ospC strains’) (Wang et al., 1999). The ospC protein is critical for establishing infection inside the vertebrate host (Schwan et al., 1995) and different ospC strains have been linked to differential B. burgdorferi dissemination and disease severity in humans. Elucidating the mechanisms driving strain frequency distribution is thus highly relevant to public health (Khatchikian et al., 2015; Wormser et al., 2008; Hanincova et al., 2013).

Multiple mechanisms simultaneously operate to maintain high levels of B. burgdorferi s.s. diversity. Balancing selection acting on B. burgdorferi outer surface proteins exposed to the host immune response is accepted as the main driving mechanism, although debate exists on whether this selection is adaptive (hosts differ in their relative fitness to different ospC types) (Hanincova et al., 2006; Brisson and Dykhuizen, 2004; Kurtenbach et al., 2006) or mediated by negative frequency dependent selection (rare amino acid types in ospC encounter weaker immune responses from the host) (Haven et al., 2011; Qiu et al., 1997). Strain coexistence can also be promoted by differential life history strategies, with varying fitness depending on the ecological context (Haven et al., 2012). By extrapolating from published data (Hanincova et al., 2008; Derdakova et al., 2004) to the early transmission period, Haven et al. (Haven et al., 2012) proposed that a “rapidly-cleared strain” of B. burgdorferi s.s. appeared to have higher early-phase transmission efficiency compared to an ‘in-host persistent’ strain. A model including competitive interactions indicated this early-phase transmission advantage would allow the rapidly-cleared strain to co-exist and even displace the in-host persistent strain under specific timing of immature tick feeding (phenology). An alternative trait-based approach to explore life-history strategies is loop analysis, where the fitness of alternative life history pathways can be directly evaluated by using the next-generation matrix to estimate the basic reproductive number, R0 (Davis and Bent, 2011; Tonetti et al., 2015).

In this study, we provide evidence that a previously unaccounted for transmission pathway, co-feeding transmission, occurs in the North American Lyme disease system (B. burgdorferi s.s./I. scapularis/Peromyscus leucopus). Co-feeding occurs when a pathogen is transmitted between two vectors (in this case I. scapularis nymphs and larvae) feeding closely in space and time via a localized rather than a systemic host infection (Higgs et al., 2005; Voordouw, 2015; Randolph et al., 1996). This transmission route minimizes pathogen exposure to the host immune system and reduces competition with other pathogens that infect the host systemically (Tonetti et al., 2015; Jacquet et al., 2015). Non-systemic transmission is the primary transmission mode for tick-borne encephalitis virus (TBEV) (Labuda et al., 1993) and was demonstrated in B. burgdorferi s.s. isolate ZS7 by I. ricinus ticks (Gern and Rais, 1996). However, the only previous study of co-feeding transmission of B. burgdorferi s.s. by its natural vector, I. scapularis, on its dominant host, the white-footed mouse (P. leucopus), found limited evidence of co-feeding transmission (Piesman and Happ, 2001). The validity of the latter study, however, is limited because larvae and nymphs were placed at the same time, therefore not allowing for the 36–48 h lag required for B. burgdorferi transmission (Randolph et al., 1996; Kahl et al., 1998; Crippa et al., 2002; Piesman, 1995). The role of co-feeding in the North American Lyme disease system is therefore ripe for re-assessment.

To examine the relevance of co-feeding in a realistic ecological setting (Randolph and Gern, 2003; Richter et al., 2002), we also measured the frequency of co-occurrence of nymphs and larvae in the same ear of field-sampled P. leucopus in a multi-site two-year study. We then extended a previously developed R0 model (Davis and Bent, 2011; Dunn et al., 2013, 2014) to include a co-feeding transmission pathway and parameterized it with our experimental and field data. While the experimental data is limited to two individual bacterial strains, our study illustrates a new mechanism than enhance coexistence of rapidly-cleared strains with in-host persistent ones. Our model also predicts the differential contribution of co-feeding to R0 across a range of tick feeding phenologies for two B. burgdorferi strains. Finally, the model can be extended to other tick-borne pathogen systems, in particular emerging pathogens, where co-feeding transmission in refractory hosts may significantly expand their host and geographic range (Jones et al., 1987; Ogden et al., 1998; Karpathy et al., 2016).

2. Methods

2.1. Laboratory methods: quantifying co-feeding transmission for two B. burgdorferi strains

2.1.1. B. burgdorferi isolates

Co-feeding transmission was assessed for two B. burgdorferi s.s. strains. Strain BBC13 (ospC type C) was isolated from a field-collected nymph from Block Island, Rhode Island, in 2013. This strain was typed as ospC type C (Tsao et al., 2013), which was the most abundant type during a two-year study in 2010 and 2011 (36.8% of fed larvae and 23.0% of nymphs) (States et al., 2014). The other experimental strain, B348 (ospC type E), was isolated from an erythema migrans lesion of a Lyme disease patient (Derdakova et al., 2004) and was shown to have low persistence in mice in previous experiments (Hanincova et al., 2008; Derdakova et al., 2004).

2.1.2. Transmission experiment design

We assessed pathogen transmission using the natural vector in the northeastern United States, I. scapularis, and the dominant reservoir host, the white-footed mouse. Six laboratory-raised (naïve) P. leucopus were purchased from an animal breeding facility (Peromyscus Genetic Stock Center, University of South Carolina) and housed at the Yale School of Public Health facilities. Three mice were infected with ospC type C strain and three with ospC type E strain (Fig. S1). On day 0 we established B. burgdorferi infection by placing eight infected laboratory-derived nymphs on the right ear of each mouse and inspected each mouse daily to confirm ticks had not moved to the left ear. On day 2 post-infestation, we placed 50 uninfected xenodiagnostic larvae on each mouse ear. On day 5, when larvae were nearly fully engorged but had not yet dropped off the mice, we removed all engorged larvae from each ear and kept them in an incubator at 25 °C and 95% humidity until they molted into nymphs. We inferred co-feeding transmission when at least one nymph molted from larvae placed on the same (right) ear as infected nymphs tested positive, while none of the larvae placed on the other (left) ear tested positive since the beginning of the experiment. We used xenodiagnosis as the standard and most sensitive method to detect disseminated infection (Marques et al., 2014; Shih et al., 1993; Kurtenbach et al., 1998). As soon as a single infected larva tested positive in the left ear, systemic transmission was assumed to have occurred. We repeated the steps of larval placement on days 7 and 14 and removal on days 10 and 17, respectively. These three periods represented the time during which non-systemic transmission via co-feeding was likely to occur. We additionally placed xenodiagnostic larvae on mice five more times on days 24, 31, 38, 52, and 101 to further assess transmission efficiency and duration of infection of the two strains. At these later points, approximately 100 clean larvae were placed on the ears of the mice and were allowed to feed to repletion, and engorged larvae were collected daily from water trays placed beneath the mouse cages. All mice were handled in accordance with Yale University IACUC procedures.

2.1.3. B. burgdorferi infection

We extracted DNA from nymphs and xenodiagnostic larvae using the Macherey-Nagel Nucleospin Tissue Kit (Macherey-Nagel, Duren, Germany). Individual ticks were homogenized in liquid nitrogen with sterile pestles. DNA extracts were screened in duplicate for B. burgdorferi DNA using quantitative PCR. Primer and probes targeted the 16S rDNA of B. burgdorferi following Hoen et al. (2009) using the Applied Biosystems 7500 real-time PCR machine (Applied Biosystems, Foster City, CA). Regression analyses using generalized estimating equations (GEE, using geepack package in R) were used to determine variations between sample groups, specifically between co-feeding transmission of strain B348 and systemic transmission of strain BBC13 during the first 17 days of the experiment and total transmission (co-feeding and systemic) of both strains for the duration of the experiment. An independent correlation matrix was included in the model to account for correlation among repeated samples. The model with the lowest QIC value was selected.

2.2. Field study: assessing potential for co-feeding transmission in nature

We estimated nymphal and larval burdens on white-footed mice and the frequency of co-occurring nymphal and larval burdens from two years of trapping small mammals at six forested sites across New England: three on Block Island, Rhode Island and three in Connecticut. From 13 May 2014 to 28 August 2014 (7 sessions) and 19 May 2015 to 27 August 2015 (7 sessions), we live-trapped small mammals twice monthly. Trapping sessions at each site consisted of three consecutive trapping nights. Connecticut trapping occurred at three forested sites: Hilltop (41.371872, −72.779478; 132 traps), Lakeside (41.371160, −72.773821; 144 traps), and Lord Cove (41.35906, −72.347; 120 traps). Block Island trapping occurred at three forested sites: Rodman’s Hollow (41.151258, −71.588489; 120 traps), Old Mill (41.163213, −71.589958; 60 traps), and West Beach (41.210015, −71.572009; 58 traps). Traps were placed in a grid formation 10 m apart.

Each Sherman box trap (30.5 cm × 7.5 cm × 7.5 cm) was baited with peanut butter, rolled oats, sunflower seeds, and two cotton balls; traps were opened at dusk and checked shortly after dawn each morning. We recorded host weight, sex, age, and carefully removed all ticks. We created three categories for nymph and larval burdens based on body location: left ear, right ear, or head and body region. Blood and tissue samples were also collected to examine infection status. After processing, mice were tagged and returned to the site of capture.

Overall larval and nymphal burdens were calculated as the daily total number of larvae and nymphs found on each newly captured mouse during each session. To assess the potential for co-feeding transmission, the daily co-aggregation by larvae and nymphs was calculated as the mean larval burden from each ear of newly-captured mice that also contained at least one nymph.

2.3. Modeling R0: Next-generation matrix including a co-feeding transmission pathway

2.3.1. R0 model specification

We expanded our previously developed model for the basic reproductive number R0 model for tick-borne pathogens (Dunn et al., 2013) to integrate co-feeding parameters derived from our laboratory and field experiments for the two B. burgdorferi strains (Table 1, Supplementary data, Supplementary text). We used a two-host next-generation matrix to represent transmission between the dominant host, P. leucopus, and larval and nymphal I. scapularis ticks. The transmission graph shows the two host types and the transmission pathways for B. burgdorferi (Fig. 1). Host type 1 represents a tick infected during its first blood meal as a larva, while host type 2 represents the vertebrate host. A larva may become infected by feeding on a systemically-infected host (k12) or via co-feeding transmission (k11), by feeding in close proximity to an infected nymph on the same mouse at the same time or soon after, before the mouse becomes systemically infected. Transmission from an infected nymph to a mouse occurs with probability k21. The next-generation matrix is:

Table 1

Definitions of parameters and point values in the R0 model

Function	Parameter	Description	Point Value				References
	dL	Days of attachment for larvae [days]	4				(Donahue et al., 1987; Brunner and Ostfeld, 2008; Diuk-Wasser et al., 2006)
	qN	Probability of nymph to mouse transmission	0.83				(Donahue et al., 1987; Gatewood et al., 2009)
	C	Proportion of ticks feeding on a competent host	0.5
	SN	Survival probability from infected larva to infected nymph	0.4				Field data; (Donahue et al., 1987; Randolph et al., 2000)
θ	ρ	Mean survival of mice [days]	133				(Harrison et al., 2011; Harrison and Bennett, 2012)
T (τ)		Systemic transmission efficiency τ days after an infected nymph has attached	Open symbols in Fig. 3				Lab data
T * (τ)		Co-feeding transmission efficiency τ days after an infected nymph has attached	Closed symbols in Fig. 3				Lab data
Phenology			CT2014	CT2015	BI2014	BI2015
Z̄N (t)	HN	Height of peak of nymphal activity	1.36	1.58	10.92	6.61	Field data
	τN	Timing of beginning of nymphal activity [days]	126	130	135	138	Field data
	μN	Time between beginning and peak of nymphal activity [days]	19.25	14.99	13.74	15.85	Field data
	σN	Shape parameter of nymphal activity	0.82	0.90	0.75	0.78	Field data
Z̄L(t)	HEt	Height of peak of total early larval activity	2.37	7.95	2.32	11.35	Field data
	τEt	Timing of peak of total early larval activity	150	148	158	160	Field data
	σEt	Shape parameter of total early larval activity	20.83	26.25	20.35	23.19	Field data
	HLt	Height of peak of total late larval activity	14.64	21.32	64.54	50.77	Field data
	τLt	Timing of beginning of total late larval activity [days]	189	186	190	191	Field data
	μLt	Time between beginning and peak of total late larval activity [days]	27.40	26.80	20.20	23.70	Field data
	σLt	Shape parameter of total late larval activity	0.48	0.28	0.50	0.54	Field data
Z¯L∗(t)	HEc	Height of peak of co-feeding early larval activity	1.93	7.44	1.59	6.90	Field data
	τEc	Timing of peak of co-feeding early larval activity	154.59	153.76	164.83	163.75	Field data
	σEc	Shape parameter of co-feeding early larval activity	16.97	19.81	17.99	23.77	Field data
	HLc	Height of peak of co-feeding late larval activity	9.58	10.89	64.33	19.13	Field data
	τLc	Timing of beginning of co-feeding late larval activity [days]	190.26	193.79	185.88	186.85	Field data
	μLc	Time between beginning and peak of co-feeding late larval activity [days]	18.24	21.72	28.19	29.51	Field data
	σLc	Shape parameter of co-feeding late larval activity	0.30	0.32	0.37	0.48	Field data

Fig. 1

Transmission graph showing three distinct transmission modes between two host types: white-footed mouse (Peromyscus leucopus) infected by an Ixodes scapularis nymph (k21), a larva infected by P. leucopus (k12), or a larva infected by a nymph through co-feeding transmission (k11). The arrows indicate these transmission modes between host types and correspond to the non-zero elements of the next-generation matrix.

The formulation of k11 was based on modeling of co-feeding by Davis and Bent (Davis and Bent, 2011), but here it was modified to incorporate the co-feeding transmission efficiency between infected nymphs and susceptible larvae as a function of days since attachment and the phenology of the two immature life stages.

The expected number of larvae an infected nymph will infect by co-feeding is given by: Where S is the proportion of fed larvae that survive and molt to become host-seeking nymphs and C is the proportion of those nymphs that feed on a competent host. The integral estimates the expected number of larvae that could become infected via co-feeding. The inner integral is a function of , the co-feeding larval phenology (see Seasonal tick burdens), T* (τ) is the probability of co-feeding transmission τ days after an infected nymph attachment, and d is the number of days of attachment for larvae. The upper bound of the inner integral, represented by l , indicates the length of the co-feeding period, which we determined experimentally. The outer integral weighs the value of the inner integral by the proportion of host-seeking nymphs a (t) that emerge and feed at time of year t (see Seasonal tick burdens).

Formulae for the remaining entries in the matrix follow Dunn et al. (Dunn et al., 2013): and Where q is the probability of transmission from an infectious nymph to a susceptible mouse, the inner integral is a function of Z̄(t), representing the larval phenology (see Seasonal tick burdens), T (τ) is the probability of systemic transmission from an infectious mouse to feeding larvae τ days post-infection and is the daily survival of mice calculated from their estimated lifespan ρ (Schug et al., 1991; Snyder, 1956). In k12, the inner integral provides the average number of larvae infected by systemic infection in the mouse after the initial period of co-feeding transmission, which starts after the lower bound l determined experimentally.

R0 is the dominant eigenvalue of K which is

2.3.2. Seasonal tick burdens

The calculation of R0 incorporates the phenology of overall tick burdens and the co-feeding phenology, which is the seasonal co-aggregation by larvae and nymphs. Following Brunner and Ostfeld (Brunner and Ostfeld, 2008) and Dunn et al. (Dunn et al., 2013), the functions Z̄ (t) and Z̄(t) were used to represent the mean overall nymphal and larval tick burdens and describes the mean co-feeding larval burden since the beginning of the year. These are formulated as:

Mean nymphal burdens over the year were fit to a lognormal curve while both total and co-feeding larval burdens were best fit by a right-shifted normal curve (‘early larval activity’) followed by a lognormal peak (‘late larval activity’). Tick activities can be categorized as nymphal (N), total early larval (Et), total late larval (Lt), co-feeding early larval (Ec), and co-feeding late larval (Lc). The parameter H denotes the height of the peak of tick activity, τ and μ determine the timing of tick activity and σ is the shape parameter for tick activity.

Parameterization was performed by fitting separate curves to tick burden data from trapping of P. leucopus from Connecticut and Block Island in 2014 and 2015. We assumed the number of larvae and nymphs X found on a host captured t days since the beginning of the year follow negative binomial random variables X = NB (m (t), α) with mean m (t) and dispersion parameter α. By maximizing the corresponding likelihood function with mean m(t) = Z̄(t), we obtained the values of parameters listed in Table 1.

2.3.3. Co-feeding contribution across tick feeding phenologies

To quantify the role of co-feeding across a range of different tick phenology scenarios and in comparison to a model not including co-feeding, we varied the proportion of larvae feeding during the early (p) and late (1 − p) larval peaks, keeping the total larval burden and the timing and the shape parameters constant. To derive realistic peak heights for both larval curves at any value of p, we fit a linear association between the observed number of co-feeding larvae and the product of nymphal and total larval feeding numbers as a measure of feeding overlap (Fig. S2). We also assumed nymphal phenology remained fixed across different larval phenologies as shown in (Diuk-Wasser et al., 2006; Gatewood et al., 2009); therefore, the linear association guaranteed the ratio of the number of co-feeding larvae to the number of total larvae was a constant for different phenologies. The heights of the co-feeding larval peaks where (H̃ (p) and H̃ (p)) for different tick phonologies (p) were then derived as: and where H̃ (p) and H̃ (p)are the calculated peak heights of overall larval burden for different tick phenologies.

2.3.4. Elasticity analysis

To assess the contribution of systemic and co-feeding transmission to R0, we performed an elasticity analysis (R software, package ‘popbio’). The elasticity of k11 and k12 represent the proportional change in R0 with proportional changes in co-feeding and systemic transmission, respectively, and were calculated as:

We assessed how R0, k11 and k12 varied across a range of tick phenologies (created by varying p) by replacing the original values of H , H , H , H with H̃(p), H̃ (p), H̃ (p) and H̃ (p), respectively.

3. Results

3.1. Quantifying co-feeding transmission for two B. burgdorferi strains

Co-feeding transmission occurred for both strains, but it was more commonly observed for mice infected with strain B348 (Supplementary data, Fig. S1). For strain BBC13 (mice #1-3), co-feeding transmission was observed only on day 5 from mouse 1, when 31.3% of larvae on the right ear tested positive for B. burgdorferi infection; no left ear ticks tested positive. Systemic transmission already occurred by day 5: from mice 2 and 3, 24.1% of left ear xenodiagnostic larvae were positive for B. burgdorferi infection. At day 10, systemic transmission occurred in all three mice: 60% of left ear ticks tested positive. On day 17, 100% of left ear larvae were positive. For strain B348 (mice #4–6), co-feeding transmission was observed for all mice on four independent instances. On day 5, mouse 4 remained uninfected while 11% of xenodiagnostic larvae on the right ears of mice 5 and 6 were infected; no left ear larvae were infected. On day 10, 57.1% of right ear xenodiagnostic larvae from mice 4 and 5 tested positive while all left ear larvae were negative for B. burgdorferi infection; however, mouse 6 exhibited evidence of systemic transmission with 2.9% of left ear larvae testing positive. By day 17, systemic transmission was observed in all mice, with 93.3% of left ear larvae testing positive for B. burgdorferi infection. After day 17, transmission efficiency to xenodiagnostic larvae declined, reaching 10.0% on day 101 (Fig. 2).

Fig. 2

Borrelia burgdorferi transmission efficiency to feeding ticks infected through co-feeding (strain B348) or systemic (both strains) transmission. The proportion of I. scapularis larvae infected via co-feeding with strain B348 during the first two time periods is represented by closed circles, the proportion of strain B348 infected through systemic transmission is represented by open circles and systemic transmission in strain BBC13 is represented by open diamonds (co-feeding was negligible for BBC13, data not shown). Co-feeding transmission is represented by larvae collected from the right ear (where nymphs were placed) that tested positive for B. burgdorferi infection while left ear larvae remained uninfected. From day 17 to the end of the experiment, only systemic transmission was assumed and the proportion of all positive ticks is represented. The dotted lines represent the assumption that transmission starts at zero and returns to zero seven days after the end of the experiment (day 108).

We observed greater transmission efficiency to xenodiagnostic larvae for BBC13 than B348 for the remainder of the experiment. By the end of the experiment (day 101), 40% of xenodiagnostic larvae were infected (Fig. 2). The GEE model showed that the total proportion of xenodiagnostic larvae infected with BBC13 (systemic transmission) and B348 (co-feeding transmission) were not significantly different before day 17. Total transmission (cofeeding and systemic) of BBC13 was significantly higher than transmission of B348 for the duration of the experiment (Fig. 2, Supplementary data, Table S1).

3.2. Field study: assessing potential for co-feeding transmission in nature

Tick burdens were higher on Block Island than Connecticut in both 2014 and 2015 (Fig. 3). The heights of the nymphal peaks were higher for Block Island in both years (H = 10.917 in 2014 and H = 6.605 in 2015) than those for Connecticut (H = 1.361 in 2014 and H = 1.575 in 2015; Fig. 3). Nymphs exhibited a single peak in early spring, on 24 May in Connecticut in both years and on 28 May 2014 and 2 June 2015 on Block Island. Overall larval burdens and co-feeding larval burdens peaked between 27 July and 3 August across all sites and years. We used the relationship between the number of co-feeding larvae and the product of nymphal and larval feeding numbers to derive the expected number of co-feeding larvae across a range of phenologies, an index of feeding overlap (p). When data from all sites and years were pooled, the Pearson correlation coefficient between the number of co-feeding larvae and the index of overlap was 0.64 and the coefficient of determination was 0.40 (Supplementary data, Fig. S2). Correlation coefficients were even higher when data were analyzed separately for different years and sites, ranging from 0.73 to 0.90 (data not shown).

Fig. 3

Ixodes scapularis tick feeding phenology at two sites for two years of our study: Connecticut 2014 (a), Connecticut 2015 (b), Block Island 2014 (c), and Block Island 2015 (d). Average larval and nymphal tick burdens are the average burdens on mice for each life stage. Co-feeding burdens represent the average larval burden per ear when nymphs were also present in the same ear.

3.3. Modeling R0: next-generation matrix including a co-feeding transmission pathway

Among the mice infected with each genotype, there was variation in the pattern of infection. For R0 modeling, we focused on those patterns most commonly observed (found in 2 out of 3 mice). Specifically, we assumed a mouse infected with strain BBC13 was already transmitting systemically by day 5, that a mouse infected with B348 experienced co-feeding transmission on days 5 and 10, and that all mice were systemically infected by day 17 (Supplementary data, Fig. S1). The R0 was greater for strain BBC13 than B348 and was greater for the Block Island than the Connecticut sites across all values of S (Fig. 4a and b). For B348, R0 estimates that included co-feeding transmission were always greater than systemic-only R0 values and this difference increased with increasing S (Fig. 4b).

Fig. 4

Estimates for the basic reproduction number R0 for B. burgdorferi strains BBC13 (a) and B348 (b) at Connecticut and Block Island sites in 2014 and 2015 across a range of estimates for tick survival multiplied by the probability of host attachment (S). We do not show co-feeding estimates for BBC13 because co-feeding was negligible for this strain. Systemic R0 estimates are represented by open symbols, co-feeding estimates by closed symbols.

For all sites and years, R0 decreased with p for strain BBC13 (Fig. 5a). The R0 (with or without including co-feeding) was relatively constant across p for strain B348 in Connecticut in both years but showed contrasting patterns on Block Island depending on whether co-feeding was included (R0 increased with p) or not (R0 did not change or slightly decreased with p) (Fig. 5b). At the limit (p = 1), the average R0 values across sites and years were almost identical for strains BBC13 and strain B348 (2.49 and 2.47, respectively). The co-feeding component of R0(k11) always increased with p (Supplementary data, Fig. S3a), while the systemic transmission component of R0(k12) did not show a consistent pattern with p (Supplementary data, Fig. S3b). To assess the relative contribution of both transmission components to R0, we performed an elasticity analysis for strain B348 (Fig. 6). The elasticity for k11 had a positive slope in relation to p, which indicates co-feeding increasingly contributed to R0 with increasing feeding synchrony. In contrast, k12 elasticity had a negative slope, indicating systemic infection contributed most to R0 when tick feeding was asynchronous.

Fig. 5

Estimates for the basic reproduction number R0 for B. burgdorferi strains BBC13 (a) and B348 (b) at Connecticut and Block Island sites in 2014 and 2015 across a range of tick feeding synchrony estimated by the proportion of early larvae (p). Systemic R0 estimates are represented by open symbols, co-feeding (B348 only) estimates by closed symbols. The open symbols with colored filling represent the actual proportion of early larvae we found at Connecticut and Block Island in 2014 and 2015.

Fig. 6

Variation across the full range of tick feeding synchrony (p) of the elasticity values for the co-feeding transmission component k11 (closed symbols) and the systemic transmission component k12 (open symbols) of the next-generation matrix for strain B348.

Finally, to assess the ecological space in which strain B348 could persist, we plotted the threshold curves (R0 > 1) for S = 0 – 0.1 and across the full range of tick feeding phenologies (p = 0 – 1). Threshold values for persistence were lower for Block Island, similar across years, and did not vary across p (Supplementary data, Fig. S4). In contrast, the greater threshold values for Connecticut varied across both years and had contrasting patterns across p.

4. Discussion

We demonstrated that co-feeding transmission can occur in the B. burgdorferi s.s/I. scapularis/P. leucopus system for a rapidly-cleared strain of B. burgdorferi. Furthermore, we assessed the ecological significance of the findings by integrating the controlled transmission experiments and two years of field-derived tick burden data into a model of R0. We also showed that co-feeding transmission increased pathogen fitness for a B. burgdorferi strain that is rapidly-cleared in the dominant mammalian host, P. leucopus. This effect is greater under simulated conditions of high tick feeding synchrony, as found in the midwestern United States and Europe (Gatewood et al., 2009; Randolph et al., 2000). Under very high synchrony, co-feeding increased the fitness of the rapidly-cleared strain (B348) to approximate R0 of the persistent strain (BBC13), increasing the likelihood that both pathogen strains coexist in the enzootic cycle.

Co-feeding has long been recognized as a key mechanism for the transmission of pathogens with a short duration of infection, such as tick-borne viruses (Randolph et al., 1996). The role of co-feeding in the maintenance of bacteria was considered more limited because a longer duration of infection enables more efficient systemic transmission. Previous R0 models of B. burgdorferi co-feeding transmission focused on persistent strains and assumed an infectious period of 120 days with constant host-to-tick transmission of 50% (Gern and Rais, 1996). Using these parameters, elasticity analyses found that B. burgdorferi strains invaded in the absence of co-feeding transmission and co-feeding only modestly increased R0 (Randolph et al., 1996; Harrison et al., 2011; Harrison and Bennett, 2012; Hartemink et al., 2008). We showed that co-feeding may have a stronger contribution to the transmission of rapidly-cleared strains such as B348. For these strains, co-feeding allowed the bacteria to be transmitted without systemic infection up to 10 days after the infectious nymph had detached from the host, which has been termed ‘extended co-feeding’ transmission (Randolph et al., 1996). Like viruses, these rapidly-cleared strains may rely more strongly on co-feeding transmission for their maintenance than previously thought (Voordouw, 2015).

Our laboratory experiment demonstrated the occurrence of co-feeding in the B. burgdorferi s.s./I. scapularis/P. leucopus transmission system, which was previously found to only occur at unrealistically high tick numbers (40 nymphs & >200 larvae) (Piesman and Happ, 2001). By placing larvae at 2, 7, and 10 after infestation with nymphs, rather than simultaneously, we were able to quantify the exact duration of the ‘extended co-feeding’ transmission period until systemic infection occurs (Randolph et al., 1996). Although the complexity of our study precluded a larger sample size, we were able to demonstrate co-feeding transmission for the rapidly-cleared strain using very strict criteria: that only larvae feeding on the ear where an infectious nymph was placed were infected and none from the other ear.

To explore how a realistic ecological context influences the role of co-feeding in transmission, we used extensive field data to parameterize the R0 model. During data collection, we specifically recorded the presence of nymphs and larvae on each ear of an individual mouse to allow for estimates of tick phenology and frequency of co-occurrence of two tick stages that would be comparable to our laboratory findings. Differences in tick burdens observed in the field also impacted the estimated role of co-feeding in R0. Greater burdens of nymphs and larvae on Block Island compared to Connecticut resulted in higher pathogen fitness across both years and all simulated phenologies. On the other hand, our model did not incorporate a potential effect of coinfection on transmission. Increased nymphal burdens may suppress susceptibility to infection if prior exposure to common strains suppressed infection with rapidly-cleared strains, as shown in (Derdakova et al., 2004; Jacquet et al., 2015).

Immature tick feeding synchrony had a complex effect on pathogen fitness. When ecological conditions were set above the threshold for persistence of both strains (above their respective value of S), R0 for B348 was either unchanged or slightly increased with greater feeding synchrony, while R0 for BBC13 decreased with increased synchrony. This discrepancy reflects the existence of two opposing effects of increased synchrony. Greater synchrony decreases the time between nymphal and larval feeding facilitating the transmission of strains with rapidly-cleared infection or relying on co-feeding transmission, such as strain B348. On the other hand, greater synchrony means an increased proportion of (naïve) larvae feed prior to (infectious) nymphs, resulting in mice becoming infected ‘too late,’ thus infecting a smaller proportion of the larval population and reducing R0. Increased fitness of rapidly-cleared strains and reduced fitness of persistent strains with increasing tick feeding synchrony would facilitated coexistence of these strains and may contribute to the higher strain diversity found in the upper Midwest of the US (Gatewood et al., 2009; Hamer et al., 2012).

Our laboratory experiments did not support the early-phase transmission advantage observed by Haven et al. (Haven et al., 2012) for rapidly-cleared strains by extrapolation of published studies. Instead, early phase transmission was similar for both strains, although based on different transmission pathways (systemic transmission for BB13 and co-feeding transmission for B348) (Fig. 2). We conclude that co-feeding provided a relative, rather than an absolute, benefit for the rapidly-cleared strain (B348). This is consistent with findings that some strains of B. afzelii - another member of the B. burgdorferi sensu lato complex, are highly efficient at all three transmission components of R0, which was reflected in their higher frequency in field collections (Tonetti et al., 2015). The greater synchrony of tick feeding in the midwestern US and Europe would further enhance the role of co-feeding by increasing the fit-ness of rapidly-cleared strains to similar levels of more persistent ones.

5. Conclusion

Our results indicate that the role of co-feeding transmission could be important in emerging areas, where it would allow B. burgdorferi s.s. to obtain some fitness benefits from incompetent vertebrate hosts (Voordouw, 2015; Randolph et al., 1996; Jacquet et al., 2015; Gern and Rais, 1996; Gern et al., 1998). Co-feeding may impact human health by contributing to the maintenance of a greater proportion of rapidly-cleared and less human virulent strains, if we assume strain competition limits the total B. burgdorferi infection intensity (Walter et al., 2016; Gatewood et al., 2009; Hamer et al., 2012). Furthermore, in the wake of the recent discovery of new I. scapularis-borne pathogens such as the Ehrlichia muris-like (EML) bacterium (Johnson et al., 2015) and B. mayonii in the midwestern United States (Pritt and Petersen, 2016), it will be increasingly important to quantify the impacts of co-feeding transmission, given the greater overlap of nymphs and larvae observed in that region (Gatewood et al., 2009; Karpathy et al. 2016). Finally, as climate change is predicted to increase feeding synchrony of the Northeastern U.S. tick populations to approximate that of the mid-western U.S. (Ogden et al., 2008), co-feeding may further increase the opportunity for transmission and maintenance of these emerging pathogens.