Transmission of West Nile and five other temperate mosquito-borne viruses peaks at temperatures between 23°C and 26°C

1Temperature is a key driver of transmission of mosquito-borne diseases. Both mosquitoes and the pathogens they transmit are ectotherms whose physiology and life histories depend strongly on environmental temperature. These temperature-dependent traits drive the biological processes required for transmission. For example, temperature-dependent fecundity, development, and mortality of mosquitoes determine whether vectors are present in sufficient numbers for transmission. Temperature also affects the mosquito biting rate on hosts and probability of becoming infectious.

2Mechanistic models based on these traits and guided by principles of thermal biology predict that the thermal response of transmission is unimodal: transmission peaks at intermediate temperatures and declines at extreme cold and hot temperatures. This unimodal response is predicted consistently across mosquito-borne diseases and supported by independent empirical evidence for positive relationships between temperature and human cases in many settings, but negative relationships at high temperatures in other studies. Accordingly, we expect increasing temperatures due to climate change to shift disease distributions geographically and seasonally, as warming increases transmission in cooler settings but decreases it in settings near or above the optimal temperature for transmission. Thus, mechanistic models have provided a powerful and general rule describing how temperature affects the transmission of mosquito-borne disease. However, thermal responses vary among mosquito and pathogen species and drive important differences in how predicted transmission responds to temperature, including the specific temperatures of the optimum and thermal limits for each vector–pathogen pair. We currently lack a framework to describe or predict this variation among vectors and pathogens.

3Filling this gap requires comparing mechanistic, temperature-dependent transmission models for many vector–pathogen pairs. However, models that incorporate all relevant traits are not yet available for many important pairs for several reasons. First, the number of relevant vector–pathogen pairs is large because many mosquitoes transmit multiple pathogens and many pathogens are transmitted by multiple vectors. Second, empirical data are costly to produce, and existing data are often insufficient because experiments or data reporting were not designed for this purpose. Here, we address these challenges by systematically compiling data and building models for understudied mosquito-borne disease systems, including important pathogens with substantial transmission in temperate areas like West Nile virus (WNV) and Eastern Equine Encephalitis virus (EEEV). Accurately characterizing the thermal limits and optima for these systems is critical for understanding where and when temperature currently promotes or suppresses transmission and where and when climate change will increase, decrease, or have minimal effects on transmission.

4In this study, we model the effects of temperature on an overlapping suite of widespread, important mosquito vectors and viruses that currently lack complete temperature-dependent models. These viruses include: West Nile virus (WNV), St. Louis Encephalitis virus (SLEV), Eastern and Western Equine Encephalitis viruses (EEEV and WEEV), Sindbis virus (SINV), and Rift Valley fever virus (RVFV) (;;;;; summarized in Table 1). All but RVFV sustain substantial transmission in temperate regions. We selected this group because many of the viruses share common vector species and several vector species transmit multiple viruses (Table 1, Figure 1). All the viruses cause febrile illness and severe disease symptoms, including long-term arthralgia and neuroinvasive syndromes with a substantial risk of mortality in severe cases. Since invading North America in 1999, WNV is now distributed worldwide and is the most common mosquito-borne disease in the US, Canada, and Europe. SLEV, EEEV, and WEEV occur in the Western hemisphere (Table 1), with cases in North, Central, and South America. For EEEV, North American strains are genetically distinct and more virulent than the Central and South American strains. An unusually large outbreak of EEEV in the United States in 2019 has yielded incidence four times higher than average (31 cases, resulting in nine fatalities) and brought renewed attention to this disease. SINV occurs across Europe, Africa, Asia, and Australia, with substantial transmission in northern Europe and southern Africa. RVFV originated in eastern Africa and now also occurs across Africa and the Middle East. These pathogens primarily circulate and amplify in wild bird reservoir hosts (except RVFV, which primarily circulates in livestock). For all six viruses, humans are dead-end or unimportant reservoir hosts, in contrast to pathogens like malaria, dengue virus, yellow fever virus, and Ross River virus, which sustain infection cycles between humans and mosquitoes. Most transmission of RVFV to humans occurs through direct contact with infected livestock (that are infected by mosquitoes), and to a lesser extent via the mosquito-borne transmission from infected vectors.

5We primarily focus on Culex pipiens, Cx. quinquefasciatus, and Cx. tarsalis, well-studied species that are important vectors for many of the viruses and for which appropriate temperature-dependent data exist for nearly all traits relevant to transmission. Although the closely-related Cx. pipiens and Cx. quinquefasciatus overlap in their home ranges in Africa, they have expanded into distinct regions globally (Figure 2; ). Cx. pipiens occurs in higher latitude temperate areas in the Northern and Southern hemisphere, while Cx. quinquefasciatus occurs in lower latitude temperate and tropical areas (Figure 2A). By contrast, Cx. tarsalis is limited to North America but spans the tropical-temperate gradient (Figure 2B). In this system of shared pathogens and vectors with distinct geographical distributions, we also test the hypothesis that differences in thermal performance underlie variation in vector and pathogen geographic distributions, since temperate environments have cooler temperatures and a broader range of temperatures than tropical environments. We also include thermal responses from other relevant vector or laboratory model species in some models: Aedes taeniorhynchus (SINV and RVFV), Ae. triseriatus (EEEV), Ae. vexans (RVFV), Cx. theileri (RVFV), and Culiseta melanura (EEEV). Additionally, we compare our results to previously published models for transmission of more tropical diseases by the following vectors: Ae. aegypti, Ae. albopictus, Anopheles spp., and Cx. annulirostris.

6We use a mechanistic approach to characterize the effects of temperature on vector–virus pairs in this network using the thermal responses of traits that drive transmission. Specifically, we use experimental data to measure the thermal responses of the following traits: vector survival, biting rate, fecundity, development rate, competence for acquiring and transmitting each virus, and the extrinsic incubation rate of the virus within the vector. We ask: (1) Do these vectors have qualitatively similar trait thermal responses to each other, and to vectors from previous studies? (2) Is transmission of disease by these vectors predicted to be optimized and limited at similar temperatures, compared to each other and to other mosquito-borne diseases in previous studies? (3) How do the thermal responses of transmission vary across vectors that transmit the same virus and across viruses that share a vector? (4) Which traits limit transmission at low, intermediate, and high temperatures? Broadly, we hypothesize that variation in thermal responses is predictable based on vectors’ and viruses’ geographic ranges.

