Mosquito husbandry

Anopheles stephensi mosquitoes (urban type form originally sourced from Walter Reed Army Institute of Research, Silver Spring, MD, USA) were reared at standard insectary conditions (27 °C ± 0.5 °C, 80% ± 5% relative humidity, and a 12 L:12D photoperiod) prior to the life table experiment, as described previously⁸. Briefly, we hatched immature mosquito larvae from eggs and placed 110 individuals into plastic trays (6 Qt., 12.4 cm × 34.6 cm × 21.0 cm) containing 500 mL of distilled water. Food (100 mg ground TetraMin fish flakes) was provided daily until most individuals reached the pupal stage. Pupae were rinsed and transferred to water-containing cups placed inside adult mosquito mesh cages for eclosion. For adult colony maintenance, An. stephensi were provided 5% dextrose and 0.05% para-amino benzoic acid (PABA) and fed whole human blood (O + , healthy male <30 years, Interstate Blood Bank, TN, USA) via water-jacketed hog intestine membrane feeders to support reproduction. All mosquito work was conducted in an arthropod containment biosafety level II facility under the purview of the University of Georgia’s institutional biosafety committee.

Experimental design

We adopted a similar experimental design as in Miazgowicz et al⁸., where we previously measured An. stephensi (urban type form) life history traits at six constant temperatures (16 °C, 20 °C, 24 °C, 28 °C, 32 °C, and 36 °C). Here, we programmed incubators (Percival; Perry, Iowa) to follow a Parton-Logan model⁷¹ for hourly diurnal temperature ranges (DTR) that are relevant for P. falciparum transmission in a natural setting (DTR of 9 °C or 12 °C) around five of the mean temperatures (16 °C, 20 °C, 24 °C, 28 °C, 32 °C ± 0.5 °C) explored previously⁸ (see Supplementary Methods). All other incubator settings (80% ± 5 RH, and 12 L:12D photoperiod) and experimental procedures were the same to allow for direct comparison between results. We were unable to include fluctuating treatments around 36 °C because the incubators did not reliably regulate their temperature above 40 °C. All experimental work for both studies was conducted during 2016-2018 at the University of Georgia (USA).

To generate a cohort of age-matched individuals, we collected pupae present at day nine post-hatch (when most immature mosquitoes reached the pupal stage) and placed them in an eclosion container within an adult cage for 24 hr. We provided a sugar solution (5% dextrose and 0.05% para-amino benzoic acid) to co-housed age-matched adults for three days prior to starting the lifetable experiment to permit mating. The lifetable experiment was initiated by providing females with an initial blood meal for 15 min, randomly sorting 300 blood-fed females into individual housing (16 oz. paper cup with mesh top), and then randomly assigning 30 individuals to each temperature treatment.

Each day until found dead, individuals were provided with a whole human blood meal for 15 minutes and inspected visually for imbibed blood to calculate the bite rate. This directly-measured bite rate was to determine if mosquitoes were biting more frequently than what would be inferred by the assumption of one bite per gonotrophic cycle. Oviposition sites (secured petri dish containing water saturated cotton and filter paper) within each individual housing were rehydrated and checked daily for eggs; if present, eggs were removed and counted. We terminated the experimental block when either all mosquitoes had died or when at most four mosquitoes remained alive at 16 °C.

The life table experiment for both fluctuation regimes was performed two independent times resulting in data from a total of 600 individuals (i.e., 60 individuals per combination of mean temperature and fluctuation regime). Life table data collected across constant temperatures from the previous study by our group consisted of 390 individuals across six constant temperatures⁸. In that experiment, sample sizes varied slightly, primarily due to an extra block of 28 °C: 60 individuals for 16 °C, 20 °C, 24 °C, 90 individuals for 28 °C, 59 individuals for 32 °C, and 61 individuals for 36 °C.

Fitting thermal performance curves (TPCs)

All analyzes were conducted in R⁷² (v4.3.1) using the ‘tidyverse’ package (v2.0.0). Figures also used ‘cowplot’ (v1.1.1), ‘patchwork’ (v1.1.2), and ‘gridExtra’ (v2.3) packages.

For each combination of trait (lifetime measures of bite rate [a], lifespan [lf], and egg production [B]) and fluctuation regime (constant, DTR 9 °C, and DTR 12 °C), we used a Bayesian framework to fit either a symmetric (quadratic: –c(T- T_min)(T- T_max)) or an asymmetric (Brière: cT(T- T_min)(T_max –T)^1/2) non-linear unimodal function to generate a TPC predicting trait values across temperature (T, in degrees Celsius). From these functions, we can compare the predicted thermal limits (T_min, T_max) and optimum temperature (T_opt) for each trait among the different DTR treatments, with c as a shape fit parameter. Both functions were restricted from becoming negative by assuming a trait value to be zero if T <T_min or T > T_max. The previous study⁸ analyzed only the constant temperature treatments and fit trait thermal responses to means from each experimental block using a truncated normal distribution. Here, we used the full dataset of three DTR treatments and fit the trait thermal responses to individual-level data, using different probability distributions for each trait based on the data type and observed distribution. For bite rate (a), we used a normal distribution truncated at zero; for lifespan (lf), we used a gamma distribution; for lifetime egg production (B), we used a negative binomial distribution (see Supplementary Methods for model specifications).

