Skip to main content

Using Bayesian state-space models to understand the population dynamics of the dominant malaria vector, Anopheles funestus in rural Tanzania

Abstract

Background

It is often assumed that the population dynamics of the malaria vector Anopheles funestus, its role in malaria transmission and the way it responds to interventions are similar to the more elaborately characterized Anopheles gambiae. However, An. funestus has several unique ecological features that could generate distinct transmission dynamics and responsiveness to interventions. The objectives of this work were to develop a model which will: (1) reconstruct the population dynamics, survival, and fecundity of wild An. funestus populations in southern Tanzania, (2) quantify impacts of density dependence on the dynamics, and (3) assess seasonal fluctuations in An. funestus demography. Through quantifying the population dynamics of An. funestus, this model will enable analysis of how their stability and response to interventions may differ from that of An. gambiae sensu lato.

Methods

A Bayesian State Space Model (SSM) based on mosquito life history was fit to time series data on the abundance of female An. funestus sensu stricto collected over 2 years in southern Tanzania. Prior values of fitness and demography were incorporated from empirical data on larval development, adult survival and fecundity from laboratory-reared first generation progeny of wild caught An. funestus. The model was structured to allow larval and adult fitness traits to vary seasonally in response to environmental covariates (i.e. temperature and rainfall), and for density dependency in larvae. The effects of density dependence and seasonality were measured through counterfactual examination of model fit with or without these covariates.

Results

The model accurately reconstructed the seasonal population dynamics of An. funestus and generated biologically-plausible values of their survival larval, development and fecundity in the wild. This model suggests that An. funestus survival and fecundity annual pattern was highly variable across the year, but did not show consistent seasonal trends either rainfall or temperature. While the model fit was somewhat improved by inclusion of density dependence, this was a relatively minor effect and suggests that this process is not as important for An. funestus as it is for An. gambiae populations.

Conclusion

The model's ability to accurately reconstruct the dynamics and demography of An. funestus could potentially be useful in simulating the response of these populations to vector control techniques deployed separately or in combination. The observed and simulated dynamics also suggests that An. funestus could be playing a role in year-round malaria transmission, with any apparent seasonality attributed to other vector species.

Background

Anopheles funestus is one of the major malaria vectors in Africa and is widely distributed across the continent [1, 2]. With the exception of Anopheles gambiae sensu stricto (s.s.), the species appears to have higher vectorial capacity than many other members of the Anopheles gambiae complex [3,4,5,6,7,8]. Anopheles funestus makes a higher contribution to transmission than An. gambiae sensu lato (s.l.) in numerous parts in sub-Saharan Africa [6, 9,10,11,12]; particularly in settings where An. gambiae abundance has plummeted due to either effective indoor-based vector control interventions [13, 14] or environmental change. It is hypothesized that An. funestus persistence despite the recent scale-up of insecticide–treated nets may have been facilitated by their earlier development of strong physiological resistance [15].

Anopheles funestus is typically grouped with An. gambiae s.l. when modelling transmission and formulating policies for malaria vector control [16,17,18]. The lack of explicit consideration of An. funestus ecology and transmission potential may be partially due to this species having been relatively neglected compared to An. gambiae. Comparatively the ecology of An. funestus s.l. is less well understood, and it is much more difficult to maintain under insectary or semi-field conditions [19]. However, this species has several unique ecological features, such as its different larval habitat and dry season persistence [20], that could give rise to distinct population dynamics and differentiate its response to core and supplementary interventions. For example, An. funestus prefers larger aquatic habitats that are semi-permanent or permanent throughout the year, and contain clear water with some emergent vegetation [20]. This differs from An. gambiae s.s. which generally prefer small temporary habitats, such as puddles, ditches or animal hoof prints [20, 21], or Anopheles arabiensis, which can breed extensively in rice fields and other sunlit open pools [1]. The use of more permanent larval habitats means that An. funestus has greater persistence through the driest periods of the year compared to An. gambiae [22], whose habitats evaporate quickly in the absence of rainfall [21, 23]. This ecological feature means that the seasonal phenology of An. funestus and its response to aquatic microclimate differs from An. gambiae [21, 22, 24]; and could thus generate differential response to seasonally-targeted interventions, such as Indoor Residual Spraying (IRS) and larviciding.

Differential use of aquatic habitats may also impact the relative importance of key intrinsic drivers of mosquito population dynamics such as density dependence. Density dependence in malaria vectors occurs during larval development as a product of competition for space and nutritional resources [25, 26]. In space-limited habitats, high larval densities can influence larval development rates and survival, but also subsequent adult fitness traits such as body size, survival, fecundity and mating success [27,28,29,30]. While there is evidence that density dependence is an important driver of An. gambiae population dynamics [25], the relative importance of this process for An. funestus is unknown. Given that larval crowding and competition are less likely within the larger habitats preferred by An. funestus, density dependence is hypothesised to may be less pronounced for this vector species. Quantifying the strength of density dependence is important to inform the ease with which vector populations can be suppressed and how quickly they can recover [26, 27, 31].

Models of vector population dynamics and their response to interventions must be parameterised by reliable estimates of their demography and fitness. For vectors in the An. gambiae complex, such estimates are often acquired from insectary and semi-field studies [32,33,34,35,36,37] as well as field studies. Similar data has been difficult to obtain for An. funestus because of its poorly understood ecology and the difficulties of creating laboratory colonies; which so far has been achieved on only two occasions [11, 19, 38]. State-space models (SSM) provide an alternative approach to indirectly estimate these parameters by fitting a population dynamics model to observed time series data [39, 40]. These models are widely used in other fields of ecology and conservation biology to investigate the population dynamics of other animals [39, 41] and guide management decisions [41]. However, these models have so far had limited update in medical entomology. Given data on population fluctuations are available, these models can infer and estimate plausible demographic rates that could generate the observed dynamics [42].

SSMs are time-series models that distinguish between two stochastic components, namely, process (i.e. biological), which captures sequential dependencies between population components (e.g. eggs, larvae, pupae) and an observation component, which captures and corrects for biases and imprecisions in the data-collection process. Prior knowledge of the model parameters is used to bolster the information content of the time series data with existing expert or laboratory data and uncertainty in estimates. Population projections are then quantified on the basis of posterior probability distributions for parameters and population states. SSMs have recently been used to elucidate the dynamics and impacts of interventions on malaria vectors in laboratory and semi field populations [32, 33], but have not yet been applied to estimate An. funestus vector demographics in the wild. Here, an innovative SSM application was developed to describe the dynamics of wild An. funestus populations in Tanzania, and use it to assess extrinsic (environmental) and intrinsic (density dependence) drivers of their fitness and demography for the first time. Empirical data from laboratory experiments on An. funestus colonization [19] were incorporated together with the wild population data to develop an SSM. Time series field data collected in 2015 [6] and 2018 south-eastern Tanzania, and corresponding environmental information were used to validate the model. Specific aims were: (1) to accurately reconstruct the population dynamics, survival and fecundity of wild An. funestus populations in southern Tanzania, (2) quantify the effects of density dependence on the dynamics, and (3) to identify and quantify seasonal variations in An. funestus demography.

Methods

Time series data on wild An. funestus populations

Indoor densities of female An. funestus s.l. adults were recorded over 12 months of entomological surveys conducted in three villages (Tulizamoyo, Ikwambi and Sululu) in Kilombero (8.1539ºS, 36.6870ºE) and Ulanga (8.3124ºS, 36.6879ºE) districts, south-eastern Tanzania from June 2018 to May 2019 (Fig. 1). The villages were selected because of the high abundance of An. funestus s.l. within which An. funestus s.s. is the dominant sibling species (93%) [1]. Annual rainfall was 1200–1800 mm, and temperature, 20–32 °C. CDC Light traps [43] were used to sample host-seeking mosquitoes from 6 pm to 6am for 5 days per week, 4 weeks a month for 12 months in 10–15 houses per village. The houses were randomly selected and consent obtained from household heads, mosquito sampling were done repeatedly in these houses. The mosquitoes were sorted by taxa and sex, and females further classified as unfed, blood-fed or gravid. Daily climatic data (rainfall and temperature) were obtained from a weather station, approximately ~ 20 km from the farthest village.

