Temperature impacts the transmission of malaria parasites by Anopheles gambiae and Anopheles stephensi mosquitoes

Extrinsic environmental factors influence the spatio-temporal dynamics of many organisms, including insects that transmit the pathogens responsible for vector-borne diseases (VBDs). Temperature is an especially important constraint on the fitness of a wide variety of insects, as they are primarily ectotherms. Temperature constrains the distribution of ectotherms and therefore of the infections that they spread in both space and time. More concretely, a mechanistic understanding of how temperature impacts traits of ectotherms to predict the distribution of ectotherms and vector-borne infections is key to predicting the consequences of climate change on transmission of VBDs like malaria. However, the response of transmission to temperature and other drivers is complex, as thermal traits of ectotherms are typically non-linear, and they interact to determine transmission constraints. In this study, we assess and compare the effect of temperature on the transmission of two malaria parasites, Plasmodium falciparum and Plasmodium vivax, by two malaria vector species, Anopheles gambiae and Anopheles stephensi. We model the non-linear responses of temperature dependent mosquito and parasite traits (mosquito development rate, bite rate, fecundity, egg to adult survival, vector competence, mortality rate, and parasite development rate) and incorporate these traits into a suitability metric based on a model for the basic reproductive number across temperatures. Our model predicts that the optimum temperature for transmission suitability is similar for the four mosquito-parasite combinations assessed in this study. The main differences are found at the thermal limits. More specifically, we found significant differences in the upper thermal limit between parasites spread by the same mosquito (An. stephensi) and between mosquitoes carrying P. falciparum. In contrast, at the lower thermal limit the significant differences were primarily between the mosquito species that both carried the same pathogen (e.g., An. stephensi and An. gambiae both with P. falciparum). Using prevalence data from Africa and Asia, we show that the transmission suitability metric S(T) calculated from our mechanistic model is an important predictor of malaria prevalence. We mapped risk to illustrate the areas in Africa and Asia that are suitable for malaria transmission year-round based temperature.

435,000 resulted in death. The African region accounts for about 91.5% of malaria cases and species more than two thermal traits were absent. Thus we focus our analysis on An. stephensi 99 and An. gambiae as these are the most complete (See SM; Appendix A1). Since even these sets 100 have data gaps in thermal traits, where there is a gap we use traits data available from the closest 101 related species based on similar biologic and ecological characteristics (see SM; Appendix A2). (1) where a is the mosquito biting rate; bc is vector competence; µ is the mosquito mortality rate; 112 P DR is the parasite development rate; EF D is the mosquito fecundity expressed as the number 113 of eggs per female per day; P EA is the proportion of eggs to adulthood; M DR is the mosquito 114 development rate; N is the density of humans or hosts; and r is the human recovery rate. Because 115 we are interested in the shape of the thermal response only, we define a suitability metric, S(T ), 116 that only incorporates the temperature dependent components, that is: where T 0 is the lower thermal limit (where the response becomes zero), T m is the upper thermal where T 0 is the lower thermal limit, T m is the upper thermal limit, and a is a constant that 129 determines the curvature at the optimum. As with the Briére function we assume the trait is 130 piecewise zero above and below the thermal limits. For concave-up symmetric responses like µ, 131 we fit a concave-up quadratic function 132 f qu (T ) = aT 2 − bT + c where a, b, and c are the standard quadratic parameters. Note that because all traits must be ≥ 0 133 we also truncate this function, creating a piecewise continuous function where f qu is set to zero if 134 the quadratic evaluates to a negative value.