7Mechanistic models allow us to incorporate nonlinear effects of temperature on multiple traits, measured in controlled laboratory experiments across a wide thermal gradient, to understand their combined effect on disease transmission. This approach is critical when making predictions for future climate regimes because thermal responses are almost always nonlinear, and therefore current temperature–transmission relationships may not extend into temperatures beyond those currently observed in the field. We use Bayesian inference to quantify uncertainty and to rigorously incorporate prior knowledge of mosquito thermal physiology to constrain uncertainty when data are sparse. The mechanistic modeling approach also provides an independently generated, a priori prediction for the relationship between temperature and transmission to test with observational field data on human cases, allowing us to connect data across scales, from individual-level laboratory experiments, to population-level patterns of disease transmission, to climate-driven geographic variation across populations. Using this approach, we build mechanistic models for 10 vector–virus pairs by estimating thermal responses of the traits that drive transmission. We validate the models using observations of human cases in the US over space (county-level) and time (month-of-onset). The validation focuses on WNV because it is the most common of the diseases we investigated and has the most complete temperature-dependent trait data. Preliminary results of this study—the thermal responses for traits and relative R0 models—were included in a review and synthesis article that was published last year. The present publication presents the complete methods and results, describes the vector and pathogen ecology in more detail, and provides original analyses of human case data.

8To understand the effect of temperature on transmission and to compare the responses across vector and virus species, we used R0—the basic reproduction number. We use R0 as a static, relative metric of temperature suitability for transmission that incorporates the nonlinear effects of constant temperature on multiple traits and is comparable across systems, rather than focusing on its more traditional interpretation as a threshold for disease invasion into a susceptible population. Temperature variation creates additional nonlinear effects on transmission that are not well-captured by R0, but could be incorporated in future work by integrating the thermal performance curves fit here over the observed temperature regime.

9The basic R0 model for mosquito-borne diseases, originally developed for malaria (Equation 1; ), includes the following traits that depend on temperature (T): adult mosquito mortality rate (µ, the inverse of lifespan [lf]), biting rate (a, the inverse of the gonotrophic [oviposition] cycle duration), pathogen development rate (PDR, the inverse of the extrinsic incubation period: the time required for exposed mosquitoes to become infectious), and vector competence (bc, the proportion of exposed mosquitoes that become infectious), where all rates are measured in inverse days. Vector competence is the product of infection efficiency (c, the proportion of exposed mosquitoes that develop a disseminated infection) and transmission efficiency (b, the proportion of infected mosquitoes that become infectious, with virus present in saliva). Mosquito density (M) also depends on temperature but is not an organism-level trait that can be measured in a laboratory setting. Two parameters do not depend on temperature: host density (N) and the rate at which infected hosts recover and are no longer infectious (r).

10Because host density (N) and recovery rate (r) are not temperature-dependent, we omit them from our model (Equation 2), which isolates the effect of temperature on transmission (see explanation of ‘relative R0’ versus absolute R0 below). As in previous work, we extend the basic R0 model to account for the effects of temperature on mosquito density (M) via additional temperature-sensitive life history traits (Equation 2): fecundity (as eggs per female per day, EFD), egg viability (proportion of eggs hatching into larvae, EV), proportion of larvae surviving to adulthood (pLA), and mosquito development rate (MDR, the inverse of the development period).

11Fecundity data were only available as eggs per female per gonotrophic cycle (EFGC; for Cx. pipiens) or eggs per raft (ER; for Cx. quinquefasciatus). Thus, we further modified the model to obtain the appropriate units for fecundity: we added an additional biting rate (a) term to the numerator (to divide by the length of the gonotrophic cycle, Equations A1 and A2) and for Cx. quinquefasciatus we also added a term for the proportion of females ovipositing (pO; Equation A2).

12We parameterized the full temperature-dependent R0 model (Equation 2) for each relevant vector–virus pair using previously published data. We conducted a literature survey to identify studies that measured the focal traits at three or more constant temperatures in a controlled laboratory experiment. From these data, we fit thermal responses for each trait using Bayesian inference. This approach allowed us to quantify uncertainty and formally incorporate prior data to constrain fits when data for the focal species were sparse or only measured on a limited portion of the temperature range (see Material and Methods for details).

13For each combination of trait and species, we selected the most appropriate of three functional forms for the thermal response. As in previous work, we fit traits with a symmetrical unimodal thermal response with a quadratic function (Equation 3) and traits with an asymmetrical unimodal thermal response with a Briére function (; Equation 4). For some asymmetrical responses (e.g. pathogen development rate [PDR] for most vector–virus pairs), we did not directly observe a decrease in trait values at high temperatures due to a limited temperature range. In these cases, we chose to fit a Briére function based on previous studies with wider temperature ranges and thermal biology theory; the upper thermal limit for these fits did not limit transmission in the R0 models, and therefore did not impact the results. Unlike in previous work, lifespan data for all vectors here exhibited a monotonically decreasing thermal response over the range of experimental temperatures available. We fit these data using a piecewise linear function (Equation 5) that plateaued at the coldest observed data point. By assuming a plateau, rather than extrapolating that lifespan continues to increase at temperatures below those measured in the laboratory, this approach is conservative, ensuring that lifespan was not a major driver of the temperature-dependence of R0 at temperatures where it was not measured and that the R0 models were instead constrained at reasonable temperatures by other traits. It is also consistent with the observed natural history of two of the vector species. To overwinter, Cx. pipiens and Cx. tarsalis enter reproductive diapause and hibernate, and Cx. pipiens can survive temperatures at or near freezing (0°C) for several months. Cx. quinquefasciatus enters a non-diapause quiescent state and is likely less tolerant of cold stress, but we wanted a consistent approach across models and other traits constrained the lower thermal limit of the Cx. quinquefasciatus R0 model to realistic temperatures. All vectors are likely to exhibit decreased lifespans at extremely low temperatures (near or below 0°C), limiting the accuracy of our inferred lifespan thermal performance curve at these temperatures.