For each trait, we selected the best-fitting functional form (quadratic or Brière) using the Deviance Information Criterion (DIC)⁵⁹. For each parameter in the mean response function (i.e., c, T_min, T_max) and the additional parameter required to specify each probability distribution (i.e., the variance for the truncated normal distribution, the rate parameter for the gamma distribution, and the r parameter for the negative binomial distribution), we assumed low-information uniform priors (T_min ~ uniform (0, 20), T_max ~ uniform (28, 45), c ~ uniform (0, 10), variance ~ uniform (0,1000), rate ~ uniform (1100), r ~ uniform (1100)) that restricted the range of parameters to biologically or statistically meaningful values (i.e., values outside of those ranges are not possible). TPCs were fitted in R using JAGS⁷³ (v.4.3.2) and the ‘R2jags’ package⁷⁴ (v. 0.7.1.1), which implements Markov Chain Monte Carlo (MCMC). MCMC results were visualized with ‘mcmcplots’ (v0.4.3). Posterior draws were obtained from three concurrent Markov chains. In each chain, a 5000-step burn-in phase was followed by 20,000 samples of the stationary chain, for a total of 60,000 posterior samples. These samples were thinned by saving every eighth sample (yielding 7500 samples) to reduce autocorrelation in the chain. For each TPC, we used the posterior distributions for the parameters to generate posterior distributions over a temperature gradient from 0–45 °C at 0.1 °C intervals, which we then used to calculate the mean, median, and 95% credible intervals.

To test for the statistical significance of fluctuation treatment, we used the Deviance Information Criterion (DIC) output from JAGS. For each trait, we compared: 1) the sum of DIC values for the three models fit separately to data from each treatment (constant, DTR 9 °C, and DTR 12 °C) and 2) the DIC of a model fit to the combined data from all treatments. Fluctuation treatment is significant if the sum of the separate models is ≥ 2 DIC units lower than the DIC value for the combined model.

Generating TPCs with rate summation

To calculate the trait thermal responses predicted by rate summation (Eq. 1) we used the 7500 posterior samples from the Bayesian fitted TPCs for each trait measured at constant temperatures. First, we used a Parton-Logan model⁷¹ to calculate a temperature profile for each mean temperature spanning 0-50 °C with 0.1 °C increments, assuming a DTR of 9 or 12 °C across a 24-hour period (see Supplementary Methods). Second, we calculated predicted trait values at each hour using the TPC for trait performance at constant temperatures. Third, a daily mean value for each trait was calculated by averaging the predicted hourly values for that trait over the 24-hour period for each mean temperature. When fluctuating temperatures extended beyond the range of our constant temperature TPCs (0 °C ≥ T ≤ 45 °C), we used the trait value predicted at the corresponding edge temperature, which was always equal or approximately equal to zero. Lastly, since rate summation was conducted for each posterior sample, we calculated the mean, median, and 95% credible interval of the resulting rate summation estimates for each mean temperature. Rate summation code also used the ‘progress’ package (v1.2.3).

Predicting thermal suitability, S(T)

Following previous work⁸, we use a modified expression for the relative pathogen basic reproductive number (relative R₀), a metric of pathogen transmission potential in a given thermal environment. This metric incorporates the thermal responses of mosquito and parasite traits to evaluate the combined effects of temperature and temperature fluctuation on the predicted thermal suitability [S(T), Eq. 2] of An. stephensi to transmit Plasmodium falciparum⁸. A scaled version of R₀(T), called S(T), is proportional to the number of new cases expected to arise from a single case assuming a fully susceptible population, and is dependent on environmental temperature, T ( °C). Further, because values for mosquito life history traits change as mosquitoes age, we have adopted the use of the S(T) expression that more precisely captures lifetime transmission potential⁸ (Eq. 2).

