Our model of mosquito–Wolbachia dynamics allows two mosquito demographic traits, namely per-capita female fecundity and larval development times, to vary with changing larval density (Fig. 1; double red lines). We express these demographic traits as functions of larval density and allow their values to differ between infected and uninfected mosquitoes (see Methods). We estimate the form of these functions using observations from two independent populations of A. aegypti housed in field-cages (see below). The other demographic traits defined in our model (Fig. 1) are assumed to be density and time independent, and are estimated from direct observations either from our work or from previous studies. Details of the models and data are provided in the Methods.
Estimating density-dependent demographic traits
We refer to our two experimental field-cage A. aegypti populations as Population A (see [17]) and Population B (see Methods). wMel Wolbachia was introduced into both populations but in different ways. Population A was initiated with a cohort of uninfected adults and then at 2 months we began regular introductions of wMel-infected adults produced from larvae which had been reared in a separate facility (see Methods). Population B was initiated with a cohort of which 40 % of the adults were infected with wMel, and this population then received no further introductions of infected mosquitoes. The parameters of the functions describing density-dependent demographic traits were estimated using Bayesian Markov Chain Monte Carlo (MCMC) methods informed by (1) our observed abundances of the juvenile mosquito life stages, and (2) our observed wMel infection frequencies in first instar larvae and pupae over time in our two populations (see Methods).
Larval development times
The mean development times of the larval cohorts are much longer at higher larval densities, with similar fitted models describing this relationship for both populations (Fig. 2a). The variation in development times is also greater at higher larval densities, with both populations again showing similar relationships (Fig. 2b). The fitted models of density-dependent larval development times explain the major features of the dynamics of mosquito numbers and Wolbachia infection frequencies as measured by daily pupal surveys (Additional file 3: Figure S1.1, Additional file 4: Figure S1.2 and Additional file 2: Text S1).
For both populations, the predicted larval development times are highly variable amongst the individuals hatched in each cohort (Fig. 3). For Population A, the predicted mean development times of infected and uninfected larvae did not differ significantly (95 % credible interval (CI) includes 0; see Methods and Additional file 6: Figure S2.1). However, for Population B infected larvae developed faster than uninfected larvae for cohorts hatched in the first 5 weeks (Additional file 6: Figure S2.1; 95 % CI > 0 for cohorts hatched in weeks 4–8). For subsequent cohorts (hatched in weeks 9–15), the predicted mean development times of infected and uninfected larvae did not differ significantly. The faster development of infected larvae in the earlier cohorts may possibly be caused by genetic differences between the infected and uninfected mosquitoes (see the Discussion).
Per-capita female fecundity
Per-capita female fecundity declined strongly (by about factor of 10) with increasing larval density, with similar fitted models describing this relationship for both populations (Fig. 2c). Models incorporating these density-dependent relationships describe the major features of the dynamics of the average number of hatched larvae per week for both populations (Additional file 5: Figure S1.3), except for a short period in Population B (weeks 11–15), which is explored in Additional file 2: Text S1. For Population B, the predicted per-capita female fecundity over time did not differ significantly between infected and uninfected adults (95 % CI includes 0; Additional file 8: Figure S3.2, Additional file 7: Figure S3.1). For Population A, our data allow estimation of the per-capita fecundity over time of uninfected, but not infected, adult females (see Methods and [17]).
Modelling density-dependent population dynamics
We incorporated these estimates of density-dependent mosquito demographic traits into our model of mosquito–Wolbachia dynamics (Fig. 1 and Methods). We first use the model to predict their values when the mosquito population is at equilibrium and Wolbachia is not present. At equilibrium, the per-capita female fecundity λ
* and the larval development time distribution \( {\tilde{T}}_L^{*} \) depend on the level of larval density-dependent competition. We define the equilibrium net larval survival (from first instar to pupal eclosion) to be \( {\theta}_L\left({\tilde{T}}_L^{*},{\mu}_L\right) \), which depends on the mortality experienced throughout the larval stage, which we assume here occurs at a constant daily rate, μ
L
. Then, at equilibrium,
$$ \frac{\lambda^{*}{\theta}_L\left({\tilde{T}}_L^{*},{\mu}_L\right){\theta}_P{\theta}_{A_G}}{\mu_A}=1 $$
(1)
where θ
P
and \( {\theta}_{A_G} \) are the probabilities of surviving through the pupal stage and the early adult stage (during which females are too young to produce eggs), respectively (see Additional file 2: Text S1 and [26]). This expression simply states that each adult female produces, on average, one adult female offspring throughout its lifetime. Therefore, when the population experiences higher juvenile or adult mortality (higher μ
L
or μ
A
), the intensity of density-dependent competition at equilibrium decreases through changes in λ
*and \( {\tilde{T}}_L^{*} \).
Field mosquito populations are expected to experience higher density-independent mortality than our experimental field-cage populations, though mortality rates of juveniles and adults in the field are uncertain [6, 18]. We assume that the field population experiences density-independent mortality at a constant daily rate, μ
I
, during both the larval and adult stages. We further assume that this mortality acts in addition to the mortality occurring in our experimental populations, so that μ
A
= μ
A
c + μ
I
and μ
L
= μ
L
c + μ
I
, where μ
A
c and μ
L
c are the daily juvenile and adult mortality rates in our field-cage populations, respectively (Fig. 1 and Additional file 1: Table S1.1). We define the intensity of density-dependent competition experienced in the field population at equilibrium relative to the equilibrium derived from the field-cage experiments as C
I
* = L
I
*/L
0
* where L
I
* and L
0
* are the equilibrium larval densities given by a fixed value of μ
I
≥ 0 and μ
I
= 0, respectively.
As we increase the level of density-independent mortality, μ
I
, the intensity of larval competition declines and the equilibrium values of both the per-capita female fecundity (Fig. 4a; red line) and the mean larval development time (Fig. 4a; blue line) vary across the range of values observed in our experimental populations. Higher mortality, μ
I
, has a stronger effect on population size than the increased density-dependent fitness, and causes a steep decline in equilibrium population size (Fig. 4a; green line).
Predicting Wolbachia dynamics following field releases
We now explore how differences in the intensity of density-dependent competition in the field population affect the dynamics of Wolbachia following release. We model a typical strategy used in actual campaigns, where a fixed number of mosquitoes infected with wMel is released every week over 3 months [6, 9] and assume that the field population is at equilibrium when releases commence. The absolute number of mosquitoes that need to be released to achieve a given infection frequency is determined by the “release ratio”, or the size of each release divided by the initial wild population size. We set as a target that the frequency of infected adults must exceed 0.6 one week after the final release [3, 6] and calculate the minimum release ratio (MRR) required to achieve the target frequency. We then obtain the time taken for the Wolbachia to become established following the final release, T
E
, given that the release ratio is equal to the MRR, and defining establishment to occur when the infection frequency exceeds 0.95.
If larvae in the field population experience significant density-dependent competition, then it is likely that the released mosquitoes, reared with plentiful food, will have higher average female fecundity. If the average fecundity of the released females is equal to that which we observed at the lowest larval density (Fig. 2c), then they will have a very strong fecundity advantage (Fig. 4b; red line), especially when the intensity of density-dependent competition in the field population (C
I
*) is high. Therefore, the MRR increases by nearly an order of magnitude across decreasing intensities of competition C
I
* (Fig. 4b; solid blue line).
However, when competition is more intense, greater absolute numbers of released mosquitoes are required to meet the target frequency (Fig. 4b; dotted blue line) even though the MRR is lower. This is because the size of the field population is much larger due to the low density-independent mortality (Fig. 4a). Further, larval development periods are longer under more intense competition, which slows Wolbachia spread [17]. Therefore, the time to Wolbachia establishment following releases (T
E
) is longer (Fig. 4c). Thus, situations where competition is intense are disadvantageous for field release strategies overall, despite released females having a much stronger fitness advantage.
Fitness disadvantages in released mosquitoes
We now consider the situation where released mosquitoes are at a fitness disadvantage compared to wild-type individuals. Specifically, we assume that released individuals are more susceptible to a chemical insecticide used widely at the release site, leading to reduced adult survival. We assume that this fitness disadvantage is independent of Wolbachia infection so the disadvantage experienced by infected individuals declines over the generations following release due to introgression of wild-type genes. We assume that resistance is encoded by one allele at a single nuclear locus (R, resistant; S, susceptible) giving three genotypes: RR, SR, and SS. Homozygote resistant individuals are unaffected by the insecticide. Homozygote susceptible individuals experience a proportional reduction in daily adult survival of 1-c
SS
, with the cost, c
SS
= 0.8, assumed to be substantial. Heterozygote survival is reduced by a fraction of this amount to 1-f
SR
c
SS
(see Methods and [27]), where 0 ≤ f
SR
≤ 1. The field population prior to releases is assumed to be entirely composed of homozygote resistant individuals and all released mosquitoes are assumed to have the homozygote susceptible genotype. As before, we assume that the per-capita fecundity of the released mosquitoes is high, equal to that observed at the lowest larval densities in the field-cages.
Releasing mosquitoes with this substantial fitness disadvantage requires much higher MRRs to meet the target Wolbachia frequency (Fig. 5a). The MRR is approximately an order of magnitude higher when heterozygotes are fully resistant (f
SR = 0) and even greater when heterozygotes are also disadvantaged (f
SR > 0). When heterozygotes are disadvantaged, the MRR is reduced less under more intense competition (Fig. 5a). This is because the Wolbachia benefits greatly from the introgression of resistant alleles through the population following releases, which causes a rise in the average fitness of its hosts. Introgression is slower when larval development periods are lengthened by more intense competition, which inhibits Wolbachia spread, particularly when heterozygotes are disadvantaged.
The time taken for the Wolbachia to become established, T
E
, is more strongly affected by density-dependent competition when released mosquitoes have a fitness disadvantage, particularly if heterozygotes are also affected (Fig. 5b). When heterozygotes are fully susceptible to insecticides (f
SR = 1), T
E
is more than 3 years when competition in the field population is intense but less than a year when competition is low. This difference arises due to the strong effects on mosquito population dynamics of making large releases of Wolbachia-infected mosquitoes, particularly the reduction in the numbers of uninfected adults during the release period because of cytoplasmic incompatibility (Fig. 5c, d; black lines). The Wolbachia frequency drops rapidly to a low level after releases end due to the high fitness disadvantage (Fig. 5c, d; red lines), and cytoplasmic incompatibility therefore becomes less effective. If the population experiences low levels of density-independent mortality, μ
I
, it can recover rapidly from the suppression because the population growth rate is much higher at low densities. The population grows to almost pre-release levels before introgression of resistance allows the Wolbachia frequency to increase, and invasion is therefore very slow (Fig. 5c; black line). However, if density-independent mortality is high and population growth is less affected by density-dependence, the population is slow to recover from the reduction in density (Fig. 5d; black line). Introgression of resistance occurs while the population size remains low, allowing much faster Wolbachia invasion (Fig. 5d).