14In the quadratic and Briére functions of temperature (T), the trait values depend on a lower thermal limit (Tmin), an upper thermal limit (Tmax), and a scaling coefficient (q). In the linear function, the trait values depend on a slope (m) and intercept (z).

15The fitting via Bayesian inference produced posterior distributions for each parameter in the thermal response functions (Equations 3, 4, 5) for each trait–species combination. These posterior distributions represent the estimated uncertainty in the parameters. We used these parameter distributions to calculate distributions of expected mean thermal performance functions for each trait over a temperature gradient (from 1°C to 45°C by 0.1°C increments). Then we substituted these samples from the distributions of the thermal responses for each trait into Equation 2 to calculate the posterior distributions of predicted R0 over this same temperature gradient for each vector–virus pair (see Material and methods and Appendix 1 for details). Thus, the estimated uncertainty in the thermal response of each trait is propagated through to R0 and combined to produce the estimated response of R0 to temperature, including the uncertainty in R0(T).

16Because the magnitude of realized R0 depends on system-specific factors like breeding habitat availability, reservoir and human host availability, vector control, species interactions, and additional climate factors, we focused on the relative relationship between R0 and temperature. We rescaled the R0 model results to range from 0 to 1 (i.e. ‘relative R0’), preserving the temperature-dependence (including the absolute thermal limits and thermal optima) while making each model span the same scale. To compare trait responses and R0 models, we quantify three key temperature values: the optimal temperature for transmission (Topt) and the lower and upper thermal limits (Tmin and Tmax, respectively) where temperature is predicted to prohibit transmission (R0 = 0).

17We fit thermal response functions from empirical data for all of the vector and virus traits that affect transmission for which data were available (Figure 1 and Appendix 1—table 1). All mosquito traits were temperature-sensitive (three main Culex species: Figure 3, Figure 4; Ae. taeniorhynchus, Ae. triseriatus, Ae. vexans, Cx. theileri, and Cs. melanura: Appendix 1—figure 1). For most species, the extensive data for larval traits (mosquito development rate [MDR] and survival [pLA]) produced clear unimodal thermal responses with relatively low uncertainty (Figure 3A,B, Appendix 1—figure 1A,B). For biting rate (a) and fecundity traits (proportion ovipositing [pO], eggs per female per gonotrophic cycle [EFGC], or per raft [ER], and egg viability [EV]), trait data were often more limited and fits were more uncertain, but still consistent with the expected unimodal thermal responses based on previous studies and theory (; Figure 3C, Figure 4, Appendix 1—figure 1C-F). However, adult lifespan (lf) data clearly contrasted with expectations from previous studies of more tropical mosquitoes. Lifespan (lf) decreased linearly over the entire temperature range of available data (coldest treatments: 14–16°C, Figure 3D; 22°C, Appendix 1—figure 1D) instead of peaking at intermediate temperatures (e.g. previously published optima for more tropical species: 22.2–23.4°C).

18The thermal responses for pathogen development rate (PDR) were similar among most vector–virus pairs (Figure 5), with a few notable exceptions: WNV in Cx. quinquefasciatus had a warmer lower thermal limit (Figure 5A); WNV in Cx. univittatus had a cooler optimum and upper thermal limit (Figure 5A); and SINV in Ae. taeniorhynchus had limited data that indicated very little response to temperature (Figure 5C). By contrast, the thermal response of vector competence (bc) and its component traits varied substantially across vectors and viruses (Figure 6). For example, infection efficiency (c) of Cx. pipiens peaked at warmer temperatures for WNV than for SINV (Figure 6A,G; 95% CIs: SINV = 14.1–30.5°C, WNV = 31.9–36.1°C), transmission efficiency (b) of Cx. tarsalis peaked at warmer temperatures for WNV and SLEV than for WEEV (Figure 6B,E,H; CIs: WEEV = 19.2–23.2°C, SLEV = 23.5–29.7°C, WNV = 23.9–29.3°C), and the lower thermal limit for vector competence (bc) for WNV was much warmer in Cx. pipiens than in Cx. univittatus (Figure 6C; CIs: Cx. univittatus = 1.5–7.1°C, Cx. pipiens = 15.0–17.9°C). Infection data (used to calculate pathogen development rate [PDR] and vector competence [bc]) for RVFV and SINV were only available in Ae. taeniorhynchus, a New World species that is not a known vector for these viruses in nature.

19Relative R0 responded unimodally to temperature for all the vector–virus pairs, with many peaking at fairly cool temperatures (medians: 22.7–26.0°C, see Table 2 for CIs; Figure 7). The lower thermal limits (medians: 8.7–19.0°C, see Table 2 for CIs; Figure 7) were more variable than the optima or the upper thermal limits (medians: 31.9–37.8°C, see Table 2 for CIs; Figure 7), although confidence intervals overlapped in most cases because lower thermal limits also had higher uncertainty (Figure 7). The Ae. taeniorhynchus models (both unnatural vector-pathogen pairs) were clear outliers, with much warmer distributions for the upper thermal limits, and optima that trended warmer as well.

20Differences in relative R0 stemmed from variation both in vector traits (e.g. in Figure 7A, with WNV in different vector species) and in virus infection traits (e.g. in Figure 7B, with different viruses in Cx. tarsalis). The upper thermal limit was warmer for WNV transmitted by Cx. pipiens (34.9°C [CI: 32.9–37.5°C]) than by Cx. quinquefasciatus (31.8°C [CI: 31.1–32.2°C]), counter to the a priori prediction based on the higher-latitude range of Cx. pipiens in North America, South America, and Europe (Figure 2). This result implies that warming from climate change may differentially impact transmission by these two vectors. Additionally, the lower thermal limit for WNV varied widely (but with slightly overlapping 95% CIs) across different vector species (Figure 7D), from 19.0°C (14.2–21.0°C) in Cx. quinquefasciatus to 16.8°C (14.9–17.8°C) in Cx. pipiens to 12.2°C (9.7–15.3°C) in Cx. tarsalis to 11.1°C (8.1–15.4°C) in Cx. univittatus (an African and Eurasian vector; Table 2). Based on these trends in the thermal limits of R0, the seasonality of transmission and the upper latitudinal and elevational limits could vary for WNV transmitted by these different species.