Fig. 1
figure 1

A map depicting the locations of various study villages where mosquito sampling was carried out in 2015–2016 and 2018–2019 (Kindly prepared by Najat Kahamba)

To complement this, additional data on An. funestus were extracted from a 2015 dataset from three other villages in Ulanga district (Mavimba, Minepa and Kivukoni) (Fig. 1) [6]. These data were collected five days per week for a period of 12 months. This data allowed us to fit the model simultaneously to multiple time series so that it could learn hierarchically from An. funestus trajectories enfolding in different years and locations. This additional data has previously been described elsewhere and used to demonstrate the epidemiological dominance of An. funestus, which now contributes > 85% of all malaria infections in the region [6].

Prior information on life-history and gonotrophic cycle stages

Female An. funestus adults collected from the same three villages in 2018 were maintained in insectary conditions for one generation to estimate baseline fitness traits as already described in Ngowo et al. [19] and in Fig. 2. Data collected from this 1st-generation laboratory progeny included: (a) proportions of eggs that hatched into larvae, larvae that transitioned to pupae, and of pupae that emerged into adults, (b) the length of the transition periods (days) between life stages [(i) eggs to 1st instar larvae, (ii) 1st instar larvae to pupae, and (iii) pupae to unfed adult female (1 day post emergence)], (c) transition period of adult females between three different stages of their gonotrophic cycle, i.e. unfed, blood-fed and gravid.

Fig. 2
figure 2

Schematic representation of the state-space population model showing different life stages compartment (circles) and flows (arrows) of Anopheles funestus. Abundance data were only available for unfed, blood-fed and gravid stages. The model assumes that once a gravid mosquito has laid eggs, they return to the unfed stage. The annotations are described in Table 1. The model incorporates six life stages (eggs, larvae, pupae, unfed, bloodfed and gravid) of An. funestus

The gonotrophic cycle starts with ‘unfed’ females who transition to ‘blood-fed’ after obtaining a blood meal. In the wild, the first gonotrophic cycle usually starts after unfed females have mated [44]; which is assumed to happen soon after emergence. In insectary experiments, females had access to males immediately on emergence. As the blood meal is digested, blood-fed females transition into the ‘gravid’ state during which eggs develop. Gravid females then oviposit their eggs into aquatic habitats and return to the ‘unfed’ stage with the cycle begins again (Fig. 2). In the wild, the rate of transition between these gonotrophic stages is governed by both intrinsic and extrinsic environmental conditions including the availability of blood-meals and oviposition sites [45]. In the insectary, the first blood-meal (arm-feeding) was offered 5 days post emergence to ensure individuals had sufficient time for mating.

Per capita fecundity was defined as the number of eggs laid per fully bloodfed adult female. The proportions surviving between life-history stages or gonotrophic stages were calculated as the inverse of the number of days required to transit from one stage to the next:

$${s}_{m}={({j}_{m})}^{1/r}$$
(1)

Here, \(m\) is the life cycle stage, \(j\) is proportion survived as the percentage of the total from preceding life stage, \(r\) is the average number of days it took to transit from one life stage to another, and \(s\) is the daily survival within the stage (Table 1).

Table 1 Priors as used in the state-space population models of Anopheles funestus and the estimated posteriors mean and 95% credible intervals

Biological process components of the Bayesian SSMs

Daily survival

The daily survival of larval stages was assumed to be the same for all instar stages. Adult survival was assumed to be the same for unfed, blood-fed and gravid females. Survival probabilities (\({s}_{l}, {s}_{a}, {s}_{f}, {s}_{v}\)) were linked to their covariates through a logit transformation of linear predictors (here, subscripts \(l, a, f, v\) refer to larvae, unfed, blood-fed and gravid, respectively). Pupal (\({S}_{p}\)) and egg (\({S}_{e}\)) survival probabilities were considered to be independent of any climatic and density-dependent covariates, and were treated as a binomial distribution, with baseline rates assigned priors as described in Eq. 1617.

A range of covariates hypothesized to be associated with the demography of An. funestus were incorporated to allow baseline larval and adult survival to vary with environmental conditions. Rainfall (current and 1 week-lagged) and temperature were incorporated into the larval survival model. Rainfall regulates the availability and permanence of aquatic habitats, thus influencing both survival and carrying capacity of larval habitats [46]. Density dependence was incorporated into the model of larval survival [25] to assess whether this could improve the fit of the adult population dynamics model. Additionally, the speed of larval development was modelled as a function of temperature based on its known importance [47, 48]. The daily survival of larvae was thus defined as a function of daily rainfall (current and lagged), daily temperature and density dependence. The daily survival rates (lowercase \({s}_{l}\left(t\right)\)) of larvae were estimated through a logit transformation of linear predictors (uppercase \({S}_{l}\left(t\right)\)).

$${s}_{l}\left(t\right) =\frac{\mathrm{exp}\left({S}_{l}\left(t\right)\right)}{1+\mathrm{exp}\left({S}_{l}\left(t\right)\right)}$$
(2)

Specifically, \({S}_{l}\left(t\right)\) is written as a function of both intrinsic and extrinsic drivers:

$${S}_{l}\left(t\right) = {\beta }_{0}-{\beta }_{1}{R}_{\left(t-1\right)}-{\beta }_{2}{D}_{\left(t-1\right)}\left(1 -\frac{{\beta }_{3}{Q}_{\left(t-1\right)}}{\mathrm{max}\left(Q\right)}\right)+ {\beta }_{4}{T}_{\left(t-1\right)}- {\beta }_{5}{{T}_{\left(t-1\right)}^{2}}+ {\varepsilon }_{l,t}$$
(3)
$${\varepsilon }_{l} \sim Normal \left(0, {\sigma }_{l,t}\right)$$
(4)

Here, \({\beta }_{0}\) is the baseline daily larval survival on the linear scale. A survival probability prior was assigned under zero rainfall and average temperature (i.e. 27 °C) and then calculated the intercept of \({\beta }_{0}\) to reflect this prior information. When there is no effect of any environmental covariates (prior takes values between 0.80 and 1, Table 1). The coefficient \({\beta }_{1}\) quantifies the effect of current rainfall \((R)\); with the envisioned scenario being that higher \(R\) (i.e. flooding) tends to wash away larvae hence reducing the baseline survival [49]. This \({\beta }_{1}\) was defined by an informative gamma prior with shape = 5.382 and rate of 46.4 (Table 2) which permits anything from no rain effects to 100% mortality. The coefficient \({\beta }_{2}\) quantifies the effect of larval density at time \(t\) on larval survival. A monotonic negative relationship was assumed based on the biologically-plausible hypothesis that larval survival is reduced at high larval density because of resource competition and intraspecific cannibalism [29, 50]. This coefficient \({\beta }_{2}\) was defined by an uninformative gamma prior with shape of 0.5 and rate of 1 (Table 2), which allows the impact of density to range from no effect to complete annihilation.

Table 2 Priors for the intrinsic and extrinsic drivers of the population dynamic as used in the state-space model of Anopheles funestus and the estimated posteriors mean and 95% credible intervals

The term inside brackets in Eq. (3) represents the fact that density dependence needs to be modulated by the availability of larval habitat. The availability of suitable aquatic habitats for oviposition will increase with rainfall; thus potentially reducing the crowding of larvae into the remaining habitats that persist during the dry season. This hypothesis has been supported for An. gambiae s.l., where their seasonal population dynamics can be explained by models incorporating a rainfall-dependent carrying capacity [25]. Here, the coefficient \({\beta }_{3}\) was a proportion that captures the potential interaction between larval habitat availability (defined as the cumulative rainfall (Q) over the past week) and larval density (\(D\)). When rain in the recent week has been the maximum observed (i.e.\(Q=\mathrm{max}(Q))\), then (\(1 -\frac{{\beta }_{3}{Q}_{\left(t-1\right)}}{\mathrm{max}\left(Q\right)}\)) would be the smallest amount of density dependency experienced by An. funestus. The prior distribution for \({\beta }_{3}\) was defined by an upward-biased beta prior with mean 0.9 and variance of 0.01 allowing \({\beta }_{3}\) to have positive impact on larvae survival.

