The effect of interspecific competition on the temporal dynamics of Aedes albopictus and Culex pipiens

Background Aedes albopictus and Culex pipiens larvae reared in the same breeding site compete for resources, with an asymmetrical outcome that disadvantages only the latter species. The impact of these interactions on the overall ecology of these two species has not yet been assessed in the natural environment. In the present study, the temporal patterns of adult female mosquitoes from both species were analysed in north-eastern Italy, and substantial temporal shifts between abundance curves of Cx. pipiens and Ae. albopictus were observed in several sites. To understand which factors can drive the observed temporal shifts, we developed a mechanistic model that takes explicitly into account the effect of temperature on the development and survival of all mosquito stages. We also included into the model the effect of asymmetric interspecific competition, by adding a mortality term for Cx. pipiens larvae proportional to the larval abundance of Ae. albopictus within the same breeding site. Model calibration was performed through a Markov Chain Monte Carlo approach using weekly capture data collected in our study sites during 2014 and 2015. Results In almost half of observation sites, temporal shifts were due to competition, with an early decline of Cx. pipiens caused by the concurrent rise in abundance of its competitor, and this effect was enhanced by higher abundance of both species. We estimate that competition may reduce Cx. pipiens abundance in some sites by up to about 70%. However, in some cases temporal shifts can also be explained in the absence of competition between species resulting from a “temporal niche” effect, when the optimal fitness to environmental conditions for the two species are reached at different times of the year. Conclusions Our findings demonstrate the importance of considering ecological interactions and, in particular, competition between mosquito species in temperate climates, with important implications for risk assessment of mosquito transmitted pathogens, as well as the implementation of effective control measures. Electronic supplementary material The online version of this article (doi:10.1186/s13071-017-2041-8) contains supplementary material, which is available to authorized users.