21Different traits determined the lower and upper thermal limits and optimum for transmission across vector–virus pairs. The lower thermal limit for transmission was most often determined by pathogen development rate (PDR; WNV and SLEV in Cx. tarsalis, WNV in Cx. quinquefasciatus) or biting rate (a; WNV in Cx. univitattus, WEEV in Cx. tarsalis, EEEV in Ae. triseriatus, RVFV and SINV in Ae. taeniorhynchus, SINV in Cx. pipiens; Appendix 1—figures 12–20), which tend to respond asymmetrically to temperature, with high optima and low performance at low temperatures. However, vector competence (bc) determined the lower limit for WNV in Cx. pipiens (Appendix 1—figure 11). The upper thermal limit was determined by biting rate (a) for the three Cx. tarsalis models and by adult lifespan (lf) for all others, although proportion ovipositing (pO) was also important for WNV in Cx. quinquefasciatus (Appendix 1—figures 11–20). In all models, lifespan (lf) and biting rate (a) had the strongest impact on the optimal temperature for transmission, with biting rate increasing transmission at low temperatures and lifespan decreasing transmission at high temperatures (Appendix 1—figures 11–20). This result is consistent with previous mechanistic models of tropical mosquito-borne diseases, despite the qualitative difference in the shape of the lifespan thermal response between those tropical mosquitoes and the more temperate mosquitoes investigated here.

22We validated the R0 models for WNV with independent data on human cases because the temperature-dependent trait data for those models were relatively high quality and because human case data were available from the Centers for Disease Control and Prevention across a wide climatic gradient in the contiguous United States. We averaged county-level incidence and mean summer temperatures across the entire period from 2001 to 2016 to estimate the impact of temperature over space, while ignoring interannual variation in disease that is largely driven by changes in host immunity and drought. We used generalized additive models (GAMs, which produce flexible, smoothed responses) to ask: does average incidence respond unimodally to mean summer temperature? If so, what is the estimated optimal temperature for transmission? Can we detect upper or lower thermal limits for transmission? Incidence of human neuroinvasive West Nile disease responded unimodally to average summer temperature and peaked at 24°C (23.5–24.2°C depending on the spline settings; Figure 8, Appendix 1—figure 24), closely matching the optima from the mechanistic models for the three North American Culex species (23.9–25.2°C; Table 2). However, the human disease data did not show evidence for lower or upper thermal limits: mean incidence remained positive and with relatively flat slopes below ~19°C and above ~28°C, although sample size was very low above 28°C and below 15°C resulting in wide confidence intervals (Figure 8, Appendix 1—figure 24).

23We used national month-of-onset data for WNV, EEEV, and SLEV to ask: is the seasonality of incidence consistent with our models for temperature-dependent transmission? The month-of-onset for cases of WNV was consistent with predicted transmission, R0(T) (Figure 9). As expected (based on previous studies and the time required for mosquito populations to increase, become infectious, and bite humans, and for humans to present symptoms and seek medical care), there was a 2-month lag between initial increases in R0(T) and incidence: cases began rising in June to the peak in August. The dramatic decline in transmission between September and October corresponds closely to the predicted decline in relative R0, but without the expected two-month lag. In general, the seasonal patterns of SLEV and EEEV incidence were similar to WNV, but differed by three orders of magnitude from ~20,000 cases of WNV to ~40–50 cases of EEEV and SLEV during the peak month (Figure 9). However, transmission of SLEV and EEEV are predicted to begin increasing 1 month earlier than WNV (March versus April, Figure 9), because the mechanistic models predict that the lower thermal limits for SLEV and EEEV are cooler than those for WNV in two of the three North American vectors (Cx. pipiens and Cx. quinquefasciatus, Figure 7). The month-of-onset data partially support this prediction, as cases of SLEV (but not EEEV) disease begin to increase earlier in the year than WNV, relative to the summer peak.

24As the climate changes, it is critical to understand how changes in temperature will affect the transmission of mosquito-borne diseases. Toward this goal, we developed temperature-dependent, mechanistic transmission models for 10 vector–virus pairs. The viruses—West Nile virus (WNV), St. Louis Encephalitis virus (SLEV), Eastern and Western Equine Encephalitis viruses (EEEV and WEEV), Sindbis virus (SINV), and Rift Valley fever virus (RVFV)—sustain substantial transmission in temperate areas (except RVFV), and are transmitted by shared vector species, including Cx. pipiens, Cx. quinquefasciatus, and Cx. tarsalis (except EEEV; Figure 1). Although most traits responded unimodally to temperature, as expected, lifespan decreased linearly with temperature over the entire temperature range of available data (>14°C) for these Culex vectors (Figure 3). Transmission responded unimodally to temperature, with the thermal limits and optima for transmission varying among some of the focal mosquito and virus species (Figure 7, Table 2), largely due to differences in the thermal responses of mosquito biting rate, lifespan, vector competence, and pathogen development rate. Human case data for WNV disease across the US exhibited a strong unimodal thermal response (Figure 8), and month-of-onset data for WNV, SLEV, and EEEV were largely consistent with the predicted seasonality of transmission (Figure 9). Thus, the mechanistic models captured geographical and seasonal patterns of human incidence, despite the complexity of the enzootic cycles and spillover into humans. Our analysis was somewhat limited by the lack of data for several trait-species combinations, or by data that were sparse, particularly at high temperatures. However, our key results—maximal transmission at intermediate temperatures—are unlikely to change, and underscore the importance of considering unimodal thermal responses when predicting how climate change will impact mosquito-borne disease transmission.