Additional covariates were incorporated to assess the role of temperature on larval survival (via the coefficients \({\beta }_{4}\),\({\beta }_{5}\)). The parameter \({\beta }_{4}\) captures the potentially positive effects of temperature on daily larval survival, which were defined by an uninformative gamma prior with mean of 1 and variance of 0.1 considering 27 °C as the optimal temperature (\(\rho )\) for maximum survival [51]. This prior allows temperature to vary from having no impact, to high positive impact on larval survival. Alternatively, the relationship between larval survival and temperature may be characterized by survival being reduced at low or very high temperature, and peaking in the middle [52]. The coefficient \({\beta }_{5}\) was incorporated to capture this potential curvilinear relationship but was dependent on \({\beta }_{4}\), to ensure that the optimum temperature was fixed at 27 °C, (\({\beta }_{5}=\frac{{\beta }_{4}}{2\rho })\). The prior of \({\beta }_{4}\) allowed extreme temperatures values away from the optimum to range from having no effect to generating 100% mortality. The parameter \({\varepsilon }_{l}\) is capturing the unexplained stochasticity associated with larval survival. This error term was defined by a normal prior with mean of 0 and a precision \({\sigma }_{l}\) from a gamma distribution with both shape and rate of 10.

The linear predictors for survival of unfed, blood-fed, gravid (uppercase \({S}_{a}\left(t\right), {S}_{f}\left(t\right), {S}_{v}(t)\)) and daily probabilities of survival (lowercase \({s}_{a}\left(t\right),{ s}_{f}\left(t\right), { s}_{v}\left(t\right)\)) were structured similarly to Eq. 3. The daily survival probabilities of adult stages were thus defined as the functions of daily temperature; such that an increase in temperature would result in an increase in the survival of all three adults stages and reduction in survival when temperature become lethal [22, 53]. The biological relationship between adult survival and temperature was assumed to be curvilinear [22, 53, 54].

$${S}_{a}\left(t\right)={\varphi }_{0}+{\varphi }_{1}{T}_{\left(t-1\right)} - {\varphi }_{2}{T}_{\left(t-1\right)}^{2}+ {\varepsilon }_{a,t}$$
(5)
$${S}_{f}\left(t\right)={\theta }_{0}+{\theta }_{1}{T}_{\left(t-1\right)} - {\theta }_{2}{T}_{(t-1)}^{2} + {\varepsilon }_{f,t}$$
(6)
$${S}_{v}\left(t\right)={\alpha }_{0}+{\alpha }_{1}{T}_{\left(t-1\right)} - {\alpha }_{2}{T}_{\left(t-1\right)}^{2}+ {\varepsilon }_{v,t}$$
(7)
$${{\varepsilon }^{*}}_{t} \sim Normal \left(0, {\sigma }_{t}^{*}\right)$$
(8)

Here, \({\varphi }_{0}, { \theta }_{0},\) and \({\alpha }_{0}\) refer to the baseline survival of unfed, blood-fed and gravid females respectively on a linear scale, under fixed temperature conditions of 27 ± 2 °C (insectary standard under which An. funestus have maximum survival [52, 54]), and assumes no blood meal limitation. The positive impact of temperature on all three life stages was represented by the coefficients \({\varphi }_{1},{ \theta }_{1}\) and \({\alpha }_{1}\) with an uninformative gamma prior with mean 12.5 and variance of 6.25. The coefficients \({ \varphi }_{2},{ \theta }_{2}\) and \({\alpha }_{2}\) correspond to the curvilinear effect of temperature on the survival of all three adult stages with their priors derived from the ratio between the linear coefficient and twice optimum temperature. This formulation ensured that the optimum temperature is fixed at (27 °C). The parameters \({\varepsilon }_{a}, { \varepsilon }_{f}, {\varepsilon }_{v}\) capture unexplained variation associated with survival during the distinct gonotrophic stages. These error terms (\({\varepsilon }^{*})\) were defined by normal priors with mean of 0 and a precision \({\sigma }^{*}\) from a gamma distribution with both shape and rate of 10 for unfed, blood-fed and gravid females.

Development between stages

The daily development probability from one life stage to the next was defined as the reciprocal of the development time (days) between the stages (assuming that all development times take longer than a day). An increase in temperature was assumed to reduce the development period of larvae [47, 53, 55, 56].

$$l\left(t\right)=\frac{\mathrm{exp}\left(L\left(t\right)\right)}{1+\mathrm{exp}\left(L\left(t\right)\right)}$$
(9)

Specifically, \(L\left(t\right)\) is written as the function of temperature covariates:

$$L\left(t\right)= {C}_{0}+{C}_{1}{T}_{\left(t-1\right)}$$
(10)

Here \({C}_{0}\) corresponds to the baseline daily development period on a linear scale defined by an informative beta prior with range defined in Eq. 1617 (Table 1). The coefficient \({C}_{1}\) explains the positive effect of temperature on larval development period, with its prior values derived from an uninformative gamma prior with mean 0.001 and standard deviation 0.001.

The development time for other life history stages (eggs and pupae) and the time between gonotrophic stages were assumed to be independent of temperature and other environmental covariates. The numbers of individuals (\({\mathcal{K}}_{m}\)) graduating from one stage to the next each day were modelled as a binomial process Eq. 11.