The parameters of S(T) include: daily per capita bite rate (a), vector competence (bc; the proportion of infectious mosquitoes), lifetime egg production (B), probability of egg-to-adult survival (p_EA), mosquito development rate (MDR), and adult mosquito lifespan (lf). Further, the S(T) formulation uses the Gompertz function over daily adult survival and the extrinsic incubation period (EIP, the inverse of the parasite development rate [PDR^-1]) to calculate the proportion of mosquitoes surviving the latency period (ϒ) as described in ref. ⁸. We fit thermal responses for these additional traits (p_EA, MDR, and bc) using previously published data measured across constant temperature gradients^19,36. For ϒ, we combined data for PDR measured across a constant temperature gradient¹⁹ with our new lifespan (lf) data in constant and fluctuating conditions, and fit a TPC for each of our three fluctuation treatments (DTR = 0, 9, and 12 °C). In all cases, we used the same methods as for the focal trait data collected here (described above), with a truncated normal distribution. We calculated thermal suitability using the full posterior distributions for each trait TPC over the temperature gradient from 0-45 °C at 0.1 °C intervals, yielding posteriors for suitability over that same gradient, with the same number of samples (7500). We then used these distributions to calculate the mean, median, and 95% credible intervals.

Absolute R₀(T) is influenced by additional factors that we do not incorporate in this study including rainfall, humidity, mosquito habitat quantity and quality, infection status, and heterogeneity in contact rates, individuals, or genotypes. Thus, we instead describe the thermal suitability of pathogen transmission, S(T), where S(T) is scaled to range between 0 and 1 at the respective minimum and maximum values for the median thermal response. We scaled all versions of the S(T) model using the maximum value from model version 1 (‘constant’, see Suitability Model Overview) in order to be able to visually compare differences in the predicted magnitude of thermal suitability between model versions. The additional R₀ parameters r (human recovery rate) and N (density of humans) are evaluated as arbitrary constants, as they are assumed to be temperature independent. Thus, a threshold of S(T) > 0 implies that the thermal conditions are suitable for the transmission of P. falciparum based solely on the temperature-dependent physiological responses of An. stephensi. Differences in the predicted critical temperatures at which S(T) reaches 0 (T_min and T_max) and 1 (T_opt) can then be compared across diurnal temperature ranges.

Sensitivity and uncertainty analysis

We performed two types of sensitivity analysis and an uncertainty analysis on each version of the suitability model to determine which traits were most important for determining the thermal optimum and limits for transmission and how each trait contributed to the uncertainty in S(T). First, we used a partial derivative approach, calculating ∂S/∂x·∂x/∂T across the temperature (T) gradient for each trait (x). This approach only works for the models without rate summation (i.e., model 1: constant and model 2: empirical fluctuating) because it uses the derivatives of the quadratic and Brière functions and their fitted parameters (T_min, T_max, and q) for each trait. Second, we held each trait constant while allowing all others to vary with temperature. Finally, we calculated the HPD interval (highest posterior density interval, the smallest interval of predicted trait value encompassing 95% of the probability density in the posterior distribution) across the temperature gradient for S(T) using the full posterior distributions for all traits (i.e. full uncertainty) and for S(T) with each trait given its mean value (i.e. removing the uncertainty for one trait at a time). We then compared the relative size of the HPD in both conditions for each trait.

Mapping thermal suitability predictions

We created maps to compare the spatial distribution of months of thermal suitability for transmission predicted by the different versions of our model, S(T). For simplicity, we only mapped model versions 1, 2, 4, and 5 (constant, empirical fluctuating, trait-level RS fluctuating – all traits, and S(T)-level RS fluctuating, respectively) for one level of DTR (12 °C) where applicable. As with previous mapping for thermal suitability of transmission^5,8,23,42,58, for each version of S(T) we determined the temperature range (at 0.1 °C resolution) where S(T) > 0.001 with a posterior probability >97.5%. This conservative threshold minimizes type I error (inclusion of unsuitable areas). Here, we also calculated the temperature range at which each model exceeded an additional threshold of suitability, S(T) > 0.5. This threshold shows where thermal suitability is relatively high (rather than simply present), and allows us to illustrate quantitative differences between model versions 4 and 5 (i.e. rate summation performed on the trait TPCs versus on the suitability TPC), which had similar T_min and T_max but different shapes otherwise. For calculating both mapping thresholds, we used the 2.5% lower CI prediction from each model scaled between 0 and 1 so that relative suitability was based on the maximum predicted suitability for that specific model (rather than model 1).

Global gridded long-term average modeled baseline monthly mean temperatures at a 5 arcminute resolution (approximately 10 km2 at the equator), were downloaded from WorldClim.org (version 1.0). The number of months (0-12) of thermal suitability under each combination of model and suitability threshold was calculated at the pixel level, and masked to countries described as the ‘endemic’ range for An. stephensi (India, Pakistan, Iran, Kuwait, United Arab Emirates, and Oman), and for all countries in the continent of Africa, where it is currently invading and establishing. Raster calculations and mapping output used the additional R packages: ‘raster’ (v3.6.26), ‘terra’ (v1.7.78), ‘sf’ (v1.0.17), ‘mapproj’ (v1.2.11), ‘mapdata’ (v2.3.1), ‘ggthemes’ (v5.1.0), and ‘rnaturalearth’ (v1.0.1).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.