25The monotonically decreasing thermal responses of lifespan (lf) within the range of the available experimental data for these more temperate mosquitoes (Figure 3D) contrast with the clearly unimodal responses of more tropical species. This contrast may reflect differing thermal physiology between species that use diapause or quiescence, two forms of dormancy, to persist over winter and those that do not. Both Cx. pipiens and Cx. tarsalis diapause, and Cx. pipiens can survive temperatures at or near freezing (0°C) for several months. Cx. quinquefasciatus enters a non-diapause quiescent state. Ae. albopictus, a species that occurs in both tropical and temperate zones, exhibits a latitudinal gradient in the United States in which more temperate populations diapause while sub-tropical populations do not. Experiments could test this hypothesis by measuring whether the functional form of the thermal response for lifespan differs between northern (diapausing) and southern (non-diapausing) US Ae. albopictus populations. Despite the difference in the shape of the thermal response, lifespan played a similarly important role here as in previous studies of mosquito-borne pathogens, strongly limiting transmission at high temperatures (Appendix 1—figures 11–20). Nonetheless, the thermal responses for lifespan here ultimately promote higher transmission at relatively cool temperatures because unlike in more tropical species, lifespan did not decline at cool temperatures within the range measured (>14°C). Given the lack of rigorous trait data, we cannot be certain of the shape of the thermal response of lifespan (lf) below 14°C, although it is almost certainly unimodal, especially at more extreme temperatures expected to be fatal even for diapausing mosquitoes (i.e. below 0°C). Our decision to assume lifespan (lf) plateaued at temperatures below the observed data was based on vector natural history and intended to be conservative. This approach ensured that lifespan was not a major driver of the temperature-dependence of R0 at temperatures where it was not measured and that R0 was instead constrained by other traits. Accordingly, our functions for lifespan (lf) do not represent the real quantitative thermal responses below the coldest observations, which limits their utility for other applications, such as predicting survival at cold temperatures and lower thermal limits on survival.

26Predicted transmission for many of the diseases in this study peaked at and extended to cooler temperatures than for previously studied diseases with more tropical distributions (see Figure 7 and Table 2 for 95% credible intervals). Here, the optimal temperatures for transmission varied from 22.7–25.2°C (excluding Ae. taeniorhynchus models, Figure 7). By contrast, models predict that transmission peaks at 25.4°C for malaria, 26.4°C for Ross River virus and dengue in Ae. albopictus, 28.9°C for Zika in Ae. aegypti, and 29.1°C for dengue in Ae. aegypti. Models for several vector–virus pairs also had cooler lower thermal limits (medians: 8.7–19.0°C) than those of diseases with more tropical distributions (medians: 16.0–17.8°C). In combination with similar upper thermal limits (see below), these patterns led to wider thermal breadths (18.2–27.7°C; Figure 7) for most of the viruses here compared to the more tropical pathogens (11.7–16.7°C), excepting WNV in Cx. quinquefasciatus (12.7°C), the vector most restricted to lower latitude, sub-tropical geographic areas (Figure 2). These results match a previous finding that temperate insects had wider thermal breadths than tropical insects, and may reflect thermal adaptation to greater variation in temperature in temperate areas compared to tropical areas. Additionally, SINV—a virus with substantial transmission at very high latitudes in Finland—had the second coolest lower thermal limit (Figure 7, Table 2). Further, lower latitude Cx. quinquefasciatus outperformed higher latitude Cx. pipiens at warmer temperatures for proportion ovipositing (pO; Figure 3A), while the reverse occurred at cooler temperatures for egg viability (EV; Figure 3C). Collectively, these results imply that, to some extent, measurements of physiological traits can predict geographic patterns of vectors or disease transmission at broad scales. However, geographic range differences (Figure 2) did not consistently predict variation in thermal responses among the Culex species in this study (e.g. biting rate [a, Figure 3C] and adult lifespan [lf, Figure 3D]), indicating that life history and transmission trait responses at constant temperatures do not always predict the geographic distributions of species. Instead, the ability to tolerate temperature extremes may limit species distributions more than their performance at average or constant temperatures. Moreover, although diseases like malaria and dengue are generally considered to be ‘tropical’, historically their distributions extended further into temperate regions. Thus, current distributions of disease may reflect a realized niche restricted by societal factors more than a fundamental niche based on ecological factors like temperature.

27In contrast to the optima, lower thermal limits, and thermal breadths, the upper thermal limits for the vector–virus pairs in this study (31.9–34.9°C, excluding Ae. taeniorhynchus models; Figure 7D, Table 2) closely matched those of more tropical diseases (31.5–34.7°C). This similarity may arise because maximum summer temperatures in temperate areas can match or even exceed maximum temperatures in tropical areas. Accordingly, there may be a fundamental upper thermal constraint on transmission that applies similarly to all mosquito-borne diseases, driven by short mosquito lifespans at high temperatures. The relatively high upper thermal limits in both Ae. taeniorynchus transmission models were driven by the thermal response of lifespan (lf), which was fit to few data points; more data are needed to determine whether it reflects the true thermal response in that species (Appendix 1—figure 1). These results indicate that as temperatures rise due to climate change, temperate diseases are unlikely to be displaced by warming alone, although they may also expand toward the poles, even as tropical diseases may expand farther into temperate zones.

28Independent human case data support unimodal thermal responses for transmission and the importance of temperature in shaping geographic patterns of mosquito-borne disease. Human cases of WNV and SINV are often positively associated with temperature. Here, we found that incidence of neuroinvasive WNV disease peaked at intermediate mean summer temperatures (24°C) across counties in the US (Figure 8) that matched the optima predicted by our models. This result adds to prior evidence for reduced transmission of WNV and other mosquito-borne diseases at high temperatures. Although we did not detect lower or upper thermal limits for West Nile neuroinvasive disease (Figure 8), this result is unsurprising based on fundamental differences between the types of temperature data used to parameterize and validate the models. The R0 model prediction is derived from data collected in a controlled laboratory environment at constant temperatures, while average incidence in the field reflects temperatures that vary at a variety of temporal scales (daily, seasonal, and interannual). Thus, we hypothesize that temperature variation over time may sustain transmission in regions with otherwise unsuitable mean summer temperatures by providing time windows that are suitable for transmission.