$${\mathcal{K}}_{m}\left(t\right)\sim Binomial({r}_{m},{ \mathcal{W}}_{m-1}(t)$$
(11)

Here the rate \({r}_{m}\) is a development probability as defined in Eq. 11 for \(m\) stage, with assigned informative prior values as described through a generic prior in Eq. 1617 (Table 1). Parameter \({\mathcal{W}}_{m-1}\left(t\right)\) refers to the number surviving the preceding life stage.

Fecundity

The number of eggs laid at each time step was drawn from a Poisson distribution whose rate was the product of per-capita fecundity (number of eggs laid by blood-fed An. funestus under insectary conditions (\({b}_{0}\)), a penalized rate for the egg survival (\({s}_{e})\), the number of gravid mosquitoes (\({V}_{(t-1)})\) and ratio of females-males (assumed to be 0.5) as assessed at the pupae stage [19].

$${b}_{t}= {\mathrm{exp}(b}_{0}+ {\varepsilon }_{b,t})$$
(12)
$$B(t)\sim Poisson\left(0.5{b}_{t}{s}_{e}{V}_{(t-1)} \right)$$
(13)

The error term \({\varepsilon }_{b}\) was defined by normal prior with mean of 0 and a precision from a gamma distribution with both shape and rate of 10.

Observation-derived components of the Bayesian SSMs

Observations of the abundance of adults (unfed, bloodfed, gravid) \(A\) at time \(t\) were modelled as a normal distribution with varying daily means \(\mathrm{\overline{a} }\) determined by the biological model and a precision \(\tau\) representing observation error. A fixed coefficient of variation (\(\xi )\) for the daily observation process was assumed and assigned an uninformative prior with values between 0.1 and 0.9 (Table 1). The CDC light trap typically samples mosquitoes from populations of unknown size, for which the daily catch rates are difficult to quantify independently. A parameter \(\vartheta\) was therefore incorporated both into the precision \(\tau\) and daily varying means to account for an observed weekly periodicity in adult abundance, which was otherwise hard to interpret. This parameter was allowed to vary both by day of the week \(j\) and between the two populations \(k\) (2015 and 2018–19 datasets). The \(\vartheta\) values were derived from a logit function \(\mathrm{ exp}\left(\rho \right)/(1+\mathrm{exp}\left(\rho \right))\), with \(\rho\) defined from the uninformative normal prior with mean and standard deviation of 0 and 10 respectively. Therefore precision \(\tau\) can be written as \(\frac{1}{{({\xi }_{t}{\overline{a} }_{t}{\vartheta }_{jk})}^{2}}\) for all the adult stages. Thereafter, the observation abundance was estimated as follows:

$$A(t)\sim Normal\left({\overline{a} }_{t}{\vartheta }_{jk}, \frac{1}{{({\xi }_{t}{\overline{a} }_{t}{\vartheta }_{jk})}^{2}}\right)$$
(14)

Trap bias

The trapping method (CDC Light traps) primarily targets unfed host-seeking mosquitoes [57]. Blood fed and gravid mosquitoes are assumed to no longer host-seek, and represent a small proportion 0.5–3% of total collections of females caught [6, 22]). To account for these biases in sampling, a new parameter of “trap-biasness”\(\omega\) was added in the observation model for both precision \(\tau\) and varying daily means \({\overline{a} }_{\eta }\). The prior values for \(\omega\) were estimated from independent studies from the same locations [6, 58], and ranged from 0.05–0.15 (Table 1), with variations between the two life stages \(\eta\). Therefore, the observation model for blood-fed and gravid (\({A}_{\eta })\) was rewritten by modifying Eq. 14 as follows

$${A}_{\eta }(t)\sim Normal\left({\overline{a} }_{\eta t}{\vartheta }_{jk}{\omega }_{\eta }, \frac{1}{{({\xi }_{t}{\overline{a} }_{\eta t}{\vartheta }_{jk}{\omega }_{\eta })}^{2}}\right)$$
(15)

Prior distributions

Since this model contains a large number of parameters, use of un-scaled informative priors restricted model convergence and mixing. A rescaled beta distributions of the informative priors [59] was opted and calibrated as follows:

$$Y=Beta \left(5, 5\right)$$
(16)
$$X={X}_{min}+Y\left({X}_{max}- {X}_{min}\right)$$
(17)

where \(Y\) is a dummy variable that takes values in the interval [0, 1] with mean of 0.5 and standard deviation of 0.15, selected to provide low likelihood at the values 0 and 1. The values of \({X}_{min}\) and \({X}_{max}\) define the range of the parameter of interest as dictated by the prior information. Since the information on priors was provided in form of mean (\(\mu\)) and standard deviation (\(\sigma\)), the values were defined as \({X}_{min},{X}_{max}= \mu \pm 2\sigma\). Survival, development period, trap-biasness, variability in daily catches and fecundity parameters were all assigned priors according to Eq. 17.

Model selection, model fitting and outputs

Model fitting was done using the R statistical software version 4.0.5 [60]. Population models were fitted using a Markov Chain and Monte Carlo sampling (MCMC) algorithm via the JAGS software [61] interfaced to R via the runjags package [62] (code provided in the Additional file 1: Fig. S1). To achieve convergence, the model with 6 chains was run in parallel for \({10}^{5}\) samples with a burn-in of \({10}^{5}\), keeping every 10th iteration for memory-saving reasons. Convergence was assessed by visual investigation of the trace plots, prior-posterior distribution using the coda package [63], effective sample sizes and the Gelman Rubin diagnostic [64]. Model comparisons were done using the deviance information criterion (DIC) [65], and the ones with the lowest DIC selected as the most preferred. The predicted and observed densities of An. funestus adult females were plotted to evaluate consistent prediction biases visually (Additional file 1: Fig. S2). Posterior means and 95% credible intervals for the key survival parameters, development period, density dependence, environmental covariates (temperature and rainfall) and fecundity were also reported to reveal different dynamical aspects of the system.

Results

Population trajectories and seasonal trends

Bayesian state-space model was used to describe the dynamics of wild populations of An. funestus. The full results, including summaries of posterior means for all the fitness and demographic parameters are reported in Table 1. The most parsimonious model (model-7, Table 3) included density dependence, and temperature and rainfall (current and lagged) impacts on larval survival, and the effect of temperature on larval development period. The only covariate that was not retained in the “model-7” was temperature impacts on adult survival. This model satisfactorily reconstructed the population dynamics of An. funestus in the study villages, with all environmental covariate relaxation applied based on DIC selection. Population trajectories were estimated for all six An. funestus life history and gonotrophic stages after accounting for potential impacts of environmental covariates and density dependence (Fig. 3).

Table 3 Model selection: Description of all models fitted with and without environmental covariates and their corresponding delta-Deviance Information Criterion ΔDIC
Fig. 3
figure 3

Reconstruction of the abundance trajectories for all the six life-stages. The red line indicates the mean posterior values and the respective 95% confidence intervals are shown in “sky-blue”. Left column (a, c, e, g, i, k) is data collected from June 2018 to May 2019 and right column (b, d, f, h, j, l) is data collected from Jan-Dec 2015. The grey area indicates the period with rainfall

These trajectories reflect the annual trend in abundance spanning in periods from low or no rainfall to high rainfall. All trajectories show a relatively high abundance of An. funestus right after the rainy season, followed by reduced but sustained abundance during the dry period for all the life stages. After accounting for observation biases during sampling, the observed abundance of unfed, gravid and blood-fed groups largely falls within the credible intervals of the predicted values (Fig. 4).

Fig. 4
figure 4

Observed vs. model estimated values for the three adult stages with data collected using CDC light trap both in May 2018 –June 2019 (left column- a, c, d) and Jan-Dec 2015 (right column- b, d, e). Red lines are the model estimated trajectories with “sky-blue” showing their 95% credible intervals. The blue circles are the observed values from the Light trap catches. Grey areas are the periods with rainfalls episodes

Survival and fecundity

Estimated An. funestus larval survival trajectories demonstrate substantial mean variability during the two seasons, with no clear pattern of seasonality (Fig. 5a and b), Table 1). Similarly, the survival trajectories of the adult stages (all gonotrophic states) were variable throughout the year, with daily survival rate ranging from 0.2 to 1.0 and not consistently differing between wet and dry seasons (Fig. 5c to h, Table 1). Per capita fecundity was estimated to be between 75 and 81 eggs per female An. funestus (Table 1). While the abundance of this species fluctuated seasonally, per capita fecundity remained consistent throughout the year (Fig. 5k and l).

Fig. 5
figure 5

Reconstruction of the survival trajectories for all the four stages (larvae, unfed, bloodfed, and gravid) which were affected by the environmental covariates. The two bottom rows show the larval development period and fecundity trends. Left column (a, c, e, g, i, k) is trajectories from June 2018 to May 2019 and right column (b, d, f, h, j, l) is from Jan-Dec 2015. Grey area is the period with rainfall. Y-axis shows the survival rates of different life stages and the bottom row (k, l) shows per-capita fecundity

Temperature was an important predictor of larval survival with a curvilinear relationship (ΔDIC = 138, Table 3, Fig. 6b), and that temperature has a positive monotonic relationship with larval development period (ΔDIC = 336, Table 3, Fig. 6a); with the larval development period estimated to last about 16 days on average. Additionally, daily rainfall was an important driver for the dynamics of An. funestus by reducing larval survival (ΔDIC = 5605, Table 3, Fig. 6c) with a negative monotonic relationship.

Fig. 6
figure 6

Relationship between environmental covariates and fitness parameters as estimated from the SSMs of population dynamic of Anopheles funestus

Effects of density dependence

Density dependence was the only intrinsic feature incorporated in this dynamic model of An. funestus. The model was able to converge efficiently without crashing when density dependency was removed, suggesting this process plays a detectable but relatively minor role population regulation when compared with extrinsic factors (ΔDIC = 222, Table 3, Fig. 6e). To verify this, a simulation was run and discovered that the estimated density dependence was actually quite low, even when simulating with a single dataset. The model fitting process also suggested the interaction parameter (\({\beta }_{3})\) between larvae density \((D)\) and one week cumulative rainfall \((Q)\) contributes to An. funestus dynamics by positively increasing larval survival (ΔDIC = 38, Table 3, Fig. 6d).

Discussion