Model calibration
The model accounts for the seven Cx. pipiens and four Ae. Albopictus, life stages of, namely eggs (Ec and Ea), the four Cx. pipiens larval instars (L1c, L2c, L3c, L4c) and one Ae. albopictus larval stage (La), pupae (Pc and Pa), non-diapausing Cx. female adults (Ac) and Ae. female adults (Aa). The ODE version of the model is based on the following system of equations:  life stages considered;    ,  , 1 , 2 , 3 , 4 ,  ,  ,  , , are the temperature dependent death rates associated with the different stages;  nEc and nEa are the number of eggs laid in one oviposition for a female of Cx. pipiens and Ae. albopictus respectively;  dAc and dAa are the length of the gonotrophic cycles of the two species;  Ka and Kc are density-dependent scaling factors driving the carrying capacity for the larval stages;  p is the probability (depending on daylight duration) that a fully developed Cx. pipiens pupa becomes a diapausing adult;  α represents the increase in mortality of Cx. pipiens larval stages due to competition with Ae.
albopictus. Its value is 0 if the daily temperature is below 15°C.;  βa and βc are the adult capture rates;  χC is a function of time defined equal to 1 when the trap is open and 0 otherwise;  Lc represents the total Cx. pipiens larval population, i.e. Lc=L1c+L2c+L3c+L4c;  Ca and Cc represents the cumulative number of captured female adult mosquitoes for Ae. albopictus and Cx. pipiens respectively.
Since only female adult mosquitoes are explicitly considered in the model, the term 1/2 in the equation for the adults accounts for the sex ratio (Delatte 2009, Vinogradova 2011. Moreover, given that diapausing Cx. pipiens females do not take blood meals before overwintering (Denlinger 2014), they are unlikely to be captured by using the considered traps. Consequently, only non-diapausing female adults are considered in the model.
We implemented model M as a discrete-time stochastic version, with time-step Δ = 1 day, in order to account for the stochastic nature of the processes. Precisely, the model is a Markov chain whose states represent the integer number of individuals in all developmental stages. Transition probabilities are built according to Poisson distributions whose means are obtained from the rate in system M. The seasonal dynamics of the mosquito population of each site is simulated from April 1 (corresponding to approximately one month before the first reported capture session) to October 31. Since, to the best of our knowledge, no data are available regarding the overwintering of Cx. pipiens and Ae. albopictus, we initialize the system with 100 non-diapausing Cx. pipiens adults and 10000 Ae. albopictus eggs. Preliminary model simulations showed no significant change of the model's behavior for different initial conditions. Mortality and developmental rates across different vector life stages have been modeled as a function of temperature as published in (Poletti 2011, Marini 2016. The probability p for a developed Cx. pipiens pupa to become a diapausing adult is a function of daylight duration as presented in (Marini 2016). The average number of eggs laid nEa and nEc per oviposition were fixed to 60 (Poletti 2011) and190 (Marini 2016) respectively. The duration of the gonotrophic cycle dAa is a function of temperature as in (Poletti 2011), while dAc is fixed to 5.54 days (Faraj 2006).
We assumed that, for each capture session, the number of captured female adult mosquitoes follows a Poisson distribution with mean obtained from the model; therefore, for each dataset, the likelihood of the observed data given a parameter set has been defined as where i runs over the number of capture sessions h, Ai (Ci) is the observed number of captured Ae. albopictus (Cx. pipiens) adults at capture session i and ̃( ) (̃( )) is the predicted number of captures of Ae. albopictus (Cx. pipiens) at capture session i simulated by the model with parameters θ = (α, β, a , ). Figure A shows an example of the temporal dynamics of larvae and adults for both species predicted by the two models in a given site and year. During spring months (April and May), the presence of Ae. albopictus larvae is limited by the relatively low temperatures and both models predict the same expansion of Cx. pipiens larvae and adults. Afterwards (beginning of June), increasing temperatures cause the rise of the Ae. albopictus population; consequently, the model with competition predicts a sharp fall of the larval (Fig. Aa) and adult (Fig. Ab) Cx. pipiens abundances. On the other hand, with the independent populations model, the decline of Cx. pipiens adults begins in late summer (August), when higher temperatures increase their mortality and progressively shortening photoperiods induce diapause in a growing number of newly emerged adult females. Ae. albopictus adults are better suited to higher temperatures and do not diapause, therefore their decline does not start until mid-September (Fig. Ab).

Model fit
The model fits the observed data quite well. In fact, about 75% of the recorded weekly captures lies in the 95% Credible Interval (CI) of model predictions with both assumptions.
In sites where the competition model was selected, the 95% CI of captured females predicted by the independent population model included the observed Cx. pipiens captures in 65% of data points overall, compared to 72% in the competition model. As expected, both models fitted equally well Ae. albopictus data, with about 80% of observations lying within the 95% CI of model predictions.

DIC and AIC analysis
We compared the goodness of fit of the model with interspecific competition against that of the model with independent populations (i.e. with α set equal to 0), using the Deviance Information Criterion (DIC): is the average value of D and var(D) is its variance.
Models with smaller DIC should be preferred. In fact, if the likelihood L is high (closer to 1) then ln is closer to 0. Moreover, var( ) increases with model complexity: in this way, the DIC penalizes models with a higher number of free parameters. We denote with DICα the value obtained with the interspecific competition model and with DIC0 the value associated with the independent populations model. Generally, model selection using the DIC criterion only requires a model to have a lower DIC than the alternative (corresponding to ΔDIC= DIC0 -DICα >0) (Spiegelhalter 2002). Considering the high stochastic noise in the capture data and the number of free parameters in our models, we conservatively restricted this criterion in such a way to minimize the risk of false positives on the existence of competition (Spiegelhalter 2002), by fixing a higher threshold on the minimum ΔDIC, i.e. ΔDIC>4. However, since the value of the threshold is arbitrary, we tested the robustness of our results by using a looser threshold, set to ΔDIC>2 as well as a different score function for model selection, namely, the Akaike Information Criterion, AIC (Burnham and Anderson 2002). AIC is defined as Where K is the number of the model parameters. Analogously to the DIC criterion, we selected models based on the value of on ΔAIC with respect to three standard threshold values (namely, 2, 4 and 7) Anderson 2002, Burnham 2011).
By loosening the DIC threshold, we included 5 additional datasets in the competition group (see Table A). The AIC yielded results very similar to the DIC, although slightly more conservative for corresponding values of the threshold. As reported in Table A

Estimates of the competition-dependent additional mortality
The ratio z = α/Kc defines the mortality rate of Cx. pipiens larvae due to each additional Ae. albopictus larva in the breeding site. Figure D shows that estimates of z for sites more strongly associated with competition (ΔDIC>4) are on average significantly higher (t-test p-value<0.05) than those associated to independent populations (ΔDIC<0): in other words, the competition model tended to be rejected when the estimated value of the competition-dependent mortality was closer to zero. The distribution of z values in sites with uncertain attribution (characterized by intermediate values of ΔDIC) was in between the two cases: this result mirrors the fact that competition has a nuanced, rather than an on/off effect on mosquito abundance.