29The temperature-dependent models also predict the general seasonality of human cases of WNV, EEEV, and SLEV (Figure 9). The 2-month lag between climate suitability and the onset of human cases, which matches previous results from other mosquito-borne diseases, arises from the time following the onset of suitable conditions required for mosquito populations to increase, become infectious, and bite humans, and for humans to present symptoms and seek medical care. Transmission of the more temperate viruses here may incur additional lags because human cases result from enzootic transmission, and multiple rounds of amplification within reservoir hosts may be required before prevalence is sufficiently high to spill over into humans. Additionally, as wild birds begin to migrate in late summer, both Cx. pipiens and Cx. tarsalis shift their feeding preferences from birds to humans, which should increase transmission to people later in the year. However, we found that cases decreased more quickly in autumn than expected from temperature effects alone. Human behavior may partially compensate for the shift in feeding preference and explain why the decrease of cases in autumn did not show the expected 2-month lag from temperature-dependent relative R0. For instance, if people wear clothing that exposes less skin and spend less time outdoors due to school schedules and changing daylight it may reduce contact with mosquitoes. Drought, precipitation, and reservoir and human immunity also strongly drive transmission of WNV and may interact with temperature. SLEV, EEEV, and WEEV are less common in nature, and thus less well-studied, but the lower thermal limits in our study support previous findings that transmission of WEEV is favored over SLEV in cooler conditions. Additionally, the seasonal patterns of incidence data (Figure 9) provide some support for the model prediction that SLEV transmission is possible at cooler temperatures than WNV by North American vectors (Table 2). By contrast, mean temperature is not associated with outbreaks of RVFV, although they are highly predictable based on precipitation driven by El Niño–Southern Oscillation cycles. Thus, disease dynamics depend on the interaction between temperature and other environmental factors, and the relative importance of temperature versus other drivers varies across systems.

30Most prior studies with mechanistic models for temperature-dependent transmission of WNV do not capture the unimodal thermal response that our mechanistic models predict and that we observe in the human case data (Table 3). Two previous models predicted that transmission of WNV would increase up to the warmest temperatures they considered, 28°C and 35°C. In both cases, the vector daily survival rates estimated from lab experiments were far less sensitive to temperature than our measure of adult lifespan (lf), and neither model was validated with field data. A third study with models for Cx. pipiens, Cx. quiquefasciatus, and Cx. tarsalis, like our study, predicted unimodal thermal responses for transmission, with very similar optima but with lower thermal limits that were ~5°C warmer, resulting in much narrower thermal breadths (Appendix 1—figure 23; ). This previous set of models was validated with annual, state-level WNV human case data (in contrast to our county-level data averaged over multiple years), and detected a positive effect of temperature, with no decline at high temperatures. The best spatial and temporal scales for validating temperature-dependent transmission models and detecting the impacts of temperature remain an open question. For instance, different approaches may be necessary to detect thermal optima and thermal limits. Critically, differences in modeling and validation approaches can lead to strongly divergent conclusions and predictions for the impact of climate change.

31Given the unimodal relationship between temperature and transmission of these temperate mosquito-borne pathogens, we expect climate warming to lead to predictable shifts in disease transmission. Warming should extend the transmission season earlier into the spring and later into the fall and increase transmission potential in higher latitudes and altitudes, although this prediction may be impacted by changes in bird migrations. However, the thermal optima for these temperate vector–virus pairs are relatively cool, so in many locations, warming could result in summer temperatures that exceed the thermal optima for transmission more frequently, reducing overall transmission or creating a bimodal transmission season. Based on the average summer temperature data (2001–2016) in our analysis (Figure 8), currently the majority of people (70%) and counties (68%) are below the optimal temperature for transmission (23.9°C, fit by the GAM). The numbers are similar when restricted to counties with observed West Nile virus cases: 69% and 70%, respectively. Thus, all else being equal, we might expect a net increase in transmission of West Nile virus in response to the warming climate, even as hot temperatures suppress transmission in some places. Still, warming is unlikely to eliminate any of these more temperate pathogens since the upper thermal limits for transmission are well above temperatures pathogens regularly experience in their current geographic ranges. More generally, our results raise concerns about the common practice of extrapolating monotonic relationships between temperature and disease incidence fit from observational data into warmer climate regimes to predict future cases.

32While the data-driven models presented here represent the most comprehensive synthesis to date of trait thermal response data and their impact on transmission for these mosquito–pathogen systems with substantial transmission in temperate regions, additional temperature-dependent trait data would increase the accuracy and decrease the uncertainty in these models where data were sparse or missing. Our data synthesis and uncertainty analysis suggest prioritizing pathogen development rate (PDR) and vector competence (bc) data and biting rate (a) data because those thermal responses varied widely among vector–virus pairs and determined the lower thermal limits and optima for transmission in many models. Additionally, vector competence (bc) and/or pathogen development rate (PDR) data were missing in many cases (WNV in Cx. quinquefasciatus and Cx. modestus [an important vector in Europe], EEEV in Cs. melanura, RVFV in vectors from endemic areas, transmission efficiency [b] for SINV) or sparse (EEEV and WNV in Cx. univittatus), as were biting rate data (Cx. univittatus, RVFV vectors). Lifespan (lf) data—key for determining transmission optima and upper thermal limits—were the missing for Ae. triseriatus, Cs. melanura, Cx. univittatus, and RVFV vectors, and at temperatures below 14°C for all vector species, so it was unclear which functional form these thermal responses should take (monotonic, saturating, or unimodal). While the other mosquito demographic traits did not determine thermal limits for transmission in models here, fecundity (typically as eggs per female per day, EFD), larval-to-adult survival (pLA), and egg viability (EV) determined thermal limits for malaria and Ross River virus. Thus, more fecundity data (missing for Cx. tarsalis, Cx. univittatus, and Ae. triseriatus; sparse for Cx. pipiens and Cx. quinquefasciatus) would also increase our confidence in the models. New data are particularly important for RVFV: the virus has a primarily tropical distribution in Africa and the Middle East, but the model depends on traits measured in Cx. pipiens collected from temperate regions and infection traits measured in Ae. taeniorhynchus, a North American species. This substitution of a mosquito species that is not a naturally occurring vector could reduce the relevance and utility of this model. RVFV is transmitted by a diverse community of vectors across the African continent, but experiments should prioritize hypothesized primary vectors (e.g. Ae. circumluteolus or Ae. mcintoshi) or secondary vectors that already have partial trait data (e.g. Ae. vexans or Cx. theileri). Although temperature itself does not predict the occurrence of RVFV outbreaks, it may affect the size of epidemics once they are triggered by precipitation. More generally, thermal responses may vary across vector populations and/or virus isolates even within the same species. Several studies have found differences in thermal performance across different populations of the same mosquito species or pathogen strains, but this variation was not systematically associated with their thermal environments of origin. Accordingly, the potential for thermal adaption in mosquitoes and their pathogens remains an open question. Regardless, more data may improve the accuracy of all the models, even those without missing data.