Advanced methods used in mainstream ecological studies was adapted and fitted a state-space model (SSM) to field and laboratory data to accurately reconstruct population dynamics of wild population of An. funestus. The SSM inferred the trajectories of multiple life-cycle and gonotrophic stages of wild An. funestus females. This allows the reconstruction of the observed trajectories of larvae and adult females for the wild An. funestus in Tanzania for the first time. This analysis indicated that the dynamics of An. funestus were best explained in a model that included density dependency, temperature (curvilinear relationship) and rainfall (negative monotonic relationship) on larval survival, and temperature on the larval development period (positive monotonic relationship). In contrast, model fit was not improved by incorporating temperature dependency into adult survival (all gonotrophic stages). Anopheles funestus abundance vary seasonally between wet and dry but demographics rates (i.e. survivals, fecundity and development period) did not vary after accounting for the impact of environmental covariates and density dependence. These results are very useful for generating hypotheses about the nature and relative magnitude of drivers of An. funestus population dynamics in the wild. This model can be extended to include a component on malaria dynamics in humans; or to compare the efficacy and effectiveness of different interventions in combination or singly. This would allow more sophisticated evaluation of suitability of An. funestus-specific interventions; including prediction of the potential combined effect of strategies that acting at different life-cycle stages and/or target different demographic processes (e.g. survival versus fecundity).

Extrinsic covariates such as rainfall and temperature were all hypothesised to be the main drivers for the dynamics of this vector species. This study supports the hypothesis that rainfall is a significant driver of the population dynamics of wild An. funestus. Overall, the abundance of all life stages were relatively higher in rainy compared to dry periods of the year as previously documented [22, 66,67,68]. Rainfall covariates were directly included in the larval survival model since it is the only stage on which rainfall was hypothesized to have a significant impact. Daily larval survival as estimated by the SSM showed high variability both within seasons and across the year. There was support for a monotonic association between rainfall and larval survival; characterized as reduction in larval survival during periods of heavy rainfall [69]. Time lags have been used to assess rainfall impacts on the dynamics of these vectors [22]. Anopheles funestus abundance have been shown to be positively associated with the cumulative lag rainfall [22]. Here one week cumulative rainfall was included in the model to account for its effect on survival. Similar to other vectors of malaria transmission such as An. gambiae, rainfall have always been considered as the main factor regulating the dynamics, despite ecological differences between the two vector species [22, 46, 70, 71].

The SSM also provided support for the hypothesis that temperature is an important driver of An. funestus dynamics; although the nature of temperature effects was complex and variable between life history stages. For example, temperature was associated with both larval survival and development, but not adult survival or fitness. Furthermore the estimated impacts of temperature on larval ecology were complex; with the SSM suggesting a curvilinear relationship with survival but a positive monotonic impact with the larval development period. These findings validate the prior studies that demonstrated that temperature had a curvilinear influence on Anopheles larval survival, with a rise in temperature above/below the optimum lowering survival [47, 48, 72]. The larval period of An. gambiae is temperature dependent [55, 72]; thus this model incorporated a positive monotonic effect such that development is fastest when temperature is high and just below maximum threshold for larval development [72, 73]. In the final model the effect of temperature on gonotrophic stages was not found to be an important driver for the dynamic of this species thus left out during model selection process.

Little is known about the effect of density dependence on An. funestus due to its ecology and reliance on the large semi-permanent and permanent breeding habitats [20, 21]. However, density dependence is already well-known to be an important driver for dynamic of other malaria vectors like An. gambiae [25, 27,28,29, 31, 67, 74] and other non-malaria vectors like Aedes aegypti [75, 76]. Variations in densities during the aquatic stages of An. gambiae s.l. have been found to affect adult fitness [28, 77]. The SSM fit better when density dependence of larval survival was included, though the relative magnitude of this process was quite small and likely to have minor impact on the overall dynamics of An. funestus populations (Fig. 6e). These findings suggest that An. funestus populations are likely to be regulated more by extrinsic than intrinsic processes. These findings corroborate the original hypotheses about density dependence having a weaker regulatory role in this species on the basis of the types of larval habitats (i.e. larger and more permanent habitat [20, 78]), which can likely sustain higher resources and thus reduce competition than in An. gambiae. This is the first report documenting the role of density dependence on the dynamics of the wild populations of An. funestus. Now that colonies are becoming more feasible, more thorough investigation on the role of density dependence in the dynamic of An. funestus is prerequisite.

In addition to highlighting potential drivers of An. funestus populations, the SSM here generated plausible estimates of key demographic and life-history process in the wild. This model estimated that An. funestus larvae takes an average of 15.6–16.1 days to grow from first instar larvae to pupae; which is relatively long compared to the other major vectors in the An. gambiae complex (9–11 days [48, 55]). This apparently longer development period of An. funestus may be a product of their adaptation to more permanent, year-round breeding habitats that are unlikely to dry up; thus reducing selection for rapid development. The SSM estimated that the daily survival of wild An. funestus larvae could be as high as 0.95, compared to the 0.83 [0.80, 0.86] mean daily survival rate of the known vector of malaria transmission An. gambiae [48, 79]. This matches observations from insectary experiments in which An. funestus larvae have higher survival than An. gambiae [19, 38, 80]. Given the apparently higher rates of survival in An. funestus than in An. gambiae, these findings suggest that more lethal intervention may be required to control An. funestus both at larvae and adults stages.

The impact of any vector control largely depends on the ecology of the specific vector species. Differences in ecology between An. gambiae and An. funestus are likely to affect the relative impact of interventions. For instance, An. gambiae prefer breeding in small and temporary habitats which dry up quickly when there is no rainfall which is opposite to An. funestus habitats. Despite the fact that An. funestus habitats are "few, fixed, and findable" and might be easily targeted for larviciding [20] during the dry season, treating habitats such as rivers or bigger ponds could pose logistical challenges. The persistence of An. funestus throughout the year even during the driest periods suggest this vector is less seasonal compared to An. gambiae s.l., which experience much more dramatic “boom and bust” dynamics in relation to seasonal rains [6, 22, 68]. If interpreted together with the observation that survival estimates were not seasonal, the model suggests that this species is likely responsible for year-round malaria transmission throughout the year, while other species, which mostly occupy temporary habitats may be responsible for any apparent seasonality in transmission.

Models of vector population dynamics can provide a useful guide for the selection of optimal vector control strategies; particular through enabling more focal investigation of the benefits of seasonal or spatial targeting and use of combined versus single interventions. Despite its complexity, this population dynamics model provides a useful framework for investigation of the stability of An. funestus populations. With additional data, this model can be further refined to include additional modifications related to vector ecology and behaviour that may impact intervention (e.g. host choice and its impacts on fitness, predation during larval or adult phase and spatial components). Such further elucidation may increase the predictive accuracy of this SSM in specific contexts, but even the more general framework developed here have flexibility to introduce stage-specific mortality effects expected from different types of vector control [16, 32, 33]. For example, this framework could be used to model the impact of combined interventions including those that target adult females (insecticide-treated nets (ITNs), IRS) and larviciding; and assessment of how mortality varies with different coverage [16, 32]. It can also be used to investigate the possible response of vector population to climate change anticipated in Tanzania and other African countries. An important limitation of this study is lack of knowledge on what percentage of An. funestus mosquito population is sampled by the trap, which is important for understanding the relative magnitude of demographic stochasticity in modelled dynamics. This highlights the need to explicitly incorporate this source of uncertainty into vector and transmission dynamics; including the need for further calibration and standardization of the efficiency and biases associated with particular mosquito trapping methods.

Conclusions

This study used Bayesian State Space Models (SSM) parameterized with empirical data to quantify key demographic and fitness processes underpinning the population dynamics of An. funestus in Tanzania. This is the first use of SSM to understand the population dynamic of the wild vector of residual malaria transmission, An. funestus in Tanzania. The model structure allowed investigation of the relative importance of seasonally-varying environmental covariates (i.e. rainfall and temperature) and density dependence; providing some support for both processes although the magnitude of the former was much greater than the latter. The ability of this model to accurately reconstruct the seasonal dynamics and demography of An. funestus indicate its value for simulating the response of these populations to vectors control measures applied either individually or in combination. Additionally, the relatively limited evidence of seasonality in key fitness and demographic rates further corroborate evidence that this vector species can facilitate efficient year-round transmission of malaria. Finally, this model also highlights the clear importance of accounting for regional and daily observation biases when modelling mosquito population dynamics.

Availability of data and materials

