Role of vector phenotypic plasticity in disease transmission as illustrated by the spread of dengue virus by Aedes albopictus

Domain: earth_science
Published: 2024-09-07
Editorial outcome: major_revision

Abstract
The incidence of vector-borne disease is on the rise globally, with burdens increasing in endemic countries and outbreaks occurring in new locations. Effective mitigation and intervention strategies require models that accurately predict both spatial and temporal changes in disease dynamics, but this remains challenging due to the complex and interactive relationships between environmental variation and the vector traits that govern the transmission of vector-borne diseases. Predictions of disease risk in the literature typically assume that vector traits vary instantaneously and independently of population density, and therefore do not capture the delayed response of these same traits to past biotic and abiotic environments. We argue here that to produce accurate predictions of disease risk it is necessary to account for environmentally driven and delayed instances of phenotypic plasticity. To show this, we develop a stage and phenotypically structured model for the invasive mosquito vector, Aedes albopictus, and dengue, the second most prevalent human vector-borne disease worldwide. We find that environmental variation drives a dynamic phenotypic structure in the mosquito population, which accurately predicts global patterns of mosquito trait-abundance dynamics. In turn, this interacts with disease transmission to capture historic dengue outbreaks. By comparing the model to a suite of simpler models, we reveal that it is the delayed phenotypic structure that is critical for accurate prediction. Consequently, the incorporation of vector trait relationships into transmission models is critical to improvement of early warning systems that inform mitigation and control strategies.

Introduction
Vector-borne diseases (VBDs) are primarily vectored by ectothermic arthropods, whose life history is sensitive to environmental variation1. Our ability to predict if vector populations can sustain pathogen transmission across the species range requires an understanding of how environmental variation and vector life-history interact2,3. There is now a rich literature exploring the mechanisms through which environmental variation alters vector trait expression, but explicit and delayed mechanisms of individual variation are generally omitted, even in extensively studied systems4. For example, adult mosquitoes experiencing hot, dry summer conditions have a shorter lifespan than adults that are subject to more favourable temperatures5. Short-lived mosquitoes are less likely to survive to complete the extrinsic incubation period and take a subsequent blood meal to transmit the pathogen. Therefore, there is environmentally driven inter-annual and inter-regional variation in the ability of mosquitoes to vector disease6.

Predictions of the relative risk of VBDs, via metrics such as the basic reproduction number, are highly sensitive to traits such as adult longevity. Seasonal and regional variation in this trait is often accounted for by parameterising functions that directly relate the current temperature to the longevity of adult mosquitoes7,8. However, adult longevity has also been shown to vary in response to an individual’s experience of competition and temperature during development. Therefore, the traits expressed by individuals in the population result from complex interactions between past environmental conditions and population states9. By ignoring the role that a vector’s historic environmental experiences have on biological traits critical to disease transmission we implicitly make the mean-field assumption, that all individuals are equally competent vectors regardless of their past experiences of the biotic and abiotic environment10. This mean-field assumption compromises our ability to assess the risk of vector-borne disease in variable environments, where seasonal patterns in environmental variables induce seasonal variation in vector population dynamics. Similarly, mean-field assumptions limit our ability to assess the relative risk of VBDs between climatic zones, such as between tropical and temperate environments, where different patterns of environmental variation occur.

The delayed effect of developmental experience on adult traits is an example of phenotypic plasticity, the ability of organisms expressing the same genotype to express different phenotypes according to the environmental conditions they are subject to11. To demonstrate how mechanisms of delayed phenotypic plasticity can be incorporated into predictions of VBDs we consider the ability of the invasive mosquito species Aedes albopictus to vector dengue virus. We use a recently derived phenotypically explicit mathematical modelling framework to represent how historical biotic and abiotic environmental conditions determine the trait dynamics of vector populations12. Ae. albopictus is a competent vector of dengue virus that is now widely distributed in temperate zones into which dengue is regularly imported13. Dengue is a viral VBD that has seen a recent dramatic increase in cases: in 2000 there were about half a million cases, which has risen to over 4.2 million in 2022, caused by a combination of climate and anthropogenic change14.

Despite predictions of broad suitability for the transmission of dengue virus by Ae. albopictus in temperate climates such as Europe by both statistical and mechanistic modelling approaches, outbreaks in these regions have so far been limited in size and duration15,16. This mismatch between observed and predicted disease incidence motivates the development of modelling approaches that are better able to reflect the currently observed global incidence of dengue. Due to the species' global range (currently between 0 and 52.5 degrees latitude13), Ae. albopictus populations are subject to a diverse range of environmental conditions. Environmental variation has been shown to induce variable trait dynamics in field populations and this may alter the ability of populations in different environments to vector disease17,18. It is our hypothesis that, by accounting for the effect of biotic and abiotic environmental variation on the phenotypic trait structure and population dynamics of Ae. albopictus, both the population dynamics of the species across its range and the current patterns of disease incidence around the globe can be better understood. Further, we propose that omitting mechanisms of phenotypic plasticity limits the ability of models to produce accurate and generalisable predictions of disease incidence.

Results
Model overview and validation

We develop a Susceptible-Exposed-Infected-Resistant (SEIR) model for the transmission of dengue virus by Ae. albopictus that incorporates a system of environmentally driven stage and phenotypically structured delay-differential equations. The model explicitly predicts mosquito population and trait dynamics. This model can represent the instantaneous effects of the current environment (e.g., temperature, precipitation, evaporation, photoperiod, and larval density) on vector traits such as development rate and through-stage survival. In addition, it captures the delayed effect of historic biotic and abiotic environments on these same traits. This allows us to represent how phenotypic plasticity in response to historic biotic and abiotic environmental variation during development alters adult longevity. We explicitly track the number of adult mosquitoes within the population that have experienced specific developmental conditions. This results in a phenotypically structured adult population that directly accounts for the effects of past environments on the ability of mosquitoes to transmit VBDs.

We use this model to consider the infection dynamics of a dengue outbreak, beginning with the introduction of humans, infected with a single serotype of dengue virus, into a completely susceptible population. This simplified representation of the urban dengue virus transmission cycle is most applicable to non-endemic temperate and sub-tropical regions, where the likelihood of multiple serotypes circulating simultaneously is reduced and the population’s prior exposure to dengue virus is limited19. Our model incorporates detailed species-specific processes, such as diapause, quiescence, the water dynamics of developmental habitats, and developmental plasticity that causes adult traits to vary in response to the environmental conditions experienced by juveniles throughout development20,21. We use data from published laboratory experiments to parametrise multi-dimensional reaction norms describing how temperature, larval density, precipitation dynamics, and photoperiod alter the traits of mosquitoes within and across developmental stages. Remotely sensed and climate reanalysis data combined with human population density data are then used to make spatio-temporal predictions of vector dynamics and VBD risk across the globe. A full description of the model is provided in the Methods section and a detailed example of the model outputs can be found in Supplementary Note 1.

We extensively validate the model predictions of mosquito population, wing length, and disease dynamics by comparing them to field surveys of mosquito populations and disease outbreaks across the species range. To validate the population dynamics we use published surveillance data, including life-stage specific population density estimates and average trait data, from 40 locations, across 13 countries and 4 continents (see Supplementary Note 1 for the full set of validations). To validate the model predictions of dengue dynamics we compare them to historical dengue outbreaks, using reports from the outbreak to select a likely introduction scenario for dengue virus into the region and to define an area over which to simulate the model. Validation data was obtained from a comprehensive search of published studies that observed the population dynamics of Ae. albopictus in the field or reported human cases of dengue. The literature search was conducted using a snowball procedure - a systematic search revealed a high number of false positives. Studies were included in our validation if the data was collected with at least a monthly resolution from a region where Aedes aegypti, a closely related vector species was absent. This exclusion criteria is necessary as Ae. aegypti is known to compete with Ae. albopictus for larval resources and is also a vector of dengue virus, factors which we do not account for here. Since substantial variation in case reporting effort through space and time is likely, our predictions of dengue transmission are not rescaled and the validations are intended to demonstrate that the model produces plausible disease dynamics (see also Supplementary Note 2). The comparisons presented in this paper represent the full set of datasets meeting the above criteria that were found in our literature search.