33Our trait-based R0 models effectively isolated the physiological effects of temperature on transmission. However, in nature many other environmental and biological factors also impact transmission of mosquito-borne disease. For example, potential factors include rainfall, habitat and land-use, reservoir host community composition, host immunity, viral and mosquito genotypes, mosquito microbiome, vector control efforts, vector behavior, and human behavior. Our analyses here suggest that temperature is important for shaping broad-scale spatial and seasonal patterns of disease when cases are averaged over time and space. Other factors may be more important at finer spatial and temporal scales, and explain additional variation in human cases. For instance, a study of WNV and two other (non-mosquito-borne) pathogens found that biotic factors were significant drivers of disease distributions at local scales, while climate factors were only significant drivers at larger regional scales. Given that our R0 models for WNV predicted very similar thermal optima across three distantly-related vector species, it is likely that our results are generalizable to other temperate locations with the same vectors (e.g. parts of Europe with transmission by Cx. pipiens) at similarly broad spatial and temporal scales, even if the other factors influencing local-scale patterns are quite different than in the US.

34As carbon emissions continue to increase and severe climate change becomes increasingly inevitable, it is critical that we understand how temperature change will affect the transmission of mosquito-borne diseases in a warmer future world. While data gaps are still limiting, the mechanistic, trait-based approach presented here is powerful for predicting similarities and differences across vectors and viruses and for making predictions for the impact of climate change. Accounting for the effects of temperature variation is an important next step for using these types of models to accurately predict transmission. In nature, mosquitoes and pathogens experience daily temperature variation that can dramatically alter performance compared to constant temperatures of the same mean. Rate summation is the most common method for predicting performance in variable temperatures based on experimental data at constant temperatures. This approach is ideal because mean temperature and daily temperature variation vary somewhat independently over space and time, and measuring vector and pathogen performance at sufficient combinations of both is logistically difficult. However, its accuracy for predicting mosquito and pathogen traits or mosquito-borne disease transmission has not been rigorously evaluated. Additionally, the potential for adaptive evolution to warmer climates is uncertain because of limited knowledge on the level of genetic variation in thermal responses for vectors or their pathogens within or between populations. Further, vectors and pathogens may experience different selective pressures, as mosquito populations may depend on either increased fecundity or longevity at high temperatures, while pathogens require longer vector lifespans. Thus, future trajectories of these diseases will depend not just on suitability of mean temperatures but also on temperature variation, thermal adaptation of vectors and viruses, land use (which governs mosquito–wildlife–human interactions), vector control activities, human and wildlife immune dynamics, and potential future emergence and spread of new vectors and viruses.

35All analyses were conducted using R 3.1.3.

36The distributions of Cx. pipiens and Cx. quinquefasciatus are georectified maps adapted from;. The northern boundary of Cx. tarsalis was taken from. For the southern boundary, we drew a convex polygon using five datasets in the Global Biodiversity Information Facility (https://www.gbif.org/).

37We found 38 studies with appropriate temperature-dependent trait data from controlled laboratory experiments. When necessary, we digitized the data using Web Plot Digitizer, a free online tool. When lifespan data were reported by sex, only female data were used. Vector competence trait data (b, c, or bc) were only included if time at sampling surpassed the estimated extrinsic incubation period (the inverse of PDR) at that temperature, which resulted in the exclusion of some studies.

38We fit trait thermal responses with a Bayesian approach using the ‘r2jags’ package, an R interface for the popular JAGS program for the analysis of Bayesian graphical models using Gibbs sampling. It is a (near) clone of BUGS (Bayesian inference Using Gibbs Sampling). In JAGS, samples from a target distribution are obtained via Markov Chain Monte Carlo (MCMC). More specifically, JAGS uses a Metropolis-within-Gibbs approach, with an Adaptive Rejection Metropolis sampler used at each Gibbs step (for more information on MCMC algorithms see ).

39For each thermal response being fit to trait data, we visually identified the most appropriate functional form (quadratic, Briére, or linear; Equations 3–5) for that specific trait–species combination. For traits with ambiguous functional responses, we fit the quadratic and Briere and used the deviance information criterion (DIC) to pick the best fit. We assumed normal likelihood distributions with temperature-dependent mean values described by the appropriate function (Equations 3–5) and a constant standard deviation (σ) described by an additional fitted parameter (τ = 1/σ2). The 95% credible intervals in Figures 3–6 estimate the uncertainty in the mean thermal response.

40We set all thermal response functions to zero when T < Tmin and T > Tmax (for Equation 3 and 4) or when T > -z/m (Equation 5) to prevent trait values from becoming negative. For traits that were proportions or probabilities, we also limited the thermal response functions at 1. For the linear thermal responses, we calculated the predicted thermal response in a piecewise manner in order to be conservative: for temperatures at or above the coldest observed data point, we used the trait values predicted by the fitted thermal response (i.e. the typical method); for temperatures below the coldest observed data point, we substituted the trait estimate at the coldest observed data point (i.e. forcing the thermal response to plateau, rather than continue increasing beyond the range of observed data).

41For the fitting process, we ran three concurrent MCMC chains for 25,000 iterations each, discarding the first 5000 iterations for burn-in (convergence was checked visually). We thinned the resultant chains, saving every eighth step. These settings resulted in 7500 samples in the full posterior distribution that we kept for further analysis.

42We used data-informed priors to decrease the uncertainty in our estimated thermal responses and constrain the fitted thermal responses to be biologically plausible, particularly when data were sparse. These priors used our total dataset, which contained temperature-dependent trait data for all of the main species in the analysis (but with the focal species removed, see below), as well as from additional temperate Aedes and Culex species.

43We fit each thermal response with a sequential two-step process, where both steps employed the same general fitting method (described above in Fitting Thermal Responses) but used different priors and data. In step 1, we generated high-information priors by fitting a thermal response to data from all species except the focal species of interest (i.e. a ‘leave-one-out’ approach). For example, for the prior for biting rate for Cx. pipiens, we used the biting rate data for all species except Cx. pipiens. For this step, we set general, low-information priors that represented minimal biological constrains on these functions (e.g. typically mosquitoes die if temperatures exceed 45°C, so all biological processes are expected to cease; Tmin must be less than Tmax). The bounds of these uniformly distributed priors were: 0 < Tmin < 24, 26 < Tmax < 45 (quadratic) or 28 < Tmax < 45 (Briére), 0 < q < 1,–10 < m < 10, and 0 < b < 250. Then in step 2, we fit a thermal response to data from the focal species using the high-information priors from step 1.

44Because we cannot directly pass posterior samples from JAGS as a prior, we modified the results from step 1 to use them in step 2. We used the ‘MASS’ package to fit a gamma probability distribution to the posterior distributions for each thermal response parameter (Tmin, Tmax, and q [Equation 3 and 4]; or m and z [Equation 5]) obtained in step 1. The resulting gamma distribution parameters can be used directly to specify the priors in the JAGS model. Because the prior datasets were often very large, in many cases the priors were too strong and overdetermined the fit to the focal data. In a few other cases, we had philosophical reasons to strongly constrain the fit to the focal data even when they were sparse (e.g. to constrain Tmax to very high temperatures so that other traits with more information determine the upper thermal limit for R0). Thus, we deflated or inflated the variance as needed (i.e., we fixed the gamma distribution mean but altered the variance by adjusting the parameters that describe the distribution accordingly). See Appendix 1 for more details and specific variance modifications for each thermal response.

45When data were missing for a vector–virus pair, we used two criteria to decide which thermal response to use as a substitute: 1) the ecological similarly (i.e. geographic range overlap) of species with available thermal responses, and 2) how restrictive the upper and lower bounds of the available thermal responses were. All else being equal, we chose the more conservative (i.e. least restrictive) option so that R0 would be less likely to be determined by trait thermal responses that did not originate from the focal species. See Appendix 1 for more information about specific models.