The R-codes used for model development are freely available in the supplementary file. The time series data is available upon reasonable request to the author.

References

  1. Gillies M, Meillon D (1968) The Anophelinae of Africa south of the Sahara (Ethiopian zoogeographical region). Johannesburg: South Afr Inst Med Res. 1968. p. 343.

  2. Sinka ME, Bangs MJ, Manguin S, Rubio-Palis Y, Chareonviriyaphap T, Coetzee M, et al. A global map of dominant malaria vectors. Parasit Vectors. 2012;5:69.

    PubMed  PubMed Central  Article  Google Scholar 

  3. Russell TL, Govella NJ, Azizi S, Drakeley CJ, Kachur SP, Killeen GF, et al. Increased proportions of outdoor feeding among residual malaria vector populations following increased use of insecticide-treated nets in rural Tanzania. Malar J. 2011;10:80.

    PubMed  PubMed Central  Article  Google Scholar 

  4. Hunt RH, Brooke BD, Pillay C, Koekemoer LL, Coetzee M. Laboratory selection for and characteristics of pyrethroid resistance in the malaria vector Anopheles funestus. Med Vet Entomol. 2005;19:271–5.

    CAS  PubMed  Article  Google Scholar 

  5. Braack LE, Coetzee M, Hunt RH, Biggs H, Cornel A, Gericke A. Biting pattern and host-seeking behavior of Anopheles arabiensis (Diptera: Culicidae) in northeastern South Africa. J Med Entomol. 1994;31:333–9.

    CAS  PubMed  Article  Google Scholar 

  6. Kaindoa EW, Matowo NS, Ngowo HS, Mkandawile G, Mmbando A, Finda M, et al. Interventions that effectively target Anopheles funestus mosquitoes could significantly improve control of persistent malaria transmission in south-eastern Tanzania. PLoS ONE. 2017;12: e0177807.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  7. Lwetoijera D, Harris C, Kiware SS, Dongus S, Devine GJ, McCall PJ, et al. Increasing role of Anopheles funestus and Anopheles arabiensis in malaria transmission in the Kilombero Valley, Tanzania. Malar J. 2014;13:331.

    PubMed  PubMed Central  Article  Google Scholar 

  8. Garrett-Jones C, Shidrawi GR. Malaria vectorial capacity of a population of Anopheles gambiae: an exercise in epidemiological entomology. Bull World Health Organ. 1969;40:531–45.

    CAS  PubMed  PubMed Central  Google Scholar 

  9. Seyoum A, Sikaala CH, Chanda J, Chinula D, Ntamatungiro AJ, Hawela M, et al. Human exposure to anopheline mosquitoes occurs primarily indoors, even for users of insecticide-treated nets in Luangwa Valley. Southeast Zambia Parasit Vectors. 2012;5:101.

    CAS  Article  Google Scholar 

  10. Riveron JM, Chiumia M, Menze BD, Barnes KG, Irving H, Ibrahim SS, et al. Rise of multiple insecticide resistance in Anopheles funestus in Malawi: a major concern for malaria vector control. Malar J. 2015;14:344.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  11. Hunt RH, Edwardes M, Coetzee M. Pyrethroid resistance in southern African Anopheles funestus extends to Likoma Island in Lake Malawi. Parasit Vectors. 2010;3:122.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  12. McCann RS, Ochomo E, Bayoh MN, Vulule JM, Hamel MJ, Gimnig JE, et al. Reemergence of Anopheles funestus as a vector of Plasmodium falciparum in Western Kenya after long-term implementation of insecticide-treated bed nets. Am J Trop Med Hyg. 2014;90:597–604.

    PubMed  PubMed Central  Article  Google Scholar 

  13. Bayoh MN, Mathias DK, Odiere MR, Mutuku FM, Kamau L, Gimnig JE, et al. Anopheles gambiae: historical population decline associated with regional distribution of insecticide-treated bed nets in western Nyanza Province, Kenya. Malar J. 2010;9:62.

    PubMed  PubMed Central  Article  Google Scholar 

  14. Killeen GF, Seyoum A, Sikaala C, Zomboko AS, Gimnig JE, Govella NJ, et al. Eliminating malaria vectors. Parasit Vectors. 2013;6:172.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  15. Okumu F, Finda M. Key characteristics of residual malaria transmission in two districts in South-Eastern Tanzania-implications for improved control. J Infect Dis. 2021;223:S143–54.

    PubMed  PubMed Central  Article  Google Scholar 

  16. Fillinger U, Kannady K, William G, Vanek MJ, Dongus S, Nyika D, et al. A tool box for operational mosquito larval control: preliminary results and early lessons from the Urban Malaria Control Programme in Dar es Salaam. Tanzania Malar J. 2008;7:20.

    PubMed  Article  Google Scholar 

  17. Killeen GF, Chitnis N, Moore SJ, Okumu FO. Target product profile choices for intra-domiciliary malaria vector control pesticide products: repel or kill? Malar J. 2011;10:207.

    PubMed  PubMed Central  Article  Google Scholar 

  18. Griffin JT, Hollingsworth TD, Okell LC, Churcher TS, White M, Hinsley W, et al. Reducing Plasmodium falciparum malaria transmission in Africa: a model-based evaluation of intervention strategies. PLoS Med. 2010;7: e1000324.

    PubMed  PubMed Central  Article  Google Scholar 

  19. Ngowo HS, Hape EE, Matthiopoulos J, Okumu FO. Fitness characteristics of the malaria vector, Anopheles funestus, during an attempted laboratory colonization. Malar J. 2020;20:148.

    Article  CAS  Google Scholar 

  20. Nambunga IH, Ngowo HS, Mapua SA, Hape EE, Msugupakulya BJ, Msaky DS, et al. Aquatic habitats of the malaria vector Anopheles funestus in rural south-eastern Tanzania. Malar J. 2020;19:219.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  21. Gimnig JE, Ombok M, Kamau L, Hawley WA. Characteristics of larval Anopheline (Diptera: Culicidae) habitats in Western Kenya. J Med Entomol. 2001;38:282–8.

    CAS  PubMed  Article  Google Scholar 

  22. Ngowo HS, Kaindoa EW, Matthiopoulos J, Ferguson HM, Okumu FO. Variations in household microclimate affect outdoor-biting behaviour of malaria vectors. Wellcome Open Res. 2017;2:102.

    PubMed  PubMed Central  Article  Google Scholar 

  23. Minakawa N, Sonye G, Mogi M, Yan G. Habitat characteristics of Anopheles gambiae s.s. larvae in a Kenyan highland. Med Vet Entomol. 2004;18:301–5.

    CAS  PubMed  Article  Google Scholar 

  24. Kabbale FG, Akol AM, Kaddu JB, Onapa AW. Biting patterns and seasonality of Anopheles gambiae sensu lato and Anopheles funestus mosquitoes in Kamuli District, Uganda. Parasit Vectors. 2013;6:340.

    PubMed  PubMed Central  Article  Google Scholar 

  25. Russell TL, Lwetoijera DW, Knols BGJ, Takken W, Killeen GF, Ferguson HM. Linking individual phenotype to density-dependent population growth: the influence of body size on the population dynamics of malaria vectors. Proc R Soc B Biol Sci. 2011;278:3142–51.

    Article  Google Scholar 

  26. Nowicki P, Bonelli S, Barbero F, Balletto E. Relative importance of density-dependent regulation and environmental stochasticity for butterfly population dynamics. Oecologia. 2009;161:227–39.

    PubMed  Article  Google Scholar 

  27. Gimnig JE, Ombok M, Otieno S, Michael G, Vulule JM, Walker ED, et al. Density-dependent development of Anopheles gambiae (Diptera : Culicidae) larvae in artificial habitats. J Med Entomol. 2002;39:162–72.

    PubMed  Article  Google Scholar 

  28. Ng’habi KB, John B, Nkwengulila G, Knols BGJ, Killeen GF, Ferguson HM. Effect of larval crowding on mating competitiveness of Anopheles gambiae mosquitoes. Malar J. 2005;4:49.

    PubMed  PubMed Central  Article  Google Scholar 

  29. Koenraadt CJM, Majambere S, Hemerik L, Takken W. The effects of food and space on the occurrence of cannibalism and predation among larvae of Anopheles gambiae s.l. Entomol Exp Appl. 2004;112:125–34.

    Article  Google Scholar 

  30. Charlwood JD. May the force be with you: measuring mosquito fitness in the field. Ecol Asp Appl Genet Modif Mosq. 2003;2:47–62.

    Google Scholar 

  31. Henderson PA, Magurran AE. Direct evidence that density-dependent regulation underpins the temporal stability of abundant species in a diverse animal community. Proc R Soc B Biol Sci. 2014;281:20141336.

    Article  Google Scholar 

  32. Viana M, Ng’habi K, Lyimo I, Ferguson HM, Matthiopoulos J, Killeen G. Mesocosm experiments reveal the impact of mosquito control measures on malaria vector life history and population dynamics. Sci Rep. 2018;8:13949.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  33. Viana M, Hughes A, Matthiopoulos J, Ranson H, Ferguson HM. Delayed mortality effects cut the malaria transmission potential of insecticide-resistant mosquitoes. Proc Natl Acad Sci USA. 2016;113:8975–80.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  34. Lyimo IN, Haydon DT, Russell TL, Mbina KF, Daraja AA, Mbehela EM, et al. The impact of host species and vector control measures on the fitness of African malaria vectors. Proc Biol Sci. 2013;280:20122823.

    PubMed  PubMed Central  Google Scholar 

  35. Mweresa CK, Omusula P, Otieno B, van Loon JJA, Takken W, Mukabana WR. Molasses as a source of carbon dioxide for attracting the malaria mosquitoes Anopheles gambiae and Anopheles funestus. Malar J. 2014;13:160.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  36. Lyimo IN, Haydon DT, Mbina KF, Daraja AA, Mbehela EM, Reeve R, et al. The fitness of African malaria vectors in the presence and limitation of host behaviour. Malar J. 2012;11:425.

    PubMed  PubMed Central  Article  Google Scholar 

  37. Ferguson HM, Ng’habi KR, Walder T, Kadungula D, Moore SJ, Lyimo I, et al. Establishment of a large semi-field system for experimental study of African malaria vector ecology and control in Tanzania. Malar J. 2008;7:158.

    PubMed  PubMed Central  Article  Google Scholar 

  38. Okoye PN, Brooke BD, Hunt RH, Coetzee M. Relative developmental and reproductive fitness associated with pyrethroid resistance in the major southern African malaria vector, Anopheles funestus. Bull Entomol Res. 2007;97:599–605.

    CAS  PubMed  Article  Google Scholar 

  39. Newman KB, Buckland ST, Lindley ST, Thomas L, Fernández C. Hidden process models for animal population dynamics. Ecol Appl. 2006;16:74–86.

    CAS  PubMed  Article  Google Scholar 

  40. Millar R. Bayesian state-space modeling of age-structured data: fitting a model is just the beginning. Can J Fish Aquat Sci. 2000;50:43–50.

    Article  Google Scholar 

  41. Viana M, Cleaveland S, Matthiopoulos J, Halliday J, Packer C, Craft ME, et al. Dynamics of a morbillivirus at the domestic–wildlife interface: canine distemper virus in domestic dogs and lions. Proc Natl Acad Sci USA. 2015;112:1464–9.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  42. Rivot E, Prévost E, Parent E, Baglinière JL. A Bayesian state-space modelling framework for fitting a salmon stage-structured population dynamic model to multiple time series of field data. Ecol Modell. 2004;179:463–85.

    Article  Google Scholar 

  43. Mboera LEG. Sampling techniques for adult Afrotropical malaria vectors and their reliability in the estimation of entomological inoculation rate. Tanzan Health Res Bull. 2005;7:117–24.

    CAS  PubMed  Google Scholar 

  44. Gillies MT. The duration of the gonotrophic cycle in Anopheles gambiae and Anopheles funestus, with a note on the efficiency of hand catching. East Afr Med J. 1953;30:129–35.

    CAS  PubMed  Google Scholar 

  45. Rúa GL, Quiñones ML, Vélez ID, Zuluaga JS, Rojas W, Poveda G, et al. Laboratory estimation of the effects of increasing temperatures on the duration of gonotrophic cycle of Anopheles albimanus (Diptera: Culicidae). Mem Inst Oswaldo Cruz. 2005;100:515–20.

    PubMed  Article  Google Scholar 

  46. Fillinger U, Sonye G, Killeen GF, Knols BGJ, Becker N. The practical importance of permanent and semipermanent habitats for controlling aquatic stages of Anopheles gambiae sensu lato mosquitoes: Operational observations from a rural town in western Kenya. Trop Med Int Health. 2004;9:1274–89.

    PubMed  Article  Google Scholar 

  47. Bayoh MN, Lindsay SW. Temperature-related duration of aquatic stages of the Afrotropical malaria vector mosquito Anopheles gambiae in the laboratory. Med Vet Entomol. 2004;18:174–9.

    CAS  PubMed  Article  Google Scholar 

  48. Lyimo EO, Takken W, Koella JC. Effect of rearing temperature and larval density on larval survival, age at pupation and adult size of Anopheles gambiae. Entomol Exp Appl. 1992;63:265–71.

    Article  Google Scholar 

  49. Beier JC, Copeland R, Oyaro C, Masinya A, Odago WO, Oduor S, et al. Anopheles gambiae complex egg-stage survival in dry soil from larval development sites in western Kenya. J Am Mosq Control Assoc. 1990;6:105–9.

    CAS  PubMed  Google Scholar 

  50. Koenraadt CJM, Majambere S, Hemerik L, Takken W. Cannibalism and predation among larvae of Anopheles gambiae s.l. Entomol Exp Appl. 2004;112:125–34.

    Article  Google Scholar 

  51. Paaijmans KP, Heinig RL, Seliga RA, Blanford JI, Blanford S, Murdock CC, et al. Temperature variation makes ectotherms more sensitive to climate change. Glob Chang Biol. 2013;19:2373–80.

    PubMed  PubMed Central  Article  Google Scholar 

  52. Paaijmans KP, Imbahale SS, Thomas MB, Takken W. Relevant microclimate for determining the development rate of malaria mosquitoes and possible implications of climate change. Malar J. 2010;9:196.

    PubMed  PubMed Central  Article  Google Scholar 

  53. Blanford JI, Blanford S, Crane RG, Mann ME, Paaijmans KP, Schreiber KV, et al. Implications of temperature variation for malaria parasite development across Africa. Sci Rep. 2013;3:1300.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  54. Paaijmans KP, Thomas MB. The influence of mosquito resting behaviour and associated microclimate for malaria risk. Malar J. 2011;10:183.

    PubMed  PubMed Central  Article  Google Scholar 

  55. Kirby MJ, Lindsay SW. Effect of temperature and inter-specific competition on the development and survival of Anopheles gambiae sensu stricto and An. arabiensis larvae. Acta Trop. 2009;109:118–23.

    PubMed  Article  Google Scholar 

  56. Paaijmans KP, Blanford S, Bell AS, Blanford JI, Read AF, Thomas MB. Influence of climate on malaria transmission depends on daily temperature variation. Proc Natl Acad Sci USA. 2010;107:15135–9.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  57. Mboera L. Sampling techniques for adult Afrotropical malaria vectors and their reliability in the estimation of entomological inoculation rate. Tanzan J Health Res. 2006;7:117–24.

    Google Scholar 

  58. Msugupakulya BJ, Kaindoa EW, Ngowo HS, Kihonda JM, Kahamba NF, Msaky DS, et al. Preferred resting surfaces of dominant malaria vectors inside different house types in rural south-eastern Tanzania. Malar J. 2020;19:22.

    PubMed  PubMed Central  Article  Google Scholar 

  59. Gelman A, Jakulin A, Pittau MG, Su YS. A weakly informative default prior distribution for logistic and other regression models. Ann Appl Stat. 2008;2:1360–83.

    Google Scholar 

  60. R Development Core Team. R: A language and environment for statistical computing. R Found Stat Comput. 2021.

  61. Plummer M. JAGS : a program for analysis of Bayesian graphical models using Gibbs Sampling JAGS : Just Another Gibbs Sampler. 2003.

  62. Denwood MJ. runjags : an R package providing interface utilities, model templates, parallel computing methods and additional distributions for MCMC models in JAGS. J Stat Softw. 2016;71:1–25.

    Article  Google Scholar 

  63. Plummer M, Best N, Cowles K, Vines K. CODA: Convergence Diagnosis and Output Analysis for MCMC. R News. 2006.

  64. Gelman A, Rubin DB. Inference from iterative simulation using multiple sequences. Stat Sci. 1992. https://doi.org/10.1214/ss/1177011136.

    Article  Google Scholar 

  65. Spiegelhalter DJ, Best NG, Carlin BP, van der Linde A. Bayesian measures of model complexity and fit. J R Stat Soc Ser C Appl Stat. 2002;64:583–639.

    Article  Google Scholar 

  66. Kreppel KS, Viana M, Main BJ, Johnson PCD, Govella NJ, Lee Y, et al. Emergence of behavioural avoidance strategies of malaria vectors in areas of high LLIN coverage in Tanzania. Sci Rep. 2020;10:14527.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  67. Smith T, Charlwood JD, Takken W, Tanner M, Spiegelhalter DJ. Mapping the densities of malaria vectors within a single village. Acta Trop. 1995;59:1–18.

    CAS  PubMed  Article  Google Scholar 

  68. Charlwood JD, Vij R, Billingsley PF. Dry season refugia of malaria-transmitting mosquitoes in a dry savannah zone of east Africa. Am J Trop Med Hyg. 2000;62:726–32.

    CAS  PubMed  Article  Google Scholar 

  69. Kelly-Hope AL, Hemingway J, McKenzie FE. Environmental factors associated with the malaria vectors Anopheles gambiae and Anopheles funestus in Kenya. Malar J. 2009;8:268.

    PubMed  PubMed Central  Article  Google Scholar 

  70. Lyons CL, Coetzee M, Terblanche JS, Chown SL. Thermal limits of wild and laboratory strains of two African malaria vector species, Anopheles arabiensis and Anopheles funestus. Malar J. 2012;11:226.

    PubMed  PubMed Central  Article  Google Scholar 

  71. Zhou G, Munga S, Minakawa N, Githeko AK, Yan G. Spatial relationship between adult malaria vector abundance and environmental factors in western Kenya highlands. Am J Trop Med Hyg. 2007;77:29–35.

    PubMed  Article  Google Scholar 

  72. Paaijmans KP, Huijben S, Githeko AK, Takken W. Competitive interactions between larvae of the malaria mosquitoes Anopheles arabiensis and Anopheles gambiae under semi-field conditions in western Kenya. Acta Trop. 2009;109:124–30.

    PubMed  Article  Google Scholar 

  73. Lyons CL, Coetzee M, Chown SL. Stable and fluctuating temperature effects on the development rate and survival of two malaria vectors, Anopheles arabiensis and Anopheles funestus. Parasit Vectors. 2013;6:104.

    PubMed  PubMed Central  Article  Google Scholar 

  74. Muriu SM, Coulson T, Mbogo CM, Godfray HCJ. Larval density dependence in Anopheles gambiae s.s., the major African vector of malaria. J Anim Ecol. 2013;82:166–74.

    PubMed  Article  Google Scholar 

  75. Yang GJ, Brook BW, Whelan PI, Cleland S, Bradshaw CJA. Endogenous and exogenous factors controlling temporal abundance patterns of tropical mosquitoes. Ecol Appl. 2008;18:2028–40.

    PubMed  Article  Google Scholar 

  76. Hancock PA, White VL, Callahan AG, Godfray CHJ, Hoffmann AA, Ritchie SA. Density-dependent population dynamics in Aedes aegypti slow the spread of wMel Wolbachia. J Appl Ecol. 2016;53:785–93.

    Article  Google Scholar 

  77. Porretta D, Mastrantonio V, Crasta G, Bellini R, Comandatore F, Rossi P, et al. Intra-instar larval cannibalism in Anopheles gambiae (s.s.) and Anopheles stephensi (Diptera: Culicidae). Parasit Vectors. 2016;9:566.

    PubMed  PubMed Central  Article  Google Scholar 

  78. Tuno N, Githeko A, Yan G, Takagi M. Interspecific variation in diving activity among Anopheles gambiae Giles, An. arabiensis Patton, and An. funestus Giles (Diptera: Culicidae) larvae. J Vector Ecol. 2007;32:112–7.

    PubMed  Article  Google Scholar 

  79. Matthews J, Bethel A, Osei G. An overview of malarial Anopheles mosquito survival estimates in relation to methodology. Parasit Vectors. 2020;13:233.

    PubMed  PubMed Central  Article  Google Scholar 

  80. Takken W, Charlwood JD, Billingsley PF, Gort G. Dispersal and survival of Anopheles funestus and A. gambiae s.l. (Diptera: Culicidae) during the rainy season in southeast Tanzania. Bull Entomol Res. 1998;88:561.

    Article  Google Scholar 