Across the species global range, the model achieves excellent predictions of the population and trait dynamics of Ae. albopictus (Fig. 1). The model’s predictions hold over a broad range of climatic regions and reflect the differences in population dynamics that are observed between temperate and tropical environments (see Monmouth, New Jersey, USA (2009), R2=0.59, and St. Paul, La Réunion (2013), R2=0.65 (Supplementary Fig. 47A)22,23. The model accurately predicts detailed within-season dynamics and consistently captures the species’ phenology (see adult numbers in Guangzhou, China (2006–2015), R2=0.6 (Supplementary Fig. 45B), and Rimini, Italy (2008), R2=0.94 (Supplementary Fig. 9A)24,25. The model additionally captures the inter-annual variation in abundance observed in multi-year datasets (see oviposition activity in Trento, Italy (2010–2020), R2=0.71 (Supplementary Fig. 14), and adult numbers in Suffolk, USA (2009–2018), R2=0.64 (Supplementary Fig. 40))26,27.

Our model predictions indicate and confirm the current limited ability of Ae. albopictus to sustain autochthonous transmission in the Alpes-Maritimes Department of France. Predictions show good resemblance to the disease dynamics observed during the multi-year epidemic on the island of Réunion and the large outbreak that occurred in Guangzhou, China (Fig. 2)28–33. The reduced resemblance of our predictions to the observed disease dynamics towards the latter half of each outbreak may be explained by the intensification of vector control activities which we do not account for here. For example, in Guangzhou in 2014, a considerable intensification of vector control activities occurred halfway through the outbreak to which, according to previous disease modelling studies, led to a substantial reduction in the final outbreak size and duration28. These validations demonstrate that, through the careful consideration of the mechanisms by which environmental variation acts on the life history of a vector species, broadly generalisable predictions of population and disease dynamics can be achieved.

The role of phenotypic plasticity in driving disease outbreaks

To examine how phenotypic plasticity contributes to the observed disease dynamics, we examine four of the outbreaks in Fig. 2 in greater detail. During each outbreak, we compare the wing lengths of all mosquitoes within the modelled population to the wing lengths of infected mosquitoes responsible for dengue virus transmission events (Fig. 3). Wing length is a trait that is observable in the field where it is known to seasonally vary in a manner that qualitatively matches with the variation we have observed in populations of Ae. albopictus and other mosquito species17,34–36. This variation is highly correlated with variation in adult longevity9. Adult longevity is used in the parametrisation of the model to link larval developmental experience to adult trait expression9. Large mosquitoes are produced as a result of low larval competition and cooler temperatures at the start of the adult season. These large mosquitoes have long lifespans which allows a greater proportion of individuals to survive through the extrinsic incubation period of the virus which is also prolonged when temperatures are low.

As each outbreak progresses the initial cohort of large mosquitoes begins to die off, being replaced by subsequent cohorts of smaller individuals that developed under higher intraspecific competition and warmer temperatures. Consequently, the wing length distribution of infected mosquitoes becomes increasingly similar to that of the entire population (Fig. 3A–D). This demonstrates that the trait dynamics of uninfected mosquitoes and mosquitoes responsible for infections are distinct, with the relationship between the two evolving through time and between outbreak locations. By considering differences between these wing length distributions we are able to directly examine the effects of phenotypic plasticity in adult longevity on disease dynamics9. For example, in Guangzhou during first quarter of the outbreak, an average of 70% of infections transmitted each day are attributable to the largest 20% of mosquitoes present in the population over this period (equivalently this is 47%, 31%, and 25% of transmission, attributable to the largest 20% of mosquitoes in Alpes-Maritimes, Reunion, and Hawai’i, respectively). More generally, across all of the outbreaks considered here we find that, in the initial stages of each outbreak, the majority of dengue virus transmission can be attributed to mosquitoes with wing lengths that are higher than the population average (Fig. 3E).

In the transmission of dengue, extended longevity is the only advantage that we have assumed that large mosquitoes have over small mosquitoes. The disproportionate amount of disease transmission events caused by large individuals during the initial stage of each outbreak is therefore attributable to the effects of delayed environmental variation on the trait structure of vector populations. This result adds significant nuance to the seminal work of Macdonald37 on malaria transmission, which demonstrated that under equilibrium conditions, adult longevity is an important factor in determining whether a VBD can spread. We show that temporal variation in adult longevity is essential for amplifying pathogen transmission in the early stages of an outbreak.

To test our hypothesis that vector trait dynamics are integral to predicting disease dynamics, we develop a suite of model variants (simplifications of our full model) that make common assumptions about how vector traits respond to environmental variation (see Supplementary Note 3). In the full model, we temporally track the number of mosquitoes within the population that express each wing length, resulting in a temporally evolving trait distribution, induced by delayed phenotypic plasticity. We consider a variant of our full model, where instead all adults are assumed to share the same invariant wing length, an assumption that is common in derivatives of the classical Ross-Macdonald transmission model, from which epidemiological metrics such as the basic reproduction number are routinely derived. We refer to this as the constant wing length model and consider its dynamics for different constant wing lengths4. Under this assumption, the only variation in adult traits occurs due to changes in the instantaneous environmental conditions, and so the effects of previous environments on the expression of traits are omitted. We also consider a model variant that uses the same delayed mechanisms of trait expression as the full model, but that represent this variation through a population average wing length rather than tracking the temporal evolution of the full trait distribution. We refer to this model variant as the unstructured model. Under this assumption, phenotypic plasticity is represented through a mean-field, where although the population's average past experience of development is accounted for, there is no explicit trait structure.

When the constant wing length model is simulated for the outbreak on Reunion (selected as a point of comparison due to its multi-year duration), under the assumption that all mosquitoes express the average wing length predicted by the full model (2.6 mm), we find that the constant wing length model predicts half as many total infections over the duration of the outbreak compared to the full model (see Supplementary Fig. 59). This demonstrates that the inclusion of mechanisms of individual variation not only allows us to attribute a disproportionate number of transmission events to large mosquitoes (Fig. 3), but also shows that the vector population’s trait dynamics have a quantifiable effect on outbreak size, and therefore that phenotypic plasticity alters disease dynamics (see Supplementary Fig. 56 for examples of the constant wing length model for other outbreaks). For the same outbreak, the unstructured model variant predicts twice as many dengue cases as the full model. The mean-field approach, characterised by the unstructured model, does not account for the non-linear effect of adult survival on vector competence, resulting in quantitatively different predictions of outbreak size (Supplementary Fig. 59). These examples show how the delayed response of vector traits to past environments through phenotypic plasticity gives rise to a trait distribution that is integral to understanding the disease dynamics we observe. Further, they show how these disease dynamics cannot be anticipated by mean-field approaches, even when they incorporate all the same mechanisms of delayed phenotypic plasticity.

Global risk

To predict transmission risk over the species global range, we derive an expression for the reproduction number, Rt, which describes the number of infections produced by the introduction of a single infected human into the population (described in the Methods and compared to the simpler model variants in Supplementary Note 3). This formulation of Rt accounts for the effect of phenotypic plasticity in the vector population trait structure, on the ability of that population to transmit disease. It produces predictions that are generalisable between climates (see Methods for details). We compute Rt across the species' global range and report the number of consecutive months per year between 2010−2020 that Rt > 1, in areas with a mean annual relative humidity greater than 55%, conditions that indicate the suitability of a region for the transmission of dengue virus and the survival of adult mosquitoes. In general, we predict less suitability for the autochthonous transmission of dengue vectored by Ae. albopictus than approaches that do not account for the delayed effects of phenotypic plasticity (See Supplementary Note 4 for a comparison between our predictions and those made by a standard formulation of the basic reproduction number R0 that uses similar environmental drivers and reaction norms). Furthermore, by comparing the average period that Rt > 1 to the average density of adult female mosquitoes during the active season (Fig. 4A, C, E) we see that abundance alone does not predict the ability of vector populations to transmit disease. Additional outputs of parameters of interest, such as the time of peak adult density and the first and last days for which transmission is predicted can be found in Supplementary Note 5.

In North America and Europe, regions that currently experience limited autochthonous transmission of dengue virus by Ae. albopictus, our approach predicts that there is currently limited transmission risk from this vector. By comparison, in areas of China where dengue outbreaks vectored by Ae. albopictus are more frequent, our approach predicts longer transmission periods over a wider area, demonstrating that our predictions reflect currently observed regional differences in dengue incidence38. By comparing the locations of historical dengue outbreaks to locations that are predicted to be at risk, we observe that despite our predictions being highly specific, they still encompass locations known to experience autochthonous dengue virus transmission33,39,40. For example, the location of 38 of the 41 clusters of autochthonous dengue transmission in Europe reported by the European Centre for Disease Control since 2010 occurred within regions that we predict are at risk of dengue transmission (Fig. 4D)13. This suggests that over the last 10 years populations of Ae. albopictus in temperate regions have had limited capability to sustain dengue transmission due to environmental constraints despite the species now wide distribution in these regions. Given the changing nature of the climate, understanding when this widely distributed species will be capable of sustaining prolonged transmission cycles in temperate regions is of critical importance that is yet to be resolved, and which is the subject of subsequent modelling work.

Methods
To predict both the trait and population dynamics of Ae. albopictus we apply the model framework developed in Brass et al.12. This framework combines a stage-structured population model of the form described in Nisbet and Gurney (1983) with an additional phenotypic structure that is able to represent the effects of within and between generation plasticity72. The phenotypic structure describes how an individual’s previous experience of its environment during development alters the traits it currently expresses. Individuals with different experiences of the past environment are grouped into different environmental classes which describe how this experience alters the traits they express now. By linking individual-level variation to population response, we are able to predict not only the population’s dynamics but also the population’s trait structure. This approach is flexible and able to incorporate the effect, both instantaneous and delayed, of environmental stressors on multiple traits, making it ideal for accounting for the complex life history of Ae. albopictus.

The model predicts the dynamics of a population of mosquitoes arising from a single water body of fixed dimensions. The model inputs are environmental variables from the location being simulated and the outputs are predictions of population and trait dynamics (Fig. 5). Adult mosquitoes oviposit eggs either onto the surface of the water or around the sides of the habitat, with the proportion of eggs being placed around the side of the habitat increasing as water level decreases73. The eggs placed into or around the water body express either a diapausing or non-diapausing phenotype which is determined by a maternal effect in response to falling temperatures and decreasing photoperiod24. Egg diapause is a form of regular dormancy that allows eggs to withstand cold temperatures and is the primary mechanism by which Ae. albopictus populations persist and overwinter in temperate climates74. Once development is complete, both diapausing and non-diapausing eggs either become quiescent, a form of irregular dormancy that allows persistence through dry periods or immediately hatch into larvae20. Quiescence continues until the dormant egg is inundated by precipitation after which it immediately hatches.

In the model we assume that the water body only varies in response to changes in temperature, the accumulation of precipitation, and through evaporation of standing water and is otherwise identical in every respect between locations. We assume that all eggs initially either express the active or diapausing phenotype. All eggs expressing the active phenotype develop at a temperature-dependent rate after which they either hatch or become quiescent according to the water level in the habitat. After completing development, diapause eggs are assumed to remain dormant until a critical temperature and photoperiod is reached after which they immediately hatch or become quiescent24. We place restrictions on both the production of diapause eggs and upon their release from diapause to ensure that they are fully developed when they hatch. The survival of eggs through the egg stage is assumed to be determined solely by temperature with this relationship differing between egg phenotypes75. We assume all eggs in the quiescent egg class are identical regardless of previous phenotype and assume that quiescent eggs hatch when the water level rises.

Larval mosquitoes compete for available resources in the aquatic habitat and are assumed to develop at a rate determined by both the current temperature and the intensity of larval competition9. The ecology of larval mosquitoes is complex, and to produce a coherent model with the data currently available it is necessary to make broad simplifying assumptions about the life history of larvae. We assume that all larvae are functionally identical, regardless of larval instar or previous experience of dormancy during the egg stage. Therefore, we only consider a single larval class76. It is assumed that there is no interspecific competition either for resources or through predation and that the only intraspecific competition is for resources and space77,78. Available resource is assumed to be consumed in its entirety and to be replenished daily, representing the product of temperature-dependent metabolic processes in the larval environment. Once larval development is complete pupation begins with the development of pupae assumed to depend solely on temperature75. The container habitats that Ae. albopictus prefer are generally small and vulnerable to flushing, a process whereby the body of water overflows and individuals are swept away21. We represent flushing by increasing the mortality of larvae and pupae whenever the height of the water in the habitat exceeds the height of the container and rainfall is sufficiently intense. These small containers are also susceptible to drying out, and whenever all water within the container evaporates, all non-quiescent juveniles are assumed to die out.

Adult mosquitoes experience developmental plasticity in response to their experience of temperature and intraspecific competition as larvae79. This is represented by tracking each individual’s experience of the biotic and abiotic environment during development and using this to determine the trait that individuals express as adults. Specifically, each individual’s experience of the average temperature and the average food available per larvae per day over the course of the larval period is used to predict that individual’s wing length. Wing length is then used to determine the fecundity and longevity of that individual as an adult80. Both fecundity and longevity are then further modified by the current temperature meaning that within the model adult traits respond to both current and historic environmental conditions43,75. The adult resources necessary to survive and complete the reproductive cycle, such as sugar and blood meals, are always assumed to be in excess and the model does not represent any effects of competition between adults81. The production of cold resistant diapausing eggs is triggered when falling temperatures and photoperiod reach a critical threshold24.

Model details

The abiotic environment. Abiotic environmental variation is incorporated into the model through temperature, photoperiod, evaporation, and precipitation from the location in which the modelled population is situated. Meteorological data is sourced from the ERA5-land climate reanalysis dataset82. Climatic variables in this dataset are available at a 0.1° × 0.1° resolution at an hourly time interval and are processed to produce a daily mean temperature and the total accumulated precipitation and evaporation each day. For use in the model each environmental variable is extended to continuous time using splines, which were then tested against the climate data to ensure there were no errors in the interpolation. Splines were created using the function Spline1D from version 0.5.2 of the Dierckx package, a wrapper from the Fortran library of the same name83.

The temperature at time t is denoted T(t), and for simplicity, we assume that the temperature of the larval pool and the air temperature are always the same84. Temperature is used throughout the model to instantaneously alter the traits (survival and development rate) expressed by individuals in each developmental stage. Temperature plays a further role in two instances of delayed phenotypic plasticity represented in the model. The average temperature an individual experiences over the larval period in combination with the average amount of food available per larva per day determines the wing lengths of emerging adults. Also included is a maternal effect where the temperature and the photoperiod at time t, denoted ψ(t), determine the proportion of produced eggs that express a diapausing phenotype. A combination of temperature and photoperiod subsequently determines when diapausing eggs are released from dormancy.

The total precipitation and evaporation each day, denoted ρ(t) and ϵ(t) respectively, are used to simulate the dynamics of the aquatic developmental habitat. The model simulates the dynamics of a single cylindrical water body with a surface area of μ mm2 and maximum volume of V mm3. Throughout the subsequent analysis these parameters are kept constant across the range of locations considered and have been selected to be consistent with productive habitats observed in the field85. Water accumulates in the habitat through precipitation and leaves through both evaporation and overspill. Overspill occurs when the amount of water in the pool exceeds the volume of the pool, at which point any additional precipitation is not added to the pool's total volume. If an overspill event happens concurrently with an intense period of precipitation, a spike in mortality of the larval and pupal classes occurs, representing individuals being flushed away21. These small containers are also susceptible to drying out, and whenever all water within the container evaporates, all non-quiescent juvenile stages (active eggs, diapausing eggs, larvae, and pupae) die-out. Habitat water volume, W(t), also governs the proportion of eggs produced that are quiescent and when they are released from dormancy.

Eggs. The number of active, quiescent and diapausing eggs present in the population at time t is denoted by Eγ(t), EQ(t) and ED(t) respectively. Active eggs develop at a temperature-dependent rate, immediately hatching once development is complete if they are at or below the current water level and moving into the quiescent class otherwise75. Quiescent eggs remain dormant until inundated by precipitation after which they immediately hatch20. Hatched eggs enter the larval class. Diapausing eggs are cold resistant and are produced in response to a maternal effect. They remain dormant until a critical temperature and photoperiod threshold is reached after which they either immediately hatch or become quiescent24. In other Aedes species, it has been demonstrated that adult body size is related to egg size and that larger egg sizes increase desiccation resistance86,87. Although weak effects of egg size on survival have been observed in Ae. albopictus, there is insufficient empirical data to justify the inclusion of such a function in this model. Therefore, it is assumed that the only maternal effect is whether an egg is active or diapausing88.

For active eggs, our model includes two state-equations that describe how the stage duration and through-stage survival of active eggs vary through time, which we denote by τEγ(t) and SEγ(t), respectively. To describe how these quantities evolve we parametrise two temperature-dependent reaction norms gEγ(T) and S^Eγ(T) which describe the development rate of active eggs and the proportion of active eggs that survive through the stage when held at constant temperature T. Note that S^Eγ(T) is distinct from SEγ(t), the former is a function of temperature describing the reaction norm, the latter is a function of time, describing time-dependent survival. We use a similar notation for survival and reaction norms of the other life stages. From the underlying Nisbet and Gurney (1983) framework we inherit a relationship between the development rate, the probability of through-stage survival and the mortality rate. We therefore use the temperature-dependent reaction norms for development and through-stage survival to define an expression for the temperature-dependent mortality rate of active eggs, δEγ(T) = −log(S^Eγ(T))gEγ(T)72. The forms of the reaction norms used for egg survival and development were selected by fitting function forms that have been used in previous models to life-history data. We trial quadratic and Briere functional forms for the development rates, and quadratic and Gaussian functional forms for through-stage survival with each functional form being truncated at a value close to 08. Functional forms were fitted using a weighted non-linear least squares approach and AIC was used for model selection, using the nls function in R89. The parameter values are reported in Supplementary Note 6, and the data used to perform this parametrisation was taken from a variety of published constant temperature laboratory experiments75,90–97. The best fitting functional form for the development rate of active eggs was a quadratic given by

gEγ(T) = max(σ11T2 + σ12T + σ13, 0.01) (2)

The best model for the through-stage survival of active eggs was a Gaussian with a functional form

S^Eγ(T) = σ22σ21√2π exp(− 1 2 (T − σ23)2 σ22) (3)

The result of this fitting is provided in Fig. 6A and B respectively. Egg quiescence is a mechanism of irregular dormancy that allows a reservoir of eggs to persist through periods when the habitat has completely dried out20. After maturing out of the active stage, a proportion of eggs are assumed to enter the quiescent egg class dependent upon the water level within the container habitat. Quiescent eggs contain fully developed pharate larvae that remain dormant inside the egg until inundated by water, and therefore within the quiescent egg class, no development occurs. The proportion of eggs that enter the quiescent class at time t is represented by the function Q(t) which is defined by

Q(t) = (V−W(t))/V, if W(t) < W(t−1) 0, otherwise (4)

This function links the proportion of eggs that become quiescent to the height of the water in the container habitat. The condition W(t) < W(t−1) ensures eggs only become quiescent if the water level on day t is lower than that on day t−1. This choice of functional form ensures that when the habitat completely dries out, W(t) = 0, all eggs become quiescent upon maturation and so Q(t) = 1. As the water level rises a lower proportion of eggs enter the quiescent class until the habitat is completely filled, W(t) = V, and therefore Q(t) = 0. Within the quiescent egg class, it is assumed that the same temperature-dependent mortality rate used for active eggs is suitable for quiescent eggs.

As egg quiescence is an irregular form of dormancy, maturation out of the quiescent class is triggered by inundation rather than the completion of development. This is represented through the function hQ(t), which describes the rate at which eggs that hatch out of the quiescent class at time t and is defined by

hQ(t) = (1 × (V−W(t)))/V, if W(t) > W(t−1) 0, otherwise (5)

When the habitat is completely full, W(t) = V, and so hQ(t) = 1 and all quiescent eggs are inundated and hatch over the course of a day. The condition W(t) > W(t−1) ensures that no eggs are released from quiescence when the water level at time t is lower than the water level at time t−1 and therefore no new eggs would become inundated.

Egg diapause is a mechanism of regular dormancy that allows the persistence of populations through cold winter months and is governed by a maternal effect24. This is represented in this model by the diapause function, D(t), which describes the proportion of adults at time t that are producing non-diapause eggs and is defined by

D(t) = (1 − 1/(1+15eψ(t,l)−ϕ(l))), if T(t) ≤ 18, and ψ(t) ≤ ψ(t−1) 1, otherwise (6)

with the functional form taken from Lacour et al.24 and where ϕ(l) is the critical photoperiod threshold required to induce the production of diapausing eggs in adults at latitude l defined by

ϕ(l) = 0.1|l| + 9.5 (7)

according to the relationship determined by Armbruster et al.74. The emergence of eggs from diapause is triggered by the current environmental conditions exceeding a critical temperature and photoperiod threshold98. We assume that once the environmental thresholds are reached, all eggs are released from diapause hatch over the course of a day. This is represented by the function hD(t) which is defined by

hD(t) = 1, if T(t) ≥ 12.5, and ψ(t) > ϕ(l), and ψ(t) > ψ(t−1) 0, otherwise (8)

The mortality rate of diapausing eggs is defined by

δD(t) = 0.01, if T(t) > −12 0.1, otherwise (9)

This choice of mortality rate allows diapausing eggs to survive exposure to extreme temperatures for short periods45. This ensures that brief periods of intense cold do not kill off the whole population whilst ensuring that sustained low-temperature conditions remain unsuitable.

Larvae. Larval Ae. albopictus consume a range of bacteria, plants and detritus. Both the time that larvae take to develop and the proportion of individuals that survive through the larval stage is assumed to depend on the intensity of intraspecific competition for resource in addition to temperature9. Predicting how much resource is available for larval Ae. albopictus in a container habitat from climate data alone would be a significant undertaking for even a limited geographical range. Resource dynamics are complex, and the amount of resources available in a container habitat varies with location, micro-climate, and community composition amongst many other factors that are not accounted for here99. Instead, we assume that the amount of food available in the container is temperature-dependent, proportional to the volume of water in the habitat and completely independent of location. The gross primary productivity of the larval habitat is then defined to be

F(t) = 10−6log10(0.45 + 0.095T(t))W(t) + fd (10)

where the logarithmic term represents the products of respiration following the work of White et al. and the constant term, fd, is a reservoir of nutrition from the decay of detritus100. It is further assumed that larvae divide the available food equally and completely (i.e., scramble competition)101. The amount of food available per larvae per day, α(t), can therefore be expressed by

α(t) = F(t)/L(t) (11)

As with active eggs, we include two state equations that describe how the stage-duration and through-stage survival of larvae changes through time, which we denote by τL(t) and SL(t), respectively. To describe how these quantities vary we parametrise two temperature and food-dependent reaction norms gL(T, α) and S^L(T, α) which describe the development rate of larvae and the proportion of larvae that survive through the stage when held at constant temperature T and resource density α. From the underlying Nisbet and Gurney (1983) framework we inherit a relationship between the development rate, the probability of through-stage survival and the mortality rate. We therefore use the temperature and food-dependent reaction norms for development and through-stage survival to define an expression for the temperature and food-dependent mortality rate of larvae, δL(T, α) = −log(S^L(T, α))gL(T, α)72. The larval mortality term, δL(T(t), α(t)), accounts for the dynamics of the container habitat and density dependence. The water level of the container habitat varies with evaporation and precipitation. Spikes in larval mortality can therefore either be induced through the habitat drying out or through the flushing of larvae when the habitat overflows. This is represented in the larval mortality rate by setting δL(t) = δd whenever the habitat dries out, (i.e., W(t) = 0) and similarly setting δL(t) = δf whenever the habitat overflows and the intensity of rainfall is sufficient to initiate flushing (i.e., W(t) = V and μρ(t) > 0.5V). The final component of the larval mortality function is an overcrowding term that increases the mortality of larvae when the larval density is high, which is defined by

δL(t) = min(−log(S^L(T(t), α(t)))gL(T(t), α(t)) + exp(exp(W(t)−1)), 0.3) (12)

Pupae. The pupal development rate and survival are assumed to depend solely on temperature as there is a lack of experimental data to quantify the effect of developmental plasticity from the larval stage. We include two state equations that describe how the stage duration and through-stage survival of pupae vary, denoted τP(t) and SP(t) respectively. To describe how these quantities vary we once again parametrise reaction norms for the development rate and survival of pupae held at constant temperature T, denoted by gP(T) and S^P(T) respectively with the forms of the reaction norms shown in Fig. 8. These reaction norms are fitted using the same procedure described for the egg stage, with a Briere function best describing the relationship between pupal development rate and temperature with form

gP(T) = max(σ31T(T−σ32)(σ33−T)1/σ34, 0.01) (13)

A quadratic functional form is chosen for the relationship between through-pupal stage survival and temperature independently of hydrological processes

S^P(T) = max(σ41T2 + σ42T + σ43, 0.01) (14)

parametrised using data from various sources75,93,95–97,106,109,121,122. The survival and mortality parameters are modified as in the larval stage to induce mortality when W(t) = 0 and to account for flushing such that

δP(t) = −log(SP(t))gP(t), if 0 < W(t) ≤ V and μρ(t) ≤ 0.5V δd, if W(t) = 0 δf, if W(t) = V and μρ(t) > 0.5V (15)

Adults. Ae. albopictus experiences developmental plasticity in response to the temperature, food availability, and conspecific density experienced during juvenile development that contributes towards temperature-dependent mortality and fecundity of adults9,80. This is included in the model by using an individual’s experience of larval competition and temperature through development to predict the wing length that the individual will express as an adult. This wing length is then used along with the current temperature to define the mortality rate of adults and the rate at which eggs are produced. This overlooks a known crowding effect whereby high densities of larvae inhibit development independently of competition for resource79. This choice is necessary as much of the experimental data used to parametrise reaction norms does not control for larval density and food availability simultaneously and therefore this relationship cannot be included in the model. Also omitted are the effects of the availability of adult resources on the production of eggs, and other varied processes demonstrably important in the life history of adults118,123. Using the same GAMM fitting procedure as described for larvae, we fit a 2-dimensional reaction norm for the wing length of emerging adults, w(Tavg(t), α(t)), as predicted by the average larval temperature and food available per larvae per day over the course of the larval period (Fig. 9A). We also fit the reaction norm for the time until 50% adult mortality, τA50(w, T(t)), as predicted by wing length and temperature (Fig. 9B). The parametrisations are carried out using data from published laboratory data from experiments that used constant temperatures5,9,79,90–94,96,108–111,113,115–118,121–143. In the model, the effects of developmental plasticity are represented by assigning each individual to an environmental class that reflects their experience of the environment through development. This is achieved by partitioning the adults into m environmental classes that define the wing length, wj, the adults are predicted to express based upon their experience of temperature and resource availability as larvae. We divide the possible wing lengths an adult could express, (wmin, wmax) into m sub-intervals. The transition function g(w) = wj if w is in the jth subinterval of (wmin, wmax).

To determine the average amount of food available for larvae maturing into adulthood at time t, we define

α(t) = 1/τL(t) ∫t t−τL(t) F(s)/L(s) ds (16)

Then, we define the transition function, w(Tavg(t), α(t)), which determines which adult class a fully developed juvenile is assigned to, with

w(Tavg(t), α(t)) = 1, if g(w(Tavg(t), α(t))) = wj 0, otherwise (17)

for j ∈ 1, ..., m.

The fecundity of an individual expressing wing length w is then given by

q(w) = 0.5exp(σ51 + σ52w) (18)

noting that a factor of 0.5 is required as only females produce eggs80. The duration of a gonotrophic cycle is assumed to be temperature-dependent and to be the reciprocal of the biting rate. This is parameterised using a Briere function taken from Mordecai et al.8 of the form

G(T(t)) = max(σ61T(T−σ62)(σ63−T)1/σ64, 0.01) (19)

from which the rate at which an individual with wing length w produces eggs at temperature T(t) can be determined. The number of eggs an adult in environmental class j produce eggs over the course of a gonotrophic cycle is defined by

qj = q(wj) (20)

Similarly, for individuals in environmental class j that express wing length wj the mortality rate of adults is defined to be

δAj(T(t)) = −log(0.5)/τA50(wj, T(t)) (21)

Population model. Following the model framework defined in Brass et al.12, based upon the work of Nisbet and Gurney72, the stage-structured model to predict the population dynamics of Ae. albopictus is

dEγ(t)/dt = REγ(t) − MEγ(t) − δEγ(t)Eγ(t) (22)
dED(t)/dt = RED(t) − MED(t) − δED(t)ED(t) (23)
dEQ(t)/dt = REQ(t) − MEQ(t) − δEQ(t)EQ(t) (24)
dL(t)/dt = RL(t) − ML(t) − δL(t)L(t) (25)
dAj(t)/dt = RAj(t) − δAj(t)Aj(t), for j ∈ 1, ..., m (26)

The variable developmental delay terms are defined such that

dτEγ(t)/dt = 1 − gEγ(t)/gEγ(t−τEγ(t)) (27)
dτL(t)/dt = 1 − gL(t)/gL(t−τL(t)) (28)
dτP(t)/dt = 1 − gP(t)/gP(t−τP(t)) (29)

The through-stage survival terms are defined such that

dSEγ(t)/dt = SEγ(t) (gEγ(t)δEγ(t−τEγ(t))/gEγ(t−τEγ(t)) − δEγ(t)) (30)
dSL(t)/dt = SL(t) (gL(t)δL(t−τL(t))/gL(t−τL(t)) − δL(t)) (31)
dSP(t)/dt = SP(t) (gP(t)δP(t−τP(t))/gP(t−τP(t)) − δP(t)) (32)

Recruitment and maturation terms defined by

REγ(t) = D(t) (∑mj=1 qjAj(t))/G(t) (33)
MEγ(t) = REγ(t−τEγ(t))SEγ(t−τEγ(t)) (34)
RED(t) = (1−D(t)) (∑mj=1 qjAj(t))/G(t) (35)
MED(t) = hD(t)ED(t) (36)
REQ(t) = Q(t)(MEγ(t) + MED(t)) (37)
MEQ(t) = hQ(t)EQ(t) (38)
RL(t) = (1−Q(t))(MEγ(t) + MED(t)) + MEQ(t) (39)
ML(t) = RL(t−τL(t))SL(t−τL(t)) (40)
RAj(t) = w(Tavg(t), α(t)) gP(t) ML(t−τP(t))SP(t−τP(t)) (41)

where C(t) is an impulse function that introduces active eggs to the system at time 0. The initial conditions are Eγ(0) = 100, ED(0) = EQ(0) = L(0) = Aj(0) = 0 for j ∈ 1, ..., m, SX(0) = exp(−τX(T(0))δX(T(0))) for X ∈ E, L, P. The initial history is for all t ≤ 0 and all simulations are performed with m = 200 unless otherwise noted. The model was simulated in Julia version 1.8 using the package DifferentialEquations144,145.

SEIR Model

The model for the population and trait dynamics of Ae. albopictus is incorporated into a susceptible-exposed-infected-resistant (SEIR) model for dengue virus vectored by Ae. albopictus (Fig. 5). We specifically consider the introduction of a single serotype of dengue virus into a completely susceptible human population. This simplified representation of the urban dengue transmission cycle is most applicable to non-endemic regions, such as Europe, the United States, and eastern Asia, where the likelihood of multiple serotypes circulating simultaneously is reduced and therefore the human population is less likely to develop serious complications19. The human population is assumed to be susceptible to infection, HS, recovered, HR, and infected, HI, 146. It is assumed that the size of the human population is constant with the exception of the entry of infected humans into the population. Human population density is estimated using the Gridded Population of the World, Version (GPWv4) estimates for 2015147. Ae. albopictus is an opportunistic feeder, and, in the absence of humans, will bite other vertebrate species. These bites cannot propagate dengue virus (outside of the geographically limited sylvatic transmission cycle148), but do allow for egg production and so can maintain the mosquito population. To represent this a population of non-human food sources is included in the model, HB, which is estimated from field observations of rates of anthropophagy and a meta-analysis of mammal densities149,150. Mosquitoes are assumed to bite at a temperature-dependent rate of b(t) bites per mosquito per day that is inversely proportional to the length of the gonotrophic cycle such that b(t) = 1/G(T(t))43. The proportion of uninfected mosquitoes that become infected after biting an infected human is denoted hv(t), a function that takes the form used in Liu-Helmersson et al. (2014)151. After the extrinsic incubation period, τEIP(t), the proportion of mosquitoes that successfully maintains infection is denoted hv(t) 43,151. After the transmission of the virus, the mosquito remains infectious for the remainder of its life, with τEIP(t) implicitly accounting for the latent period between exposure and infection without the need for the explicit inclusion of an exposed class. In this model, large mosquitoes are assumed to have no direct advantage in the act of transmission over small mosquitoes and do not bite more aggressively nor transmit disease more competently. The only advantage large mosquitoes have over small is in their extended longevity, which allows a greater proportion of mosquitoes contracting an infection to survive to pass that infection on, especially under stressful environmental conditions.

As the population model predicts the number of mosquitoes arising from a single habitat and we now wish to predict the disease transmission dynamics for a human population, it is necessary to reconcile both populations onto the same spatial scale. Mark-recapture studies have demonstrated that adult female Ae. albopictus can disperse up to a kilometre from their release site152. Further, analysis of dengue epidemics has shown that disease dynamics can best be understood in 2km × 2km grids153. For these reasons the model is simulated at a spatial scale of 4km2, using the number of humans per 4km2 for the human population in combination with an estimate of the number of larval habitats per 4km2 for the mosquito population which is halved to represent that only female mosquitoes bite. The number of productive larval habitats per 4km2, κ, is assumed to be fixed between locations and within seasons and estimated from field, surveys of larval habitats in Emilia-Romagna, Italy, Hawai’i, USA, and Maryland, USA25,154–156. The structured human population is represented by

dHS(t)/dt = −R(t) (42)
dHI(t)/dt = R(t−τIIP) − R(t−τIIP−τREC) + C(t) − C(t−τREC) (43)
dHR(t)/dt = R(t−τIIP−τREC) + C(t−τREC) (44)

where

R(t) = vh(t)b(t)2κ ∑mj=1 Ij(t) HS(t)/HT (45)

and HT(t) = HS(t) + HI(t) + HR(t) + HB and C(t) is an impulse defined as in the mosquito population model that initiates the transmission cycle through the introduction of infectious humans at a predetermined time that varies from outbreak to outbreak. Note that the inclusion of non-human food sources introduces a dilution effect, as when HB is large in relation to HT(t) the value of R(t) will be low. This means that when the human population size is small, and therefore most blood meals are assumed to come from hosts that cannot become infected with dengue virus, the model predicts that only a low number of total bites by infected mosquitoes result in susceptible humans becoming infected. The omission of an explicit state-equation for humans who have been exposed to infection, but who are not yet infectious is a consequence of the intrinsic incubation period being implicitly accounted for in the delay structure through the inclusion of the term τIIP.

To represent infections in the mosquito population, infectious adult classes are added, one for each environmental class, denoted Ij(t). The number of mosquitoes in environmental class j that become infected at time t is denoted RIj(t), and defined by

RIj(t) = hv(t)b(t)Aj(t) HI(t)/HT (46)

and the number of mosquitoes in environmental class j that become infectious at time t is denoted MIj(t) and defined by

MIj(t) = RIj(t−τEIP(t))SEIPj(t) (47)

where SEIPj(t) is the proportion of individuals in environmental class j that survive the infectious period. The dynamics of infectious mosquitoes can therefore be described by

dIj(t)/dt = MIj(t) − δIj(t)Ij(t), for j ∈ 1, ..., m (48)

dτEIP(t)/dt = 1 − gEIP(t)/gEIP(t−τEIP(t)) (49)
dSEIPj(t)/dt = SEIPj(t) (gEIP(t)δEIP(t−τEIP(t))/gEIP(t−τEIP(t)) − δEIP(t)) (50)

where gEIP(t) is defined as in Brady et al.43. The equation for the rate of change of adults in environmental class j is modified to include a term representing mosquitoes becoming exposed to infected humans and thus leaving the uninfected adult classes such that

dAj(t)/dt = RAj(t) − MIj(t) − δAj(t)Aj(t), for j ∈ 1, ..., m (51)

Where MIj(t) describes the rate that adults become infective. The mosquito dynamics of the SEIR model are therefore given by

dEγ(t)/dt = REγ(t) − MEγ(t) − δEγ(t)Eγ(t) (52)
dED(t)/dt = RED(t) − MED(t) − δED(t)ED(t) (53)
dEQ(t)/dt = REQ(t) − MEQ(t) − δEQ(t)EQ(t) (54)
dL(t)/dt = RL(t) − ML(t) − δL(t)L(t) (55)
dAj(t)/dt = RAj(t) − MIj(t) − δAj(t)Aj(t), for j ∈ 1, ..., m (56)
dIj(t)/dt = MIj(t) − δIj(t)Ij(t), for j ∈ 1, ..., m (57)

The variable delay terms satisfy

dτEγ(t)/dt = 1 − gEγ(t)/gEγ(t−τEγ(t)) (58)
dτL(t)/dt = 1 − gL(t)/gL(t−τL(t)) (59)
dτP(t)/dt = 1 − gP(t)/gP(t−τP(t)) (60)
dτEIP(t)/dt = 1 − gEIP(t)/gEIP(t−τEIP(t)) (61)

The through-stage survival terms satisfy

dSEγ(t)/dt = SEγ(t) (gEγ(t)δEγ(t−τEγ(t))/gEγ(t−τEγ(t)) − δEγ(t)) (62)
dSL(t)/dt = SL(t) (gL(t)δL(t−τL(t))/gL(t−τL(t)) − δL(t)) (63)
dSP(t)/dt = SP(t) (gP(t)δP(t−τP(t))/gP(t−τP(t)) − δP(t)) (64)
dSEIPj(t)/dt = SEIPj(t) (gEIP(t)δEIP(t−τEIP(t))/gEIP(t−τEIP(t)) − δEIP(t)) (65)

The initial conditions are the same as those for the disease free mosquito model, with the additional history HI(0) = HR(0) = 0 for t ≤ 0 and HS(0) is the initial human population density as estimated from the Gridded Population of the World, Version (GPWv4) for 2015147. The mosquito population is assumed to be initially entirely uninfected and mosquitoes are introduced to the system as described in Equation (33). Refer to the table in Supplementary Note 6 for an overview of all parameters and notation.

Rt derivation

To quantify how disease risk changes between regions we derive an expression for the reproduction number, Rt, which describes the number of infections produced by the introduction of a single infected human into the population, accounting for the potential changes in the population structure as the outbreak progresses158. We assume that a single infected human is introduced precisely at the end of the intrinsic incubation period and so can be bitten for the full duration of the infectious period. We note that for a mosquito to become infectious at time t the infecting bite must have occurred at time t−τEIP(t). Therefore, the number of mosquitoes that become infectious at time t per day per 4km2 can be expressed as

∑mj=1 (gEIP(t−τEIP(t))SP(t−τEIP(t))) b(t−τEIP(t))hv(t−τEIP(t))2κAj(t−τEIP(t))SEIPj(t) (66)

The infected human is assumed to remain infectious for τREC days, and so the total number mosquito per 4km2 that become infected as a result of the introduction of an infected human at time t−τEIP(t) is the case reproduction number for human to mosquito transmission and can be expressed as

∑mj=1 ∫t+τREC t (gEIP(s−τEIP(s))SP(s−τEIP(s))) b(s−τEIP(s))hv(s−τEIP(s))2κAj(s−τEIP(s))SEIPj(s) ds (67)

The number of susceptible humans that a single infectious mosquito bites per day at time t can be expressed as

b(t)vh(t)HS(t)/HT (68)

and we can therefore approximate the mosquito to human case reproduction number, the total number of new human infections caused by a mosquito that became infectious at time t, by

∫t+1=t δAj(t)b(s)vh(s)HS(s)/HT ds (69)

Note that this is only an approximation as the lifespan of the infected mosquito does not account for the number of bites that occur in the time between t and t+1=δAj(t). The human-to-human case reproduction number is the square root of the product of the human-to-mosquito and mosquito-to-human case reproduction numbers157. Hence, the number of new infections caused by the introduction of a single infectious human at time t−τEIP(t) can be approximated by

Rt−τEIP(t) = (∑mj=1 ∫t+τREC t (gEIP(s−τEIP(s))SP(s−τEIP(s))) b(s−τEIP(s))hv(s−τEIP(s))2κAj(s−τEIP(s))SEIPj(s) ds) × ∫t+1=t δAj(t)b(u)vh(u)HS(u)/HT du (70)

From Equation (70), we can compute Rt0, where t0=t−τEIP(t), directly from the numerical output. We note that in Equation (70) the mosquito-to-human case reproduction number appears inside the integral for the human-to-mosquito case reproduction number. This accounts for the fact that depending on the human population and environmental conditions at the time a mosquito becomes infected, the resultant number of mosquito-to-human cases will vary. We test the validity of this approximation in two ways. In Supplementary Note 4 we demonstrate that for a population held at constant temperature when Rt = 1 − ϵ case numbers decrease and that when Rt = 1 + ϵ case numbers increase for a suitable ϵ. Further, Supplementary Fig. 70 shows the average total annual number of dengue cases over the same spatial region and temporal extent that is considered in Fig. 4.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors intuitively make the case for temporal effects of the environment on phenotypic plasticity. Their observations on the shortcomings of previous models are very astute – I would describe their model as being much more biologically-informed than existing models, and hence intuitively more likely to describe real-world nuances which are missed by more basic analyses. These hypotheses seem to pan out when validating their model against field data. However, I do have some criticisms of the data they have used which will require some additional systematic analyses and in some cases describing the unavoidable gaps in the data. I focus mainly on validation of the model and data utilization in my comments as I believe these to be the weakest parts of the work. Overall though, and assuming my concerns can be adequately addressed, I believe this to be an excellent and very substantial advancement in disease risk prediction and commend the authors.\n114-129 and figure 1: I am not sufficiently familiar with the individual filed observations and data produced. What the authors show is excellent consistence between their model and the field data from the listed papers. However, there lacks a convincing described methodology to show me that these studies are representative/exhaustive. I am certainly not arguing that the authors have cherry-picked only the studies whose data align with their predictions, but the authors do not do a sufficient job in convincing beyond reasonable doubt that this is not the case. I would like to see a paragraph detailing essentially a mini-systematic review where they have unbiasedly identified all papers with relevant data, and then systematically compared their predictions to all of them, rather than just the displayed sub-set. And this should be broken down by the type of data in fig 1 (i.e. oviposition activity, number of larvae, number of adults…). From what is written, it could be that this is an exhaustive list, a representative list, or just a cherry-picked list of neatly aligning data – convince the reader that this is one of the first two, not the latter.\n137-149 and figure 2: Same comment as for figure 1 really. There have been a lot of dengue outbreaks, why are these four specifically chosen and represented here? Presumably you can only include albopictus-driven outbreaks, so this may reduce the number considerably, but there must be more than these presented. Again, I’d like to see a systematic methodology of selection of the examples, and ideally a comparison of model to outbreak for all useable outbreaks (I’m aware this may be too many to be feasible, in which case there needs to be a very convincing ‘representative vs cherry-picked’ explanation). And the same for vector control activities, in Guangzhou 2013-2014, you very neatly show the vector control coincides with a ‘missing’ observed dengue peak that your model predicts, is this an exception or does this pattern follow for a large proportion of ‘missing’ observed peaks – for this I’m fully aware that data on activities may be lacking, so a discussion paragraph here would suffice.\nGeneral statement – you assess the model predictions of outbreaks by comparing to real-world albopictus-driven outbreaks, by necessity of course. However, might this lead to a bias as albo-only regions may be different in terms of climate/land use to albo + aegypti regions (e.g. you have two small island examples). Not much can be done about this, but I suggest you comment on the effect on the utility and validation of your model. How can we know it is generalizable to aegypti endemic regions and e.g. mid-continental tropical regions?\nFigure 4 – this figure gives me more confidence that the earlier studies were representative rather than being cherry-picked. But again the full analysis/validation is limited to Europe where there is albo but no aegypti. Again I understand the need for this, but Europe is not a representative example of general albo-driven DENV transmission – the model could be working very well here, but less well in e.g. more tropical regions (where there is more DENV). Are there enough data/resolution and land area to do these analyses for La Reunion for example, it seems strange that it is not included here? Same for Japan, I would presume there is sufficient data there – even though it’s largely sub-tropical.\nI would like to see some discussion of the lack of data availability and the effect on the model. E.g. Figure 6 showing reaction norms for eggs vs temperature is a bit data-lacking, a couple of additional points could drastically change the shape of your curves. This is far more of a problem for Figure 8, especially B… Of course you have to work with the data that are available, so not a criticism of your work, but these shortcomings could have large effects on the model, especially at the higher end of the temperature range (incidentally, where you have fewer field-data for validation as noted before... This should be specifically addressed in your discussion).\nMinor\nParagraph lines 59-70 (and a few later examples such as line 176) – strictly, it is ‘dengue virus’ which is vectored by Ae. albopictus, not ‘dengue’, which refers to the disease itself in humans.\n(Remarks on code availability):\nBeyond my skillset", "reviewer_2": "In this manuscript, the authors use a SIR transmission model with mosquito population dynamics to investigate the impact of phenotypic plasticity in vector traits on dengue transmission in the context of Aedes albopictus. The mosquito population dynamics sub-model is more complex than those typically encountered in the epidemic modelling literature and has been introduced in a previous work. An important feature of the present model is its ability to yield a realistic dynamics of mosquito wing-length, which the authors use as a proxy for mosquito size and lifespan. This trait varies over time due to climatic and ecological factors, resulting in subtle dynamical changes in epidemic risk. Overall, I believe that the manuscript is well written and that the implications from this model are potentially very important. I particularly appreciated several of the considerations offered in the text. However, I believe that the implications for epidemic modelling are not sufficiently articulated in the text. Please find my detailed comments below:\n- The authors focus on the case of dengue and Aedes albopictus. Not being an expert in mosquitos, I find it unclear whether similar reasonings apply to other host-pathogen systems as well. Is phenotypic plasticity equally variable/relevant to epidemic spread in other mosquito species? I would thus suggest the authors to comment on the relevance of their findings to other systems. Eventually, the title could be made more specific as in its current form, the manuscript focuses heavily on a single application and it remains unclear how these result generalise.\n- The model appears to be accurate at reproducing mean wing length data whenever such data is available (Fig 1). However, Fig 1 and Supp. plots show that such data is often available within short time intervals. Moreover, wing-length does not vary a lot during these periods, whereas the model seems to suggest that this trait can be much more variable over a longer time window. Given that dynamical variation in wing-length underpins much of downstream results, I wonder if the authors could provide further evidence to support the “U” shape displayed by wing-length? At present, I fear that the model is extrapolating and I am not sure whether these extrapolations are justified or not.\n- The authors show, rather convincingly, that different modelling assumptions yield quantitative differences in epidemic dynamics. In particular, the authors show that dengue incidence and attack rate can vary considerably depending on how heterogeneity in wing-length is accounted for (or not accounted at all). However, while I do not have any doubts that the presented model is more biologically plausible than other alternatives, I’m not entirely convinced about its advantages in epidemiological modelling (e.g. in real-time tracking, planning for vector control). The practical superiority of the full model is also difficult to assess without proper model comparison: it could well be that simpler models fit the data equally well or better, but we don’t know. Here I’m not suggesting to perform a proper model selection exercise, but it would be good if the authors expanded on the implications of their eco-epidemiological model over simpler models. For example, the consideration that epidemic risk varies over time due to wing-length (and that mosquito abundance alone is not sufficient to predict transmission) is extremely interesting in my opinion. Potential avenues of discussion include implications for vector control and disease persistence during inter epidemic periods. Also, what are the consequences for R0 estimation? If I had to fit a simpler model to data (case data, serological data), would I get a wildly different estimate of R0 with respect to using the full model? How would these estimates then affect what we know about disease management?\n- How do estimates of suitability compare with another widely used suitability index, “Index P” (Obolski et Al, Methods in Ecology and Evolution, 2019; Nakase et Al, Scientific data, 2023)? Does the present index provide further insights on why Italy experienced an unprecedented dengue outbreak in 2023 (please refer to Branda et Al., Dengue virus transmission in Italy: surveillance and epidemiological trends up to 2023, medRxiv)?\nMINOR COMMENTS\n- L114 onwards: please briefly list which traits are being predicted.\n- Fig 2: make x axis labels uniform (2018-2019-…), increase line width…\n- L139 why complex?\n- Please increase x,y tick label size where this is too small to read.\n- Methods L23: “both phenotypes”: what does this refer to?\n- Methods L29: “of” should be “or”.\n- Methods Eq.(3): perhaps a typo in “1(…)”?\n- Methods L134: please add further details on how expressions like that for delta(T) are derived.\n- Methods L329: I would recommend the authors to label their model as a SEIR instead of a SIR, since the it does indeed account for the intrinsic incubation period. It also feels a more appropriate labelling choice for a disease such as dengue.\n- L351: “an infectious humans”\n- Methods L404: do you mean $R_{t’}?$\n- Methods L406: “This Is accounts”\n- Methods why relative humidity is not included??\n- Perhaps show population dynamics of mosquitoes and compare it with other population dynamics models?\n- Supplementary and Main: “Emilia-Romagna” and “Catania” are the correct names.\n- Supplementary P8 “Through through”\n- Supplementary P9 “decease”\n- Supp Fig 7A: It is not clear to me if the y axis denotes total numbers of proportions of adults.\n- Supp Fig 10 legend “Simluated”\n- S. 1.1.2 “year 2017”. Perhaps it should be 2012.\n- Supp Fig 14: the y label is cropped out.\n- S2.4: is it reasonable to assume that the population in Reunion is fully susceptible against DENV, given previous circulation of the virus in that area?\n- S3.2 “iin”\n- In the maps, both in the main and Supp material, it would be good to highlight the locations considered in this study."}, "author_response": "We thank the reviewers for their careful considerafion of our manuscript and perfinent suggesfions.\nIn response to these we have made changes that have materially improved the clarity of the\nmanuscript. In the response that follows the reviewer’s comments are in red and our point-by-point\nresponses are in black. Where substanfial changes to the manuscript have been made we provide\nthe altered text here, a line number where the change can be found, and addifionally highlight all\nchanged text in red in the manuscript.\nReviewer #1 (Remarks to the Author):\nThe authors intuifively make the case for temporal effects of the environment on phenotypic\nplasficity. Their observafions on the shortcomings of previous models are very astute – I would\ndescribe their model as being much more biologically-informed than exisfing models, and hence\nintuifively more likely to describe real-world nuances which are missed by more basic analyses.\nThese hypotheses seem to pan out when validafing their model against field data. However, I do\nhave some crificisms of the data they have used which will require some addifional systemafic\nanalyses and in some cases describing the unavoidable gaps in the data. I focus mainly on validafion\nof the model and data ufilizafion in my comments as I believe these to be the weakest parts of the\nwork. Overall though, and assuming my concerns can be adequately addressed, I believe this to be\nan excellent and very substanfial advancement in disease risk predicfion and commend the authors.\nWe thank the reviewer for their kind comments and useful crificism which we have now addressed\nthrough addifional clarificafions and substanfial new discussion paragraphs.\n114-129 and figure 1: I am not sufficiently familiar with the individual filed observafions and data\nproduced. What the authors show is excellent consistence between their model and the field data\nfrom the listed papers. However, there lacks a convincing described methodology to show me that\nthese studies are representafive/exhausfive. I am certainly not arguing that the authors have cherry-\npicked only the studies whose data align with their predicfions, but the authors do not do a sufficient\njob in convincing beyond reasonable doubt that this is not the case. I would like to see a paragraph\ndetailing essenfially a mini-systemafic review where they have unbiasedly idenfified all papers with\nrelevant data, and then systemafically compared their predicfions to all of them, rather than just the\ndisplayed sub-set. And this should be broken down by the type of data in fig 1 (i.e. oviposifion\nacfivity, number of larvae, number of adults…). From what is wriften, it could be that this is an\nexhausfive list, a representafive list, or just a cherry-picked list of neatly aligning data – convince the\nreader that this is one of the first two, not the lafter.\nThe list of studies included in our validafion were the result of an extensive literature search and to\nthe best of our knowledge we believe it is representafive of the literature. Our only\ninclusion/exclusion criteria were that samples needed to be at a least a monthly temporal resolufion\nin a locafion without Aedes aegypfi. We made addifional efforts to search the Japanese language\nliterature where many inifial studies on Aedes albopictus were conducted. Despite the importance of\nmosquito-borne disease, there is yet to be a well populated database for vectors that reports the\ndata needed to validate populafion dynamics at the temporal resolufion required. Much of the\nsurveillance data that does exist is difficult to access, as many countries collect mosquito surveillance\ndata at the state level and then do not publicly release it as it is a mafter of health policy. The\npopulafion dynamics that are publicly available in the academic literature are often not digifised,\nserving as secondary parts of larger studies looking into other quesfions and so are only parfially\nreported.\nThe studies we did find were found by inifially searching google scholar and web of science for terms\nsuch as “Aedes albopictus” and then snowballing from these to find addifional references. We\ngenerally encountered a high number of false posifives for more specific search strings, such as\n“Aedes albopictus seasonal abundance” and did not ulfimately use a systemafic approach to gather\nvalidafion data. Despite these limitafions we think the work undertaken here is sfill valuable as there\nis currently no centralised global repository for vector populafion dynamics data (although there are\nongoing efforts to change this) and so this represents the most comprehensive collecfion of\npopulafion dynamics for Aedes albopictus to date.\nThe inclusion criteria for our validafion dataset have now been made explicit in the text (L116-135).\n“We extensively validate the model predicfions of mosquito populafion, wing length, and disease\ndynamics by comparing them to field surveys of mosquito populafions and disease outbreaks across\nthe species range. To validate the populafion dynamics we use published surveillance data, including\nlife-stage specific populafion density esfimates and average trait data, from 41 locafions, across 14\ncountries and 4 confinents (see Supplementary Informafion 1 for the full set of validafions). To\nvalidate the model predicfions of dengue dynamics we compare them to historical dengue\noutbreaks, using reports from the outbreak to select a likely introducfion scenario for dengue virus\ninto the region and to define an area over which to simulate the model. Validafion data was obtained\nfrom a comprehensive search of published studies that either observed the populafion dynamics of\nAe. albopictus in the field or reported human cases of dengue. The literature search was conducted\nusing a snowball procedure - a systemafic search revealed a high number of false posifives. Studies\nwere included in our validafion if the data was collected with at least a monthly resolufion from a\nregion where Aedes aegypfi, a closely related vector species was absent. This exclusion criteria is\nnecessary as Ae. aegypfi is known to compete with Ae. albopictus for larval resources and is also a\nvector of dengue virus, factors which we do not account for here. Since substanfial variafion in case\nreporfing effort through space and fime is likely, our predicfions of dengue transmission are not\nrescaled and the validafions are intended to demonstrate that the model produces plausible disease\ndynamics (see also Supplementary Informafion 2). The comparisons presented in this paper\nrepresent the full set of datasets meefing the above criteria that were found in our literature search.”\nFigure 1 has been addifionally adjusted to show all locafions from which we were able to access\npopulafion dynamics data.\n137-149 and figure 2: Same comment as for figure 1 really. There have been a lot of dengue\noutbreaks, why are these four specifically chosen and represented here? Presumably you can only\ninclude albopictus-driven outbreaks, so this may reduce the number considerably, but there must be\nmore than these presented. Again, I’d like to see a systemafic methodology of selecfion of the\nexamples, and ideally a comparison of model to outbreak for all useable outbreaks (I’m aware this\nmay be too many to be feasible, in which case there needs to be a very convincing ‘representafive vs\ncherry-picked’ explanafion).\nAs the reviewer notes, we can only include outbreaks in this validafion where Aedes albopictus was\nthe sole vector implicated. Further, as our goal was to predict the temporal dynamics of human\ndengue cases, we also require that cases be reported at a sufficient temporal resolufion from a well-\ndefined region in space that does not consolidate data from several locafions. The subset of\noutbreaks shown in Figure 2 is the totality of the data we were able to find that meets these criteria\nwhich is now stated in the text as part of our response to the previous comment (L116-135).\nAnd the same for vector control acfivifies, in Guangzhou 2013-2014, you very neatly show the vector\ncontrol coincides with a ‘missing’ observed dengue peak that your model predicts, is this an\nexcepfion or does this paftern follow for a large proporfion of ‘missing’ observed peaks – for this I’m\nfully aware that data on acfivifies may be lacking, so a discussion paragraph here would suffice.\nData on vector control acfivifies is generally lacking. Guangzhou 2013-2014 is unusually well-studied\n(because of its remarkable magnitude) and so there is plenty of published informafion regarding the\nchanges in vector control strategies in the region. Although each of the other outbreaks was also\nfollowed by an increase in vector control acfivifies there is not enough informafion to aftribute the\n“missing” peaks to this increase (although it is our suspicion that this is the case). As suggested this is\nnow summarised in the new discussion paragraph addressing vector control (L312-331).\n“When designing vector control campaigns to implement intervenfions such as the sterile insect\ntechnique, gene-drive or transgenic control, it is common to rely on mathemafical models of vector\npopulafion dynamics \\citep{Ogunlade2023}. These models are used before implemenfing control\nmeasures to inform the required level of populafion suppression and the intensity of control\nmeasures necessary to achieve this target \\citep{Oliva2021, Bier2022, Powell2022}. Similarly, after\ncontrol measures have been implemented, mathemafical models are a pracfical method to assess\nhow effecfive an intervenfion has been by predicfing how an outbreak would have progressed in the\nabsence of control \\citep{Bonds2012}. Without this addifional analysis field-trials can only directly\nlink the efficacy of their intervenfion to a reducfion in mosquito abundance and not necessarily the\nreducfion in human cases of disease that they aim to achieve, two factors we show here are not\nnecessarily directly linked \\citep{Balatsos2024,Zhang2024}. The 2014 outbreak of dengue in\nGuangzhou provides a well-studied example showing the potenfial of our approach as we and\nprevious models are able to use mathemafical models to assess the efficacy of the implemented\ncontrol measures (see Figure~2) \\citep{Tang2016,Lin2016}. However, for the majority of the\noutbreaks we considered, this sort of retrospecfive analysis is not possible due to a lack of detail in\nthe reporfing of the nature, fiming, and intensity of control measures implemented. Further, to make\nreliable predicfions to assist in the design of vector control it is imperafive to use models that have\nbeen independently validated against real-world data from the target system \\citep{Brady2016}.\nExtensive and independent validafion of mathemafical models of systems of VBDs is currently not\nstandard, despite the ready availability of surveillance data and is a crifical step in ensuring these\npredicfions are robust.”\nGeneral statement – you assess the model predicfions of outbreaks by comparing to real-world\nalbopictus-driven outbreaks, by necessity of course. However, might this lead to a bias as albo-only\nregions may be different in terms of climate/land use to albo + aegypfi regions (e.g. you have two\nsmall island examples). Not much can be done about this, but I suggest you comment on the effect\non the ufility and validafion of your model. How can we know it is generalizable to aegypfi endemic\nregions and e.g. mid-confinental tropical regions?\nThe restricfion of the model to regions where Aedes albopictus is the only vector is a substanfial\nlimitafion inherent in our validafion approach. It is unfortunately necessary as we cannot account for\nthe effect that compefifion between albopictus and aegypfi will have on both populafion and disease\ndynamics. To address this our group is well underway in the development of a complementary model\nfor Aedes aegypfi that uses this same approach and is producing promising inifial results. Combining\nthese two models to understand dengue transmission dynamics in regions where both vectors\ncoexist will be a substanfial challenge but one that we hope will address the reviewer’s quesfions in\nthe long term. However, in the absence of the model for Aedes aegypfi there are sfill conclusions we\ncan make about the model’s generalisability which we now examine in a new discussion paragraph\naddressing model generalisability (L332-349).\n“Dengue transmission frequently occurs in regions where Ae. albopictus co-occurs with the primary\nvector of dengue, Ae. aegpyfi, and with which it competes for resources in larval habitats\n\\citep{Juliano1998}. This complicates our ability to predict the effect of trait dynamics on\ntransmission efficacy as the relafionship between adult wing length, temperature, larval density, and\nfood availability is asymmetrically altered by the degree of interspecific compefifion that individuals\nof either species experience \\citep{Juliano2004,Murrell2008,Lizuain2022}. Further, in endemic\nregions co-circulafing serotypes and previous exposure to infecfion complicate transmission\ndynamics by introducing an immune structure to the human populafion \\citep{Gibbons2007}. An\nextension of our model to these areas is possible but would addifionally require the development of\na comparable model for Ae. aegypfi, to quanfify the effect that interspecific compefifion has on the\nlife-history of both species and then to extend the disease dynamics to consider mulfiple co-\ncirculafing serotypes of dengue virus. However, even in these more complex transmission scenarios\nthe underlying interacfion between environmental drivers and vector trait expression is driven by the\nsame mechanisms of phenotypic plasficity that we explore here, as they are a fundamental feature\nof vector ecology. By developing our understanding of the drivers of disease dynamics that arise in\nregions with a single vector transmifting a single serotype of dengue virus, we are now well\nposifioned to begin addressing the more complex dynamics underlying pafterns of dengue\ntransmission in other regions.”\nFigure 4 – this figure gives me more confidence that the earlier studies were representafive rather\nthan being cherry-picked. But again the full analysis/validafion is limited to Europe where there is\nalbo but no aegypfi. Again I understand the need for this, but Europe is not a representafive example\nof general albo-driven DENV transmission – the model could be working very well here, but less well\nin e.g. more tropical regions (where there is more DENV). Are there enough data/resolufion and land\narea to do these analyses for La Reunion for example, it seems strange that it is not included here?\nSame for Japan, I would presume there is sufficient data there – even though it’s largely sub-tropical.\nIt is possible to perform the analysis we do in Figure 4 because in Europe Aedes albopictus is present\nand Aedes aegypfi is not, a limitafion that we now discuss in the new discussion paragraph on\nlimitafions. Although La Reunion has both highly spafially and temporally resolved case data the\nERA5 environmental data we use to drive the mosquito populafion dynamics is not high enough\nresolufion to perform a similar analysis, as the whole island consfitutes around 20 pixels. Although\nthis could potenfially be remedied through the inclusion of microclimafic models this is beyond the\nscope of what we aimed to achieve here.\nWe also presumed that we would be able to validate the model against Japanese dengue outbreaks.\nHowever, the only instance of autochthonous transmission we could find evidence for is the 2014\nTokyo outbreak that we have modelled (see Furuya (2015)). Similar efforts in the USA found very\nlimited evidence autochthonous transmission of dengue with the only instances outside of the\nHawai’i outbreak we model being limited to Puerto Rico, and small outbreaks in Florida and Texas\nwhere Aedes aegypfi was the vector implicated.\nI would like to see some discussion of the lack of data availability and the effect on the model. E.g.\nFigure 6 showing reacfion norms for eggs vs temperature is a bit data-lacking, a couple of addifional\npoints could drasfically change the shape of your curves. This is far more of a problem for Figure 8,\nespecially B… Of course you have to work with the data that are available, so not a crificism of your\nwork, but these shortcomings could have large effects on the model, especially at the higher end of\nthe temperature range (incidentally, where you have fewer field-data for validafion as noted before...\nThis should be specifically addressed in your discussion).\nThe availability of laboratory data is now addressed in a new data gaps discussion paragraph (L289-311):\n“The ability of mechanisfic modelling approaches to incorporate detailed vector ecology allows them\nto produce generalisable predicfions of disease risk. However, incorporafing these details requires\nsubstanfial amounts of life-history data that does not yet exist for many vectors, favouring the use of\nless detailed mechanisfic models than we develop here\n\\citep{Mordecai2019,Brady2014,Obolski2019}. Even for well-studied species such as Ae. albopictus,\nthere are data gaps in our current knowledge. Some of these gaps are life-stage specific, for example,\nwhen rearing adults for infecfion trials it is common to record the combined life-history traits of both\nlarval and pupal life-stages using a single generic class for juveniles, data which less detailed\nmechanisfic approaches can readily ufilise but that which we cannot. Although there are a sufficient\nnumber of studies exploring the effects of compefifion on larval life-history to produce well\nparametrised larval reacfion norms, there is very liftle available experimental data describing the life-\nhistory of pupae which are rarely specifically reared for study. As each life-stage may be affected\ndifferently to the same environmental stressor, life-stage specific gaps in our knowledge of vector\nlife-history potenfially mean that environmental limitafions on these life-stages are overlooked in our\npredicfions of vector dynamics and consequently disease transmission dynamics. Similarly, although\nthere has been a substanfial effort to determine the lower thermal limits of the species' reacfion\nnorms due to the species northward expansion, there is a more limited amount life-history data\navailable to describe the species responses to the extreme high temperatures which are more\nfrequent in the regions where the majority of dengue virus transmission occurs\n\\citep{Thomas2012,Damtew2023}. Although these limitafions have not prevented our model from\nproducing accurate predicfions of populafion dynamics in the regions we consider, unless addressed\nthey may limit the generalisability of our predicfions of populafion dynamics to future climates, a\nfactor that could be mifigated through modificafions to standard experimental protocols.”\nMinor\nParagraph lines 59-70 (and a few later examples such as line 176) – strictly, it is ‘dengue virus’ which\nis vectored by Ae. albopictus, not ‘dengue’, which refers to the disease itself in humans.\nThis is now corrected, and “dengue” and “dengue virus” are now used precisely throughout.\n(Remarks on code availability):\nBeyond my skillset\nReviewer #2 (Remarks to the Author):\nIn this manuscript, the authors use a SIR transmission model with mosquito populafion dynamics to\ninvesfigate the impact of phenotypic plasficity in vector traits on dengue transmission in the context\nof Aedes albopictus. The mosquito populafion dynamics sub-model is more complex than those\ntypically encountered in the epidemic modelling literature and has been introduced in a previous\nwork. An important feature of the present model is its ability to yield a realisfic dynamics of\nmosquito wing-length, which the authors use as a proxy for mosquito size and lifespan. This trait\nvaries over fime due to climafic and ecological factors, resulfing in subtle dynamical changes in\nepidemic risk. Overall, I believe that the manuscript is well wriften and that the implicafions from this\nmodel are potenfially very important. I parficularly appreciated several of the considerafions offered\nin the text. However, I believe that the implicafions for epidemic modelling are not sufficiently\narficulated in the text. Please find my detailed comments below:\n- The authors focus on the case of dengue and Aedes albopictus. Not being an expert in mosquitos, I\nfind it unclear whether similar reasonings apply to other host-pathogen systems as well. Is\nphenotypic plasficity equally variable/relevant to epidemic spread in other mosquito species? I\nwould thus suggest the authors to comment on the relevance of their findings to other systems.\nThis same mechanism of phenotypic plasficity altering vector competence has been proposed in\nother vector species (see Cator et al. (2023)). We further discuss the generalisability of our findings\nto Aedes aegpyfi, the primary vector of dengue, in a new discussion paragraph. (L332-349)\n“Dengue transmission frequently occurs in regions where Ae. albopictus co-occurs with the primary\nvector of dengue, Ae. aegpyfi, and with which it competes for resources in larval habitats\n\\citep{Juliano1998}. This complicates our ability to predict the effect of trait dynamics on\ntransmission efficacy as the relafionship between adult wing length, temperature, larval density, and\nfood availability is asymmetrically altered by the degree of interspecific compefifion that individuals\nof either species experience \\citep{Juliano2004,Murrell2008,Lizuain2022}. Further, in endemic\nregions co-circulafing serotypes and previous exposure to infecfion complicate transmission\ndynamics by introducing an immune structure to the human populafion \\citep{Gibbons2007}. An\nextension of our model to these areas is possible but would addifionally require the development of\na comparable model for Ae. aegypfi, to quanfify the effect that interspecific compefifion has on the\nlife-history of both species and then to extend the disease dynamics to consider mulfiple co-\ncirculafing serotypes of dengue virus. However, even in these more complex transmission scenarios\nthe underlying interacfion between environmental drivers and vector trait expression is driven by the\nsame mechanisms of phenotypic plasficity that we explore here, as they are a fundamental feature\nof vector ecology. By developing our understanding of the drivers of disease dynamics that arise in\nregions with a single vector transmifting a single serotype of dengue virus, we are now well\nposifioned to begin addressing the more complex dynamics underlying pafterns of dengue\ntransmission in other regions.“\nEventually, the fitle could be made more specific as in its current form, the manuscript focuses\nheavily on a single applicafion and it remains unclear how these result generalise.\nWe agree that the fitle should be more specific and have changed the manuscript fitle to:\n“Phenotypic plasficity in vector traits drives trends in global disease incidence: the spread of dengue\nvirus by Aedes albopictus”\n- The model appears to be accurate at reproducing mean wing length data whenever such data is\navailable (Fig 1). However, Fig 1 and Supp. plots show that such data is often available within short\nfime intervals. Moreover, wing-length does not vary a lot during these periods, whereas the model\nseems to suggest that this trait can be much more variable over a longer fime window. Given that\ndynamical variafion in wing-length underpins much of downstream results, I wonder if the authors\ncould provide further evidence to support the “U” shape displayed by wing-length? At present, I\nfear that the model is extrapolafing and I am not sure whether these extrapolations are jusfified or not.\nIn general, the relafionship between high larval compefifion/temperatures and a reducfion in wing-\nlength is well quanfified through laboratory experiments in Aedes albopictus and a range of other\nmosquito species (see for example Ezeakacha and Yee (2019)). The “U” shape is an emergent\ndynamical behaviour that arises from the reacfion norms we have stafisfically fifted as a\nconsequence of seasonal variafion from low->high->low larval compefifion and low -> high -> low\ntemperatures. This is generally the behaviour that is observed in field populafion of larvae and\ntherefore qualitafively is as we would expect.\nHowever, we agree with the reviewer that as this is such a crifical part of our argument that we\nshould back this up with addifional field data. Unfortunately, the availability of seasonal wing-length\ndata for Aedes albopictus is generally limited and the few datasets we present are the extent of the\ndata we are aware of. However, for other mosquito species we can observe similar U-shaped trends\nin seasonal wing-length. For example, Foley et al. (2020) demonstrate a U-shaped wing-length\ndynamics in four different Anopheles vectors of malaria in Korea, the same paftern of variafion is\nobserved in various Culex species in Japan, the USA, and Australia (where the relafionship becomes\nan “n” to match the Austral winter) , and for Aedes species in the UK. A subset of these examples is\nnow presented in the text to provide addifional context for our results. (L165-169)\n“Wing length is a trait that is observable in the field where it is known to seasonally vary in a manner\nthat qualitafively matches with the variafion we have observed in populafions of Ae. albopictus and\nother mosquito species \\citep{Willis1994,Packer1989,Foley2020,Kassim2013}. This variafion is highly\ncorrelated with variafion in adult longevity \\citep{Ezeakacha2019}.“\n- The authors show, rather convincingly, that different modelling assumpfions yield quanfitafive\ndifferences in epidemic dynamics. In parficular, the authors show that dengue incidence and aftack\nrate can vary considerably depending on how heterogeneity in wing-length is accounted for (or not\naccounted at all). However, while I do not have any doubts that the presented model is more\nbiologically plausible than other alternafives, I’m not enfirely convinced about its advantages in\nepidemiological modelling (e.g. in real-fime tracking, planning for vector control). The pracfical\nsuperiority of the full model is also difficult to assess without proper model comparison: it could well\nbe that simpler models fit the data equally well or befter, but we don’t know. Here I’m not suggesfing\nto perform a proper model selecfion exercise, but it would be good if the authors expanded on the\nimplicafions of their eco-epidemiological model over simpler models. For example, the considerafion\nthat epidemic risk varies over fime due to wing-length (and that mosquito abundance alone is not\nsufficient to predict transmission) is extremely interesfing in my opinion. Potenfial avenues of\ndiscussion include implicafions for vector control and disease persistence during inter epidemic\nperiods. Also, what are the consequences for R0 esfimafion? If I had to fit a simpler model to data\n(case data, serological data), would I get a wildly different esfimate of R0 with respect to using the\nfull model? How would these esfimates then affect what we know about disease management?\nThe constant wing-length model developed in Supplementary Informafion 3 is the closest analogy to\nthe simpler Ross-Macdonald style models that our framework can make. A comparison between the\nRt esfimates made in the full model with phenotypic plasficity and by the simpler constant wing\nlength model variants are shown in Supplementary Figure 56, which is now menfioned in the main\ntext (L225-226). We further show that the different model variants make different predicfions of Rt\nwhich then consequently result in different dengue dynamics. As the role of environmental variafion\nis so fundamental in driving the trait and disease dynamics it is difficult to make any more specific\nstatements (although this is an interesfing avenue for future explorafion).\nWe directly compare the spafial distribufion of risk predicted by our approach and a commonly used\nsimpler model (Mordecai et al. 2017) in Supplementary Informafion 4 (referenced in the main text\n(L247-250)). In general our approach predicts risk in fewer places for a shorter period of fime than\nprevious approaches. This is exemplified In Europe, where we can be certain that dengue is caused\nby albopictus, and we predict risk over a more limited geographic region than the previous approach.\nDespite this our model predicts 38 of the 41 clusters of autochthonous dengue transmission that\nhave occurred in Europe in the 21st century (L264-266).\nFurthermore, we now discuss the use of our accurate temporal predicfions of both vector and\ndengue dynamics, in applicafion to vector control in a new discussion paragraph (L312-331).\n“When designing vector control campaigns to implement intervenfions such as the sterile insect\ntechnique, gene-drive or transgenic control, it is common to rely on mathemafical models of vector\npopulafion dynamics \\citep{Ogunlade2023}. These models are used before implemenfing control\nmeasures to inform the required level of populafion suppression and the intensity of control\nmeasures necessary to achieve this target \\citep{Oliva2021, Bier2022, Powell2022}. Similarly, after\ncontrol measures have been implemented, mathemafical models are a pracfical method to assess\nhow effecfive an intervenfion has been by predicfing how an outbreak would have progressed in the\nabsence of control \\citep{Bonds2012}. Without this addifional analysis field-trials can only directly\nlink the efficacy of their intervenfion to a reducfion in mosquito abundance and not necessarily the\nreducfion in human cases of disease that they aim to achieve, two factors we show here are not\nnecessarily directly linked \\citep{Balatsos2024,Zhang2024}. The 2014 outbreak of dengue in\nGuangzhou provides a well-studied example showing the potenfial of our approach as we and\nprevious models are able to use mathemafical models to assess the efficacy of the implemented\ncontrol measures (see Figure 2) \\citep{Tang2016,Lin2016}. However, for the majority of the\noutbreaks we considered, this sort of retrospecfive analysis is not possible due to a lack of detail in\nthe reporfing of the nature, fiming, and intensity of control measures implemented. Further, to make\nreliable predicfions to assist in the design of vector control it is imperafive to use models that have\nbeen independently validated against real-world data from the target system \\citep{Brady2016}.\nExtensive and independent validafion of mathemafical models of systems of VBDs is currently not\nstandard, despite the ready availability of surveillance data and is a crifical step in ensuring these\npredicfions are robust.”\n- How do esfimates of suitability compare with another widely used suitability index, “Index P”\n(Obolski et Al, Methods in Ecology and Evolufion, 2019; Nakase et Al, Scienfific data, 2023)? Does the\npresent index provide further insights on why Italy experienced an unprecedented dengue outbreak\nin 2023 (please refer to Branda et Al., Dengue virus transmission in Italy: surveillance and\nepidemiological trends up to 2023, medRxiv)?\nOur understanding is that “Index p” is a metric that esfimates the vector competence porfion of a\ntradifional R0 equafion (i.e. R0 with no mosquito abundance), and so what we say about mean-field\nbased approaches to esfimafing disease risk also directly applies to index P including the comparison\nin Supplementary Informafion 4 (see also comments in the main text L45-59 & 217-234).\nAlthough 2023 is outside our study fime-frame the areas of Italy where transmission was observed in\n2023 correspond to areas that we consistently predict are suitable for (relafively) small dengue\noutbreaks (See Figure 4). Our research group is currently exploring whether or not climafic changes\ncan explain the recent intensificafion of European outbreaks, but this is out with the scope of this\npiece of work. Our inifial impression is that in Europe, although the transmission season is limited,\nthe outbreak size is quite environmentally sensifive (see Alpes-Marifimes in Figure 2) and so it is\ncertainly possible that 2023 was excepfionally suitable.\nMINOR COMMENTS\n- L114 onwards: please briefly list which traits are being predicted.\nChanged to “wing-length dynamics” as this is the only trait predicted by the model that is directly\nvalidated against field-data. We do model other traits such as stage specific development rates and\nthrough-stage survival proporfions but as these are not roufinely measured in the field we cannot\nvalidate these predicfions beyond the cumulafive effect each of these have on the populafion’s\noverall dynamics. This is consistent with other environmentally sensifive VBD models (e.g., Mordecai,\nCaldwell, Metelmann etc).\n- Fig 2: make x axis labels uniform (2018-2019-…), increase line width…\nIn Figure 2 the legends have been made uniform and line-width increased.\n- L139 why complex?\nRemoved complex. As Reunion is the only locafion we consider where the outbreak is sustained\nconfinuously for more than a year, this outbreak is epidemiologically more complicated (resistant\npopulafions, interepidemic periods of low transmission etc.). However, we feel this is adequately\nconveyed by “mulfi-year” alone, as from the perspecfive of our model this outbreak is treated in the\nsame way as all of the other outbreaks (and indeed this is the point).\n- Please increase x,y fick label size where this is too small to read.\nAxis labels and legend sizes have been increased throughout the text.\n- Methods L23: “both phenotypes”: what does this refer to?\nChanged to “both diapausing and non-diapausing eggs” to increase clarity.\n- Methods L29: “of” should be “or”.\nTypo corrected.\n- Methods Eq.(3): perhaps a typo in “1(…)”?\nWe assume that hatching takes a day to complete and that only a proporfion of quiescent eggs\nundergo this process. Therefore h_q(t) is a rate but (1 – (V-W(t)/V)) looks inifially dimensionless\nbecause there is an implicit 1/day. The 1(…) was an aftempt at dimensional clarity but admiftedly\nlooks odd. We have changed the equafion to 1x(…) to increase clarity.\n- Methods L134: please add further details on how expressions like that for delta(T) are derived.\nThe relafionship between the mortality rate from through-stage survival and stage durafion is\ninherited from the underlying Nisbet and Gurney framework for stage-structured delay differenfial\nequafions as a consequence of how they “lump” the von-Foerster equafion into age classes. We feel\nthat the argument demonstrafing that this form of delta(T) is valid is too technical to reproduce here,\nbut we now briefly explain the provenance of the equafion in the text and provide a reference that\nsupplies the full mathemafical details. (L519-523)\n“From the underlying Nisbet and Gurney (1983) framework we inherit a relafionship between the\ndevelopment rate, the probability of through-stage survival and the mortality rate. We therefore use\nthe temperature dependent reacfion norms for development and through stage survival to define an\nexpression for the temperature dependent mortality rate of acfive eggs, $\\delta_{E_\\gamma}(T) = -\n\\log(\\hat{S}_{E_\\gamma}(T))g_{E_\\gamma}(T)$ \\citep{Nisbet1983}.”\n- Methods L329: I would recommend the authors to label their model as a SEIR instead of a SIR, since\nthe it does indeed account for the intrinsic incubafion period. It also feels a more appropriate\nlabelling choice for a disease such as dengue.\nWe agree and the model is now referred to as an SEIR model throughout and text added to make it\nclear how the intrinsic incubafion period is accounted for without the need for an explicit exposed\nclass. (L748-751)\n“The omission of an explicit state-equafion for humans that have been exposed to infecfion, but who\nare not yet infecfious is a consequence of the intrinsic incubafion period being implicitly accounted\nfor in the delay structure through the inclusion of the term $\\tau_{IIP}$.”\n- L351: “an infecfious humans”\nThis error has been corrected.\n- Methods L404: do you mean $R_{t’}?$\nThis error has been corrected as suggested.\n- Methods L406: “This Is accounts”\nThis error has been corrected.\n- Methods why relafive humidity is not included??\nHumidity is undoubtedly a crifical factor in determining the life-span of adult mosquitoes, and so\npotenfially crifical in determining the potenfial of vector populafions to transmit disease. However,\nhumidity is so rarely varied in laboratory seftings, that we were unable to rigorously parametrise\nreacfion norms describing the relafionship between adult longevity and humidity (see Schmidt et al.\n(2018)). This data gap is not unique to Aedes albopictus and is a barrier to considering the effects of\nrelafive humidity on the populafion dynamics of mosquitoes in general (Brown et al. (2023)). We are\ncurrently involved in an effort to design experimental work to obtain these parameters, but this is a\nsubstanfial undertaking.\n- Perhaps show populafion dynamics of mosquitoes and compare it with other populafion dynamics\nmodels?\n- Supplementary and Main: “Emilia-Romagna” and “Catania” are the correct names.\nTypo corrected in the main text and throughout the Supplementary Informafion.\n- Supplementary P8 “Through through”\nThis error has been corrected.\n- Supplementary P9 “decease”\nThis error has been corrected.\n- Supp Fig 7A: It is not clear to me if the y axis denotes total numbers of proporfions of adults.\nBoth Figure and legend updated for clarity, we now plot the number of infected adults per 4km2.\n- Supp Fig 10 legend “Simluated”\nThe Figure legend has been corrected.\n- S. 1.1.2 “year 2017”. Perhaps it should be 2012.\nThis error has been corrected to 2012.\n- Supp Fig 14: the y label is cropped out.\n- S2.4: is it reasonable to assume that the populafion in Reunion is fully suscepfible against DENV,\ngiven previous circulafion of the virus in that area?\nAlthough DENV had previously circulated at low levels in Reunion, the outbreak we consider here\nwas the first large outbreak on the island since 1978 (see Vincent et al. (2023)). In these previous\nyears DENV-2, the serotype responsible for the majority of cases in the outbreak we consider, was\nrarely detected (hftps://www.ncbi.nlm.nih.gov/pmc/arficles/PMC9333208/ ). Given both of these\nfactors, we feel it is reasonable to make the assumpfion that the populafion is fully suscepfible to\nDENV-2, as the proporfion of the populafion with resistance at the start of 2017 was likely low. This\nis now summarised in the supplementary text (S.2.4.).\n“Previous to the 2017-2021 outbreak the last large outbreak of dengue in on Reunion occurred in\n1978 \\citep{Vincent2023}. Further, in the years preceding this outbreak the dominant circulafing\nserotype was DENV-I. We therefore assume that the populafion of Reunion is inifially completely\nsuscepfible to the introduced dengue serotype.”\n- S3.2 “iin”\nThis error has been corrected.\n- In the maps, both in the main and Supp material, it would be good to highlight the locafions\nconsidered in this study.\nThe locafions of all field datasets considered are now provided in Figure 1."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have satisfied my initial concerns with the paper. In particular, demonstrating that they had not cherry-picked examples, and the addition to the discussion regarding limitations in Ae. albopictus-only areas needed for validation. I recommend the manuscript for publication.", "reviewer_2": "I am pleased with the authors' response and I would like to commend them for the huge amount of work that went into this manuscript."}, "author_response": ""}]

Decision letter