46When there was more than one option for how to parameterize a model (e.g. vector competence data for WEEV in Cx. tarsalis were available in two forms: separately as b and c, and combined as bc), we calculated R0 both ways. The results were very similar, except for the model for RVFV with lifespan data from Cx. pipiens lifespan in place of Ae. taeniorhynchus (Appendix 1—figure 22). See Appendix 1 for sensitivity and uncertainty methods and Appendix 1—figures 11–20 for results.

47We obtained county-level neuroinvasive WNV disease data from 2001 to 2016 for the contiguous US (n = 3109) through the CDC’s county-level disease monitoring program. Data were available as total human cases per year, which we adjusted to average cases per 1000 people (using 2010 US county-level census data) to account for population differences. We averaged cases across years beginning with the first year that had reported cases in a given county to account for the initial spread of WNV and the strong impact of immunity on interannual variation. Ninety-eight percent of human cases of WNV in the US occur between June and October (data described below), and cases of mosquito-borne disease often lag behind temperature by 1–2 months. Thus, we extracted monthly mean temperature data between the months of May–September for all years between 2001 and 2016 and averaged the data to estimate typical summer conditions for each county. Specifically, we took the centroid geographic coordinate for every county in the contiguous US with the ‘rgeos’ package and extracted corresponding historic climate data for monthly mean temperatures (Climate Research Unit 3.1 rasters) from 0.5°2 cells (approximately 2500–3000 km2) using the ‘raster’ package. The monthly mean temperatures in this climate product are calculated by averaging daily mean temperatures at the station level (based on 4–8 observations per day at regular intervals) and interpolating these over a grid.

48We fit a generalized additive model (GAM) for average incidence as a function of average summer temperature using the ‘mgcv’ package. We used a gamma distribution with a log-link function to restrict incidence to positive values and capture heteroskedasticity in the data (i.e. higher variance with higher predicted means), adding a small, near-zero constant (0.0001) to all incidence values to allow the log-transformation for counties with zero incidence. GAMs use additive functions of smooth predictor effects to fit responses that are extremely flexible in the shape of the response. We restricted the number of knots to minimize overfitting (k = 7; see Appendix 1—figure 24 for results across varying values of k). For comparison, we also used the ‘loess’ function in base R ‘stats’ package to fit locally estimated scatterplot smoothing (LOESS) regressions of the same data. LOESS regression is a simpler but similarly flexible method for estimating the central tendency of data. See Appendix 1—figure 25 for LOESS model results. See Appendix 1—figure 26 for non-binned county-level data.

49We calculated monthly temperature-dependent relative R0 to compare with month-of-onset data for neuroinvasive WNV, EEEV, and SLEV disease aggregated nationwide (the only spatial scale available) from 2001 to 2016, using the same county-level monthly mean temperature data as above. For WNV, we used the subset of counties with reported cases (68% of counties). For SLEV and EEEV we used all counties from states with reported cases (16 and 20 states, respectively). We calculated a monthly R0(T) for each county, and then weighted each county R0(T) by its population size to calculate a national monthly estimate of R0(T). For WNV, the county-level estimates of R0(T) used models for three Culex species (Cx. pipiens, Cx. quinquefasciatus, and Cx. tarsalis) weighted according to the proportion of WNV-positive mosquitoes reported at the state level, reported in. SLEV and EEEV both only had one R0 model. The estimated monthly temperature-dependent relative R0 values and month-of-onset data were compared visually.

50All data and code are available on Github in the following repository: https://github.com/mshocket/Six-Viruses-Temp (; copy archived at https://github.com/elifesciences-publications/Six-Viruses-Temp). All data and code are also available in the Dryad Data Repository.