Download references

Acknowledgements

We are grateful to member of IHI-VectorSphere insectary’s Joseph Mgando, Dickson, Mwasheshi, Godfrey Matanila, Neema Nombo, Rukiyah Mohammed, Johnson Kyeba and Khamis Bwanally for their contribution in different occasion during mosquito rearing, colonization and collection of wild population data. The earlier version of the model development benefited from the contributions of Luca Nelli, Fergus Chadwick, Heather McDevitt, Crinan Jarrett, Yacob Haddou and Prashanth Selvaraj. We also thank the three anonymous reviewers for their contributions.

Funding

This work was supported by Howard Hughes Medical Institute (HHMI)-Gates Foundation (Grants: OPP1099295) and Gates Foundation (Grants: INV-002138) both awarded to FOO and HMF which support HSN.

Author information

Authors and Affiliations

Authors

Contributions

HSN, FOO, HMF and JM designed the study. HSN, IHM and EEH conducted the lab experiment and collected data used in model development. HSN and JM conducted the formal model development, fitting and validation. HSN drafted the manuscript. FOO, HMF and JM reviewed the manuscript and provide supervision. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Halfan S. Ngowo.

Ethics declarations

Ethics approval and consent to participate

This study was approved by Ifakara Health Institute Review Board (Ref: IHI/IRB/No: 007-2018) and the Medical Research Coordinating Committee (MRCC) at the National Institute for Medical Research-NIMR (Ref: NIMR/HQ/R.8a/Vol.IX/2895). Individual verbal and written consent was also obtained from household owner where CDC light traps were places for collecting mosquitoes.