Bayesian Fitting of Thermal Traits
We fit each unimodal thermal response for all traits for each mosquito species (An. stephensi and 137 An. gambiae) or parasite strain (P. falciparum and P. vivax) with a Bayesian approach using the 138 JAGS/rjags package (Plummer, 2016) in R (R Development Core Team, 2017). We defined an 139 appropriate likelihood for each trait (e.g., binomial likelihoods for proportion data, truncated 140 normal for continuous numeric traits) with the mean defined by the either a Briére function (for 141 asymmetric relationships) or quadratic (symmetric relationships). For all traits we chose relatively 142 uninformative priors that limit each parameter to its biologically realistic range. More specifically, 143 we assumed that temperatures below 0 • C and above 45 • C are lethal for both mosquitoes and 144 malaria parasites (Lyons et al., 2012;Mordecai et al., 2017). Based on these assumptions we set 145 uniform priors for the minimum temperature (T 0 ) between 0 -24 • C and for the the maximum unimodal thermal response we calculated the posterior mean and 95% highest probability density SM; Appendix A4 for the traits used for each mosquito/parasite set). We used these samples of 162 S(T ) across temperature to calculate the posterior mean and the 95% HPD of the overall thermal 163 response of suitability, as well as for the critical thermal minimum, thermal maximum, and 164 optimal temperature for transmission suitability by the two mosquito species for each parasite. 165 We also used the posterior samples of the extrema of S(T ) (lower thermal limit, upper thermal 166 limit, and optimum) to assess the significance of the differences in these summaries between the were random draws from the sample distribution, we would expect to see approximately even 176 numbers of differences above and below zero, and the mean would be close to zero. We 177 calculated the mean of the differences (primarily to identify which of the baseline or comparative 178 was on average higher, with a positive value indicating that the baseline was larger). We then 179 calculated the proportion of differences that skewed towards the baseline. If the proportion of 180 differences in the direction of the baseline is greater than or equal 0.95, we classify this as a 181 significant difference in the statistic or there is strong evidence that the values are different. If the 182 proportion of that the differences skewed in the direction of the baseline is between 0.6 and 0.94,

183
we said that there is some evidence that the difference is significant. Otherwise we said that there 184 is not a significant difference, because the difference between these two groups could be 185 explained by chance alone. 186 We also performed an uncertainty analysis for S(T ) to quantify the contribution of each trait to 187 the overall uncertainty in suitability. We calculated the uncertainty due to each thermal trait on 188 S(T ) as follows. First we set all traits but one focal trait to their respective posterior means. We  logarithm transformed versions of both of population density and GDP. 218 We categorized each location in our dataset with a binary response variable to indicate whether or 219 not malaria had been observed at that location in our prevalence dataset, that is we defined the 220 response y i for the i th location as We fit a series of logistic regression models with different subsets of predictors that include S(T ), log population density, and log GDP, or a combination of these (See SM; Appendices A9 and A10 for full models). The general logistic GLM is defined by the mean equation: where η i is the linear predictor, with β 0 as the intercept, β 1 , . . . , β n are the regression parameters,  Smyth, 1996). 234 We then chose between all candidate models in terms of model probabilities (P r) that are based no residual deviance and its (D 2 ) will be equal to 1. The (D 2 ) is the GLM analogue to (R 2 ) in 242 linear models (Dunn and Smyth, 2018).

243
Following methodology for spatial validation of suitability prediction models in ( few overall data points, resulting in significant overall uncertainty, but lack of higher temperature 274 data leads to especially large uncertainty in the upper thermal curve/limit. Bite rate (a) similarly, 275 has more uncertainty at the upper thermal limit, although the effect is less pronounced because 276 more total data are available.

277
In contrast, most An. stephensi traits exhibit considerably less uncertainty around the mean. An. stephensi with P. falciparum is more uncertain at the minimum thermal limit because data 283 were only available above 21 • C ( Figure 2A). Similarly PDR for P. falciparum is uncertain near 284 the maximum thermal limit because there are only three data points above 32 • C ( Figure 2C). In 285 contrast P. vivax related traits are better constrained. .7 • C respectively ( Figure 3). However, these median ranges mask a great deal of uncertainty.

298
The suitability metric for all four mosquito parasite pairs is predicted to peak (i. Appendices A7, B12c, and B15c). The other comparisons at the maximum thermal limit are not 308 significant (See SM; Appendices B12 and B15).

309
The predictions for the lower thermal limits were much more similar to each other. Our results

310
show that the only significant differences are between the transmission of P. falciparum by  Appendices B5 and B6). This is a common pattern as S ∝ µ −3 , so small changes in µ when µ is 318 small (that is near optimal temperatures for mosquito longevity) will have an out-sized impact on 319 S. This component is thus almost completely responsible for the location and height of the peak 320 of suitability.