Consent for publication

Permission to publish was obtained from Medical Research Coordinating Committee (MRCC) at the National Institute for Medical Research (Ref. No: NIMR/HQ/P.12 VOL XXXIII/116).

Competing interests

The authors declare no conflicts of interest.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1.

SSMs model development R codes. Figure S2. Goodness-of-fit: Observed versus predicted unfed, bloodfed and gravid densities across all populations. Adjusted R-squared, intercept and slope values are from a linear model of the predicted against observed values. Dotted lines correspond to 1:1 line. Left column (a,c,d) is data collected from June 2018 to May 2019 and right column (b,d,e) is data from Jan-Dec 2015. Grey area is the period with rainfall. Figure S3.1. Prior (orange histogram) and posterior (blue histogram) distribution of the main baseline and observational parameters in the state-space model. Figure S3.2. Prior (orange histogram) and posterior (blue histogram) distribution of the main environmental covariates parameters in the state-space model.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Ngowo, H.S., Okumu, F.O., Hape, E.E. et al. Using Bayesian state-space models to understand the population dynamics of the dominant malaria vector, Anopheles funestus in rural Tanzania. Malar J 21, 161 (2022). https://doi.org/10.1186/s12936-022-04189-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12936-022-04189-4

Keywords

  • Anopheles funestus
  • State space model
  • Population dynamic
  • Seasonality
  • Abundance
  • Density dependence