321
In contrast, the traits that drive uncertainty in the temperatures around the thermal limits varies 322 between each mosquito-parasite pair and is sensitive to the amount and quality of data available for each. For example, our suitability metric for P. falciparum in An. stephensi seems the most 324 well-resolved of all of the combinations. In this case, uncertainty in µ dominates across almost all 325 temperatures, and it is only at the high temperature end, above ≈ 37 • C, that most of the 326 uncertainty is caused by something else, specifically egg to adult survival (P EA

Model Validation
We assessed whether our suitability metric is consistent with patterns spatially explicit prevalence 347 data for P. falciparum in Africa and Asia and P. vivax in Asia in two ways. For each set of parasite 348 prevalence data we evaluated seven models that included a combination of the logarithm of 349 population density, logarithm of GDP, and the probability that the suitability metric is greater than 350 zero (i.e., S(T ) > 0 , SGTZ). This included three "null" models that incorporated socio-economic 351 factors but not the suitability metric. We chose as the best model the one that has the lowest value 352 of the Bayesian Information Criterion (BIC), and also evaluated the relative model probability 353 based on BIC. We also used a visual, qualitative approach, that examines histograms of the This indicates that the suitability metrics are also likely consistent with the observed pattern of 384 prevalence in both places, although further data to confirm this would be helpful. Our maps indicate that in Africa, approximately 15% of the continental land area will have 395 temperatures suitable year-round for the transmission of P. falciparum malaria by An. stephensi 396 mosquitoes, and 44% for at least six months of the year ( Figure 4A). We found that An. stephensi, using data from historical and newly published studies. In doing so, we are able to examine the extent to which predictions may vary between mosquito species and between parasite 422 types. We also identified persistent data gaps that must be addressed to further improve these 423 models and allow more precise comparisons between mosquito/parasite complexes.

424
Our results suggest that there is little difference in the optimal temperature for malaria

438
In contrast, we find that there are differences between the temperatures that could limit suitability 439 for the two focal parasites to be spread by different mosquitoes. For example we found evidence 440 that there are differences in the upper thermal limit between P. vivax and P. falciparum when 441 spread by An. stephensi, and between P. falciparum when spread by the two mosquitoes. There is 442 also evidence that the lower suitability threshold differs by mosquito species when transmitting 443 the same pathogen. However, there is a great deal of uncertainty in the estimates of these limits 444 due to poor data availability, especially for traits of An. gambiae and for P. falciparum.

445
As a result of the observed difference in the upper and lower thermal limits across 446 mosquito-parasite systems, we also observe differences in the predicted temperature ranges for suitability, with greater temperature ranges for pathogens transmitted by An. stephensi than 448 An. gambiae (Figure 3). Based on available trait data and previous studies, An. stephensi 449 mosquitoes seem to have greater plasticity in thermal tolerance than An. gambiae mosquitoes 450 (Bayoh, 2001;Kirby and Lindsay, 2004;Miazgowicz et al., 2020). Our predictions largely agree 451 with these broader mosquito thermal studies.

452
The larger thermal breadth for transmission of both P. falciparum and P. vivax by An. stephensi 453 than An. gambiae has a knock-on effect for the spatial extent of suitability for transmission. Our  An. stephensi with P. falciparum. These data would improve certainty in these models, especially 490 at the thermal limits.

491
Our approach has some important limitations, some of which could be addressed by extending the 492 mechanistic models. One limitation is that we use constant temperatures in our models, and did

Model
Formula for linear predictor BIC Model pr. Table A9: Logistic regression models used for validation exercise -to test if the probability S(T )>0 alone or in combination with socio-economic variables is a good predictor of malaria prevalence caused by P. falciparum in Africa and Asia. For Africa and Asia we have 8343 and 3001 prevalence records respectively with 7934 P. falciparum positive prevalence cases and 3410 P. falciparum negative prevalence cases. The logistic model is defined so that for each response, y i , P r(y i = 1) = θ = logit −1 (η i ) where η i is the linear predictor as given in the table above.
Observations at data point i are then Binomial random variables with "success" probability θ and sample size n i . Abbreviations used on the   Figure B4: Comparison of the thermal response for parasites traits between mosquito/parasite systems: An. gambiae/P. falciparum (An. gamb P. falc), An. gambiae/P. vivax (An. gamb P. vivax), An. stephensi/P. falciparum (An. step P. falc), and An. stephensi/P. vivax (An. step P. vivax). A) Parasite development rate and B) Vector competence. Min: Minimum temperature, Opt: optimum temperature, and Max: Maximum temperature.  Figure B8: Plots of the probability that S(T ) >0 versus the probability of malaria prevalence by Plasmodium falciparum transmitted by An. gambiae in Africa and by An. stephensi in Asia from presence/absence model Pfalc PA5 (Table 8; SM), for different levels of population density across different rows and different levels of per capita Gross Domestic Product (GDP) across different columns. Tick-marks are human malaria prevalence data. Lines and shaded areas: mean and 95% confidence interval from the fitted GLM. Figure B9: The natural log of the probability of S(T )>0 times the per capita Gross Domestic Product (GDP) predict the probability of prevalence of P. vivax transmitted by An. stephensi in Asia. Figure B10: Plots of the probability that S(T ) >0 versus the probability of malaria prevalence by P. vivax transmitted by An. stephensi in Asia from presence/absence model Pviv PA5 (Table 9; SM), for different levels of population density across different rows and different levels of GDP across different columns. Tick-marks are human malaria prevalence data. Lines and shaded areas: mean and 95% CI from GLM. Figure B11: Randomized quantile residuals extracted in R using the qresid function in the package statmod (Smyth et al., 2019) for the best fitted models (Pfal PA5 and Pviv PA4) shown in tables 8 and 9 in SM. Randomized Quantile Residuals are interpreted as standards residuals, and should be normally distributed if the assumptions of the underlying model are appropriate for the data. A) Fitted values plotted versus residuals for P. falciparum presence (red) and absence (black) prevalence cases. B) Q-Q plot for the quantile residuals for P. falciparum. C) fitted values plotted versus residuals for P. vivax presence (red) and absence (black) of prevalence cases. D) Q-Q plot for the quantile residuals for P. vivax.

B.3 Comparison of the thermal responses for parasite and mosquito traits
B.6 Summary of the sample posterior distributions, S(T ) Figure B12: Histograms of the S(T ) posterior distributions at the maximum temperature limit. A) An. stephensi with P. falciparum versus An. stephensi with P. vivax. B) An. gambiae with P. falciparum and An. gambiae with P. vivax. C) An. gambiae with P. falciparum versus An. stephensi with P. falciparum. D) An. gambiae with P. vivax versus An. stephensi with P. vivax. Figure B13: Histograms of the S(T ) posterior distributions at the optimum temperature limit. A) An. stephensi with P. falciparum versus An. stephensi with P. vivax. B) An. gambiae with P. falciparum and An. gambiae with P. vivax. C) An. gambiae with P. falciparum versus An. stephensi with P. falciparum. D) An. gambiae with P. vivax versus An. stephensi with P. vivax. Figure B14: Histograms of the S(T ) posterior distributions at the minimum temperature limit. A) An. stephensi with P. falciparum versus An. stephensi with P. vivax. B) An. gambiae with P. falciparum and An. gambiae with P. vivax. C) An. gambiae with P. falciparum versus An. stephensi with P. falciparum. D) An. gambiae with P. vivax versus An. stephensi with P. vivax.
B.7 Comparison of the mosquito/parasite posterior distributions, S(T ) Figure B15: Sample differences of the S(T ) posterior distributions at the maximum temperature limit. A) An. stephensi with P. falciparum versus An. stephensi with P. vivax B) An. gambiae with P. falciparum versus An. gambiae with P. vivax C) An. stephensi with P. falciparum versus An. stephensi with P. falciparum, and D) An. stephensi with P. falciparum versus An. gambiae with P. vivax. Figure B16: Sample differences of the S(T ) posterior distributions at the optimum temperature limit. A) An. stephensi with P. falciparum versus An. stephensi with P. vivax B) An. gambiae with P. falciparum versus An. gambiae with P. vivax C) An. stephensi with P. falciparum versus An. gambiae with P. falciparum, and D) An. stephensi with P. vivax versus An. gambiae with P. vivax. Figure B17: Sample differences of the S(T ) posterior distributions at the minimum temperature limit. A) An. stephensi with P. falciparum versus An. stephensi with P. vivax B) An. gambiae with P. falciparum versus An. gambiae with P. vivax C) An. stephensi with P. falciparum versus An. gambiae with P. falciparum, and D) An. stephensi with P. vivax versus An. gambiae with P. vivax. B.9 Proportion of suitable months for malaria transmission in Africa and Asia Figure B19: The proportion of number of suitable months for transmission S(T )>0 in Africa with the presence of A) P. falciparum malaria prevalence in area suitable for An. stephensi B) P. vivax malaria prevalence in area suitable for An. stephensi C) P. falciparum malaria prevalence in area suitable for An. gambiae D) P. vivax malaria prevalence in area suitable for An. gambiae.