Population of the temperate mosquito, Culex pipiens, decreases in response to habitat climatological changes in future.

Predictions of the temporal distribution of vector mosquitoes are an important issue for human health because the response of mosquito populations to climate change could have implications for the risk of vector-borne diseases. To elucidate the effects of climate change on mosquito populations inhabiting temperate regions, we developed a Physiology-based Climate-driven Mosquito Population model for temperate regions. For accurately reproducing the temporal patterns observed in mosquito populations, the key factors were identified by implementing the combinations of factors into the model. We focused on three factors: the effect of diapause, the positive effect of rainfall on larval carrying capacity, and the negative effect of rainfall as the washout mortality on aquatic stages. For each model, parameters were calibrated using weekly observation data of a Culex pipiens adult population collected in Tokyo, Japan. Based on its likelihood value, the model incorporating diapause, constant carrying capacity, and washout mortality was the best to replicate the observed data. By using the selected model and applying global climate model data, our results indicated that the mosquito population would decrease and adults' active season would be shortened under future climate conditions. We found that incorporating the washout effect in the model settings or not caused a difference in the temporal patterns in the projected mosquito populations. This suggested that water resources in mosquito habitats in temperate regions should be considered for predicting the risk of vector-borne diseases in such regions. ©2017. The Authors.

Introduction

Climate change has several impacts on ecosystems, and the spread of vector‐borne diseases that threaten human lives and health is of considerable concern [McMichael et al., 2003]. Changes in temperature, rainfall, and humidity influence the scale of mosquito‐borne diseases in complex ways [Smith et al., 2014]. Therefore, adaptation to and preventive measures to mitigate possible situations arising from climate change are urgent issues worldwide [Smith et al., 2014]. However, appropriate preparation for future climate conditions is difficult because risks related to climate change differ according to a climatic zone and area characteristics.

The responses of mosquito populations to climate change represent a potential risk of vector‐borne diseases, especially in temperate regions. Although a few indications of reduced pandemic or epidemic vector‐borne disease risks owing to improved hygiene and medical care are noted in temperate areas, cities within these areas might have epidemics, because of the presence of the vector and high population density. In addition, the distribution area and growth rates of disease‐vector mosquitoes are likely to change in such cities owing to rapid temperature increases [Fujibe, 2009]. Thus, cities in temperate areas have high risk of spreading vector‐borne diseases under climate change conditions.

In this scenario, research forecasting the temporal distribution of mosquitoes in relation to the timing of mosquito control or focusing on disease preservation for calculating the risk of vector‐borne disease is considerably important. A statistical approach, i.e., the correlation between vector presence and environmental information, can illustrate vector's current geographical distribution and anticipate that under future climate conditions [Hongoh et al., 2012; Fischer et al., 2014; Ren et al., 2016; Samy et al., 2016]. This approach is reasonable to describe macroscale spatial vector distributions [Ren et al., 2016]. However, knowledge regarding process‐based growth models are required to understand the direct influence of climate on vector presence and to reveal the detailed temporal distribution of the vector, because the response of population growth to climate change is complex and cannot be predicted using a statistical model [Ogden et al., 2005; Hongoh et al., 2012]. Although few studies have modeled mosquito population growth under climate conditions in the future, Morin and Comrie [2013] predicted intra‐annual variations in Culex quinquefasciatus populations from several climatic regions by using the Dynamic Mosquito Simulation Model and future climate data. Their results suggested that under future warmer climate conditions, C. quinquefasciatus inhabiting the subtropical region would appear earlier and disappear later in the year, and population density would be lower in the summer than that under current climate conditions. However, Morin and Comrie [2013] targeted only the southern area of the United States, not providing general data on how mosquito populations respond to climate changes in other areas. Thus, broadening our knowledge regarding how climate change affects temperate mosquitoes is necessary.

Modeling mosquito population dynamics in temperate climate areas can generate detailed temporal variation data under the current climate conditions, but this has not been actively performed. Most studies concerning mosquito population dynamics focused on Aedes aegypti [Focks et al., 1993; Otero et al., 2006], Anopheles gambiae [Depinay et al., 2004; White et al., 2011], and C. quinquefasciatus [Morin and Comrie, 2010] from tropical and subtropical regions because they are the vectors of ongoing infectious diseases in these areas. In the present study, we focused on mosquito population dynamics in temperate regions. For constructing the population dynamics model for temperate mosquitoes, the distinctive competence (diapause and effects of rainfall and temperature) of these mosquitoes needs to be considered.

In temperate regions, diapause is a key factor that maintains mosquito populations throughout the year. Nevertheless, incorporating this into a model is difficult because of the lack of field observations or laboratory data. Cailly et al. [2012] improved Anopheles spp. population model by incorporating the effect of diapause based on limited observation data; further, the diapause effect was incorporated into models for Aedes albopictus [Tran et al., 2013; Jia et al., 2016] and Culex pipiens [Ezanno et al., 2015; Marini et al., 2016]. However, the population dynamics incorporating diapause are not completely understood because few studies focus on diapause characteristics in each mosquito species, including photoperiod and temperature thresholds [e.g., Oda, 1971; Hawley et al., 1987; Pumpuni et al., 1992; Lacour et al., 2015]. Moreover, these characteristics depend on not only the mosquito species but also latitude and/or growing environment [Oda et al., 1988]. Thus, the mosquito population model incorporating the effect of diapause can be improved.

To predict mosquito population dynamics under future climate conditions, the relationships between mosquito population development and temperature and rainfall must be taken account when constructing the model. Mosquito growth rate and survival are highly affected by thermal conditions [Ikemoto, 2005; Loetti et al., 2011; Oda et al., 1999]. In addition to temperature, water plays a pivotal role in controlling habitat conditions for juveniles [Dieng et al., 2012]. Several researches suggested that larval abundance largely depends on rainfall because it creates or expands their habitat [Mori and Wada, 1978; Bomblies et al., 2008; Galardo et al., 2009]. Conversely, recent studies have revealed the negative effects of heavy rainfall for eggs and larvae, which might be washed out from their aquatic habitats owing to overflow [Paaijmans et al., 2007; Dieng et al., 2012]. Furthermore, Collins et al. [2013] suggested that rainfall quantity and pattern might change in the future, altering juveniles' survival. Therefore, considering water resources and temperature at a habitat is essential for estimating population dynamics under climate change.

Although the influence of habitat's water resources on mosquito populations is unclear, some models consider it as a primary factor determining population growth [Shaman et al., 2002; Depinay et al., 2004]. Most studies regarding population models assume that the quantity of soil moisture or rainfall positively affects the carrying capacity of populations during their aquatic stage [Morin and Comrie, 2010; White et al., 2011; Tran et al., 2013]. In contrast, Cheng et al. [2016] considered the spillover effect from heavy rain. Although some studies based on mosquito population observations corroborated that water conditions highly affect mosquito populations [Mori and Wada, 1978; Bomblies et al., 2008; Galardo et al., 2009; Paaijmans et al., 2007; Dieng et al., 2012], the mechanism by which it operates remains poorly understood. Thus, modeling this effect is a new challenge in mosquito population dynamics.

To reveal the consequences of altered mosquito population dynamics owing to climate change in temperate regions, we modeled mosquito populations by simultaneously using physiological ecology and climatology perspectives. We named this new model as the Physiology‐based Climate‐driven Mosquito Population model for temperate regions (PCMP). To construct our model based on realistic population dynamic data, we used available data from several studies concerning field observations and experimentation of physiological responses to the changes in temperature. In the present study, we developed a population dynamics model for temperate mosquitoes focusing on three factors: the effect of diapause, the positive effect of rainfall for the larval carrying capacity, and the negative effect of rainfall on aquatic stages (i.e., the washout mortality). We chose C. pipiens, which is the vector of the West Nile virus and is broadly spread across medium‐ and high‐latitude areas, as a typical case of a temperate mosquito. By modifying the several determinant factors, i.e., by using different combinations of the three factors (hypothesis), we produced temporal variations in mosquito populations and calibrated model parameters for each hypothesis by using weekly mosquito observation data. By comparing the outputs of each model validation, we identified the factors contributing the most to mosquito population dynamics. Next, we investigated how the washout mortality changed population dynamics by comparing mosquito populations based on the model incorporating the effect of washout mortality to that not incorporating this effect, at all population stages within a year. Lastly, we estimated the dynamics of mosquito populations in the future by using global climate model (GCM) data with two representative concentration pathway (RCP) scenarios. In addition, we revealed how water conditions affect mosquito population projections.

Data and Methods

Field Observations and Meteorological Data

Tsuda and Hayashi [2014] collected weekly from May 2003 to December 2013 at the National Institute of Infectious Diseases (35°42′15″N, 139°43′00″E), Tokyo, Japan. We used those mosquito field observation data for parameter calibration and model validation. Adult mosquitoes were captured from one site, and two traps were set at a different height from the ground; the ground height of one trap was about 7.5 m (near a tree canopy) and that of the other was 1.5 m (near shrubs) [Tsuda and Hayashi, 2014]. Culex pipiens gr., which includes C. pipiens pallens and C. pipiens form molestus, were collected over a 24 h period every week throughout the year by using dry ice traps [Tsuda and Hayashi, 2014]. The total number of observation days was 548 because observations were not conducted on some days. We used the total number of mosquitoes collected from May 2003 to December 2010 (n = 399) for calibration and from January 2011 to December 2013 (n = 149) for validation.

Current meteorological data were obtained from the Japan Meteorological Agency database [Japan Meteorological Agency, 2016]. Air temperature, precipitation, relative humidity, solar radiation, wind speed, and cloud amount data collected at the meteorological station (35°41′30″N, 139°45′00″E, 25 m in altitude, 4.4 km away from the mosquito observation site) that was nearest to the mosquito observation field at Tokyo were used. Actually, the mosquito habitat's condition differs within Tokyo urban area depending on microclimate of the observed location. However, observations at the same time at multiple points have not been obtained. Therefore, we assumed that the meteorological data on the weather station could represent the typical weather condition of the urban area of Tokyo. Photoperiod was calculated using the National Oceanic and Atmospheric Administration solar calculator [National Oceanic and Atmospheric Administration, 2016].

For estimating mosquito population dynamics under future climate condition, we used the climate projections of the Model for Interdisciplinary Research on Climate (MIROC5) [Watanabe et al., 2010]. MIROC well reproduced precipitation from the seasonal rain front in East Asia [Kusunoki and Arakawa, 2012]; its latest version, MIROC5, was used in this study. The average daily values of air temperature, relative humidity, solar radiation, and cloud amount, and daily cumulative precipitation were converted to the units used for the data collected at the study point by using the Inverse Distance Weighted interpolation method, because the original data were under the mesh climate value format. We used the historical scenarios for the period 1990–2009 and RCP 2.6 and RCP 8.5 scenarios for the period 2080–2099. These scenarios consider different emission concentrations: RCP 2.6 was set as the low emission scenario maintaining temperature rise below 2°C, and RCP 8.5 as the high emission scenario.

Model Framework

The overall flow of the simulations performed in this study is shown in Figure 1. First, we calculated the microclimate conditions at mosquitoes' aquatic habitat for accurately describing their growth during the water temperature‐dependent stages. The mosquito population model was then developed using those data and other simple meteorological data. In the model, the population size was calculated based on the development and mortality rates depending on the daily thermal and water conditions.

Figure 1

The model framework used in this study that includes the three models’ conceptual diagrams. (right) The mosquito population model is the main model simulating population dynamics, (left) whereas the other two models are needed to reproduce the typical habitat of juvenile phases by using observed meteorological data. Solid arrows outside the boxes represent the input and output data of each model. F and F correspond to mean long‐wave radiation from air and water, respectively.

Calculation of Habitat Environmental Conditions

The mosquito life cycle consists of four stages: three aquatic (eggs, larvae, and pupae) and one aerial (adults; Figure 1). The eggs of C. pipiens are deposited in masses on the water surface of puddles, ponds, rice fields, etc. Larvae hatch in the water and develop into pupae, which then metamorphose into adults. Female adult mosquitoes suck animal blood to gestate eggs and to move to the oviposition site. Most adults emerged in late summer shift to diapause adults with undeveloped follicles that do not suck blood and stay motionless in warm spots over the winter [Oda et al., 1988].

For characterizing the aquatic growth environment, we calculated water temperature, soil moisture ratio, and runoff by using basic climate factors (Figure 1). The daily water temperature was estimated using a simple energy balance model following the method of Ohta et al. [1993]. They developed a physical model for describing the thermal environment of ponded shallow water based on the assumption that water temperature is initially affected by the partition of solar energy between air, water, and soil (Figure 1 and supporting information Text S1). Kashiwada and Ohta [2010] applied this method to simulate mosquito's development because the environment of mosquito juveniles' habitat is almost identical to that of ponded shallow water. Furthermore, Ohta and Kaga [2012] used the water balance model proposed by Tao et al. [2003] for calculating soil moisture in mosquito habitat. Soil water content was calculated as the difference between water input (precipitation and snow melt) and output (evapotranspiration), and runoff was defined as the overflowing water, assuming soil as a bucket (Figure 1). In the present study, soil water content and runoff were calculated following Ohta and Kaga [2012] (Figure 1 and Text S2).

Population Dynamics Model

In the PCMP model, we assumed C. pipiens gr. had five life stages: three aquatic [eggs (E), larvae (L), and pupae (P)] and two aerial stages [active adults (A) and diapausing adults (A dia)]. We formulated the population dynamics model as follows (All parameters used in this section are listed in Table S1–S3.):

where d (i = E, L, and P), m (i = E, L, and P), m , K, o , and z (j = 1, 2) represent the temperature‐dependent development rate, temperature‐dependent mortality rate, mortality rate that reflects washout individuals due to heavy rain, carrying capacity of larvae, temperature‐dependent oviposition rate, and parameters for adult diapausing, respectively. Adult mosquitos shift from A to A dia with z 1 and from A dia to A with z 2. These parameters are day length dependent as described below. The model describes the female population only. We assumed that 1% of active adults were captured in traps at observation.

The relationships between temperature and development/mortality rates were calculated, and their fitting curves were obtained by applying the temperature‐dependent functions and parameters from Ewing et al. [2016] (Text S3 and Figure S1). Because larval and pupal mortality rates were combined as the reciprocal of longevity, we split them based on the length of each growth stage. Development/mortality rates were calculated using daily air or water temperature [T or T , respectively (°C)], depending on the growth stage (Figure 1).

Our model implemented diapause, mosquitoes' essential characteristic for surviving winter season in temperate regions, as parameters z 1 and z 2 , which represent the probability of shifting from active to diapause adults and that of diapause breaking, respectively. To simplify the model mechanism, we assumed that z 1/z 2 depended only on photoperiod, which regulates diapause. These parameters were formulated as sigmoid curves as follows:

where α and β (i = 1, 2) are the slope and the half‐saturation point of the sigmoidal curve for the photoperiod threshold, respectively, and D week is the average sunlight duration (h) during the most recent 7 days. We assumed that diapausing adults never died during diapause, for simplicity.

The effect of habitat water conditions on the mosquito population was considered in larval carrying capacity and washout mortality. We assumed that soil moisture altered larval carrying capacity because it determines larvae habitat space. In the model, K is defined as

where κ max and W represent the maximum carrying capacity and the soil water content (mm), respectively, and W* refers to the soil water‐holding capacity (mm), which reflects the effects of soil texture, soil organic content, and plant root depth, obtained from Dunne and Willmott [1996]. Soil moisture ratio (%) is derived from W/W*. Washout mortality (m ) was assumed as the additional mortality of aquatic stages (eggs, larvae, and pupae) by spilling over from water surface owing to soil water runoff caused by heavy rain and was defined as

where r and γ represent soil runoff (mm) and daily runoff threshold (mm) causing washout mortality, respectively.

Model Calibration

Based on the several factors controlling mosquito population dynamics, we constructed eight different models. In the present study, we assessed the combination of the effects of diapause and habitat water conditions, which affected larval carrying capacity and washout mortality. Each factor implemented in the models is shown in Table 1. Herein, models are designated by three capital letters (e.g., DCW model): the first denotes the effect of diapause, where “D” means the model incorporates diapause and “N” means it does not; the second denotes larval carrying capacity, where “C” means carrying capacity was set to a constant value (K = κ max) throughout the simulating period and “V” means it varied according to soil water content; the last refers to washout mortality, where “W” means the model incorporates washout mortality and “N” means it does not. For each model, we conducted parameter calibration by using observation data and the simulated annealing method (see below). We then compared the reproducibility of each model by using validation data and selected the most suitable one based on its goodness of fit.

Table 1

Log Likelihoods of Each Model Obtained During the Validation Period by Using the Calibrated Parametersa

Carrying Capacity	Washout Mortality	Diapause
		N	D
C	N	−4241	−2766
C	W	−3987	−2712
V	N	−4613	−3155
V	W	−4292	−3082

Capital letters stand for constant (C) or variable (V) carrying capacity and including washout mortality (W) or diapause (D) in the model or not (N).

For parameter calibration, we optimized unknown parameters (indicated in Greek letters) by using simulated annealing (SA) [Kirkpatrick et al., 1983] with Markov Chain Monte Carlo (MCMC) sampling. SA is a general‐purpose optimization technique applying MCMC sampling, which is a method of obtaining a posterior distribution of a parameter. By gradually lowering the “temperature state” from high to low (annealing treatment), reaching the ground state (mode of posterior distribution) is possible while conducting a global search within the parameter space. “Temperature” refers to the concept of statistical physics, and “high temperature state” refers to the state that the difference in likelihood between different parameter values is scaled small (the algorithm of SA was described in detail in Kirkpatrick et al. [1983] and Goffe et al. [1994]). Although the three factors considered in the models are indispensable to explain the dynamics of mosquito populations in temperate areas, their detailed mechanisms have not been investigated because of the difficulty of observing their action in real populations. In parameter optimization procedures, as the objective function, we set the likelihood of the observed data over calibration periods, assuming that for each week, the number of observed trapped adult mosquitoes followed a Poisson distribution with weekly averaged population obtained from the model. Optimal parameters were then determined by heuristic estimation. All unknown parameters were kept constant during the simulation year.

Model parameters were calibrated using weekly observation data from May 2003 to December 2010 (n = 399), and the likelihood of each model obtained during the validation period (January 2011 to December 2013) was calculated using those calibrated parameters. Simulations started 1 year before the calibration period with an initial population comprising 100 active and 100 diapausing adults, and first year simulations were iterated 10 times so that they were not influenced by the initial state. We used the Euler method for the numerical calculation with a time step of 1/100 days. During the calibration and validation periods, populations were continuously calculated. The population dynamics of each growth stage and the microclimate data for 2011 were output to confirm the effects of diapause and water condition and to reveal the relationships between climate conditions and population dynamics.

Mosquito Population Dynamics Under Future Climate Conditions

Predictions were performed based on parameters calibrated as explained above. Historical and future mosquito population dynamics under each scenario were estimated for 19 year periods, excluding the first year for iterations. Wind speed was assumed as constant throughout the simulation periods (2 m s−1) and day length data were from 2010. We subtracted weekly averaged captured active adult population from that of a week before to define the active period, which was assumed to start when the difference became smaller than −1 after winter. Similarly, the active period was assumed to finish when the difference become smaller than 1 in autumn, after the active adult population fell below half of the peak population.

Results

Model Selection

Table 1 shows the likelihood of each model, revealing DCW as the best. Models incorporating diapause had higher log likelihoods than that without diapause (DXX versus NXX in Table 1, with X representing any condition in the other two factors). Similarly, models considering washout mortality had higher likelihoods than that not considering it (XXW versus XXN), and models with constant carrying capacity also performed better than that with variable carrying capacity (XVX versus XCX). We also calculated r value for the model from the relation of the estimated average weekly population against the observed population (Text S4). As the result, DCW model had the highest r value and was statistically significant with a significance level of 5% (Table S4). Thus, including diapause and washout mortality effects improved model accuracy, whereas variations in larval carrying capacity caused by rain had a negative effect. The time series trends in field‐observed and simulated mosquito populations are shown in Figure 2. The appearance/disappearance patterns estimated by the DCW model reproduce those of field‐observed data, although some points were underestimated during the active season (Figure 2a). The population dynamics simulated by the three different models are shown in Figure 2b, and the optimized parameters of each model are listed in Table 2. To reveal the effect of the calibrated parameters on the mosquito population dynamics, we performed sensitivity test for the parameters; their results are shown in supporting information (see Text S5). The major difference among models was reflected in the carrying capacity parameters, since it takes higher values with increasing model factors. In DCN and DCW models, population dynamics showed a unimodal pattern within each year, whereas the model not considering the effect of diapause (NCN) resulted in a bimodal pattern. Conversely, population peaks were higher in models including diapause than in NCN. Although DCW and DCN revealed similar patterns, DCW performed better than the model without washout mortality, i.e., DCN, based on their likelihoods in weeks 22–25 and 30–41 (Table 3).

Figure 2

Population dynamics of active adult mosquitoes estimated using the three models and observed in the field. Model outputs from 2003 to 2010 were used for parameter calibration and those from 2011 to 2013 were used for model validation (Table 1). (a) The solid line represents the active adult population obtained from the DCW model, and the inverted triangles represent weekly field data. (b) Differences among model outputs as time series trends; the grey solid line is the same as the solid line in Figure 2a, which is the model selected to replicate the field‐observed data, i.e., the DCW model. The black solid and dotted lines represent the active adult population dynamics obtained from NCN and DCN model estimations, respectively. The three capital letters defining model acronyms represent diapause inclusion (D) or not (N), constant carrying capacity (C), and washout mortality inclusion (W) or not (N).

Table 2

Calibrated Parameters for the Three Models Considering Constant Carrying Capacity (C) and Including Diapause (D) and Washout Mortality (W) or Not (N)a

Model	κ max (100–106)	α 1 (10−6–106)	α 2 (10−6–106)	β 1 (7–18)	β 2 (7–18)	γ (0–200)
NCN	101.07	–	–	–	–	–
DCN	101.37	105.66	100.35	13.90	16.71	–
DCW	101.50	105.08	100.37	13.84	16.61	0.01

κ max, maximum carrying capacity; α 1 and α 2, slopes of the sigmoidal curve for the photoperiod threshold; β 1 and β 2, half‐saturation points of the sigmoidal curve for the photoperiod threshold; γ, daily runoff threshold causing larval mortality. Numbers below parameters represent the search range for parameter optimization.

Table 3

Log Likelihoods Calculated Every 4 Weeks During the Period of the Appearance of Active Adultsa

Model	Week
	18–21	22–25	26–29	30–33	34–37	38–41
DCN	−54.1	−631.6	−237.2	−162.1	−162.3	−138.9
DCW	−55.9	−582.5	−232.1	−137.7	−146.7	−128.2

Both models incorporated diapause and constant carrying capacity (DC), differing only in the inclusion of washout mortality (W) or not (N). Values are the sums of those obtained for the same 4 weeks of each year (2011–2013).

Mosquito Population Dynamics Under Current Climate Conditions

The active adult populations from model outputs began to increase in early May, peaked at the end of July, and were present until November or December (Figure 2a). Field observation data revealed that active adults were present from the end of March until November, peaking in mid‐July [Tsuda and Hayashi, 2014]. Although the number of field‐observed individuals varied widely throughout the year, the patterns evidenced in model simulations (Figure 2b) were almost identical to those observed during field observation (Figure 2a). However, the number of individuals in peak populations was lower in model simulations than in field observations throughout the calibration and validation periods.

Climatic feature variations during 2011 were also evaluated in the present study (Figure 3). The field observation point was characterized by a cold winter and warm summer (the lowest average daily temperature is about 2°C and the highest is about 32°C), and yearly cumulative precipitation ranged from 1330 to 1860 mm, occurring more often from June to October due to the seasonal rain front crossing Japan in June and September (Figure 3a). Soil water runoff mainly occurred in June and September along with continuous rainfall or large amount of daily precipitation (Figure 3b), and thus, the calculated soil water content, i.e., soil moisture ratio, never reached 0% throughout the year (Figure 3b).

Figure 3

Climatic features at the field observation point in Tokyo during 2011. (a) Air temperature (T ), water temperature (T ), and precipitation (bars), (b) and soil moisture ratio (black line) and runoff from soil (bars), calculated from other meteorological data. T and T are shown as daily averages; precipitation, soil moisture ratio, and runoff are shown as daily values.

A part of the simulated population for each growth stage is shown in Figure 4. Although fluctuations at each stage were synchronized, they were larger for aquatic stages than for adult stages. Active adults emerged in early April; their population grew from early May to late June and then temporarily decreased from mid‐July to early August. Diapausing adults existed throughout the year, but their numbers were reduced as those of active adults increased in May. The diapause period, when z 1 was set to >0, occurred from the 219th day of the year (early August) to the 130th day of the following year (early May), and the diapause‐breaking period, when z 2 was set to >0, occurred from the 34th day (early February) to the 317th day (early December). The diapausing adult population was always larger than the active adult population. No significant differences were noted between diapausing and diapause‐breaking periods across simulation years, because diapause depends only on sunlight duration.

Figure 4

Simulated population dynamics for (a) eggs, (b) larvae, (c) pupae, and (d) active and diapausing adults in 2011. Black and grey lines represent the output of the DCW and DCN models, respectively. Blue bars in Figures 4a–4c represent the day of washout. Both models incorporate the effect of diapause (D) and constant carrying capacity (C), differing only in the incorporation of the washout effect (W) or not (N).

Washout occurred mainly in June and September along with soil water runoff (Figures 3b and 4a–4c). Populations at the stages affected by washout were reduced to almost half. Washout occurrences throughout the year ranged from 9 to 38. In the model considering the washout effect (DCW), heavy precipitation in May, June, and September reduced eggs and pupae populations (Figures 4a and 4c), and consequently, egg and pupa populations were larger in the DCN than in the DCW model output. However, during the summer, i.e., in July, August, and September, larvae populations showed an opposite pattern to egg and pupa populations (Figure 4b). For adult populations, the two models had similar patterns, although the estimated number of active adults differed between the two models (Figure 4d): the active adult population obtained by the DCW model (black line) increased less than the population obtained by the DCN model (grey line) in June, but increased more from late July to mid‐August.

Historical and Future Predictions Based On GCM Data

Simulations performed using GCM data (Figure 5) revealed that during the historical period (1991–2009), C. pipiens were always present from the end of April to the end of January (Figures 5a and 5c), similar to the trend obtained with current meteorological data. Population peaks in the historical simulation were larger, and active periods were longer than those obtained with current climate data. However, in future predictions (2081–2099) under the RCP 2.6 and RCP 8.5 scenarios, mosquito populations decreased throughout the simulated periods (Figures 5b and 5d), the active period shortened (Figure 6), and population peaks occurred earlier than in historical simulations (Figure 5, Table 4, and Figure 6). These trends were more evident in RCP 8.5 simulations, where active adults disappeared considerably earlier than in other simulations, although they appeared at a similar timing (Figures 5a and 5b and 6).

Figure 5

Active adult populations simulated using global climate model data. (a and b) The results of the DCW model and (c and d) the results of the DCN model. Figures 5a and 5c represent historical (1991–2009) time series, and Figures 5b and 5d show the projections for future population dynamics (2081–2099) under two different scenarios: blue lines, representative concentration pathways (RCP) 2.6; red lines, RCP 8.5. Both models incorporate the effect of diapause (D) and constant carrying capacity (C), differing only in the incorporation of the washout effect (W) or not (N).

Figure 6

Length of adults' activity periods under the three climate scenarios investigated in this study. The grey, blue, and red lines indicate the results of historical (1981–1999), representative concentration pathways (RCP) 2.6, and RCP 8.5 future projections (2081–2099), respectively.

Table 4

Simulated Appearance Periods of Active Adults Based on Global Climate Model Data by Using Models Incorporating Diapause and Constant Carrying Capacity (DC) and Differing Only in the Inclusion of Washout Mortality (W) or Not (N)a

Model	Scenario	Start (week)	End (week)	Peak (week)	Peak Captured Population (number of individuals)
DCN	Historical	18.0	42.9	28.7	92.8
	RCP 2.6	19.0	39.5	27.9	61.2
	RCP 8.5	20.0	36.0	27.1	28.6
DCW	Historical	19.0	41.1	29.7	63.9
	RCP 2.6	19.8	38.7	28.3	37.6
	RCP 8.5	20.1	35.6	27.0	19.6

Values are the averages from 19 year simulations under each scenario: Historical, (1981–1999) and future (2081–2099) representative concentration pathway (RCP) projections.

Future projections based on the DCN model (Figures 5c and 5d) showed an overall increase in the number of mosquitoes compared to that obtained with the DCW model in all scenarios (Figures 5a and 5b and Table 4). In addition, variations in population growth rates became larger when the washout effect was added to all scenarios (Figures 5a and 5b). The active period is shortened by delaying active adults' appearance and by accelerating their disappearance in both the model assumptions (Table 4). However, peak population decrease in relation to historical simulations did not differ between the two models for both RCP scenarios. Conversely, shifts in peak timing, which were accelerated in future predictions, were slightly larger in the DCW model than in the DCN model (Table 4).

Discussion

By projecting the temporal distribution of mosquitoes, we revealed the responses of mosquitoes to climate change to better understand the risk management of vector‐borne diseases in the future. The PCMP model incorporating diapause and water conditions (i.e., DCW model) reproduced the seasonal trends observed weekly at the study point throughout the simulated period, i.e., under current temperate climate conditions. Incorporating diapause greatly improved model predictions (Table 1). Considering water conditions, the effect of soil moisture content, reflected in larval carrying capacity, contributed little to improve model accuracy, whereas washout mortality was partly effective to reproduce mosquito population dynamics in the model incorporating it. This result strongly supports those of Paaijmans et al. [2007] and Dieng et al. [2012], who emphasized that rain spillover affected mosquito populations. The effect of soil moisture on larval carrying capacity was not strong, contrary to the positive effect referred in many studies [Morin and Comrie, 2010; Tran et al., 2013; Jia et al., 2016]. The increasing and decreasing population levels obtained were similar to that reported in other studies on population dynamics models [Ezanno et al., 2015; Marini et al., 2016]. However, a marked decrease in active adult mosquitoes was noted for populations under future climate conditions (Figure 6), remarkably differing from the predictions obtained in previous studies [Marini et al., 2016]. This difference sheds new light on the risk management of vector‐borne diseases in the future because the shift observed in the vector prevalence period leads to fundamental changes in adaptation and preventive measures.

The appearance/disappearance of mosquitoes throughout the year is regulated by diapause in nature [Oda, 1968]; therefore, the diapause effect introduced in the DCW model should control the period of adult mosquito activity. The NCN model revealed a bimodal pattern in population dynamics, as shown in Morin and Comrie [2013], and adults' active periods during the year were longer than those observed in the field. This pattern is unnatural for temperate mosquito populations, and thus the model failed to reproduce it (Figure 2b). The population obtained from the model without diapause re‐increased rapidly in autumn, i.e., when the temperature became cooler, allowing the population to grow after its decrease during the high summer temperature (Figure 3a). However, in nature, the population decrease is observed in temperate regions [Maeda et al., 1978; Tsuda and Hayashi, 2014]. By incorporating the effect of diapause into the model, the population increase during autumn is suppressed because some active adults shift to the diapause state, and thus, the model output becomes consistent with the population trends observed in the field (Figure 2). The increasing trend observed in simulations for diapausing adults during autumn was consistent with that shown by Ezanno et al. [2015] for the seasonal distribution of C. pipiens in France by using a population dynamics model incorporating diapause. Thus, incorporating the effect of diapause is essential for modeling mosquito population dynamics in temperate regions, even under current climate conditions.

In addition to the effect of diapause, our results suggest that water conditions in the mosquito habitat strongly influence mosquito populations in temperate regions. With regard to washout mortality, the population dynamics changed at each growth stage, since the washout effect suppressed the increase in eggs and pupae during May and June (Figures 4a and 4c). However, this was not the case for the larval stage, since, in the model incorporating the washout effect, the larval population became larger than that in the model not considering the washout effect (Figure 4b). This was because of the differences in larval carrying capacity, which was independently optimized for each model setting, i.e., whether washout mortality was incorporated or not (Table 2). Differences between active adult population dynamics in the two models reflected the overall effect of adding rainfall, a negative factor to the model (Figure 4d and Table 1). Therefore, the likelihoods in weeks 22–25 and 30–41 by DCW model were evaluated to be higher than those by the DCN model (Table 3). In fact, Paaijmans et al. [2007] reported high losses of Anopheles gambiae larvae due to rainfall in western Kenya, i.e., in the tropical region. Based on the observations in the Alpes‐Maritimes in France, i.e., under current temperate climate conditions, Lacour et al. [2015] noted that the rainfall at the end of A. albopictus active season was important to kill juveniles and remaining adults, though their target species has a different life history from that of C. pipiens. With respect to soil moisture content, most studies considering the effects of habitat's water resources [Morin and Comrie, 2010; Tran et al., 2013; Jia et al., 2016] revealed that high moisture had a positive effect on mosquito populations. In the present study, the no positive effect of variable carrying capacity might be attributed to the continuous presence of water resources in the studied habitat. The habitat where the eggs, larvae, and pupae of temperate C. pipiens populations developed was never exposed to drought, unlike the mosquito populations occurring in semiarid regions or in any areas having a strong dry period where water is a critical limiting factor for population growth (Figure 3b). Thus, the effect of rainfall on mosquito populations might be positive or negative depending on the climatic region. Therefore, fine‐scale tailored predictions need to be performed for each region. If the present method applies to another location with different climate or regional characteristics, following the procedure from calibration and model selection is required since the present method involves site‐specific parameterization.

By using two mosquito population models to replicate the current observation data and applying GCM data, our results indicated that under future climate conditions, active mosquito populations will decrease and that the active season will be shortened. This is because of an increase in the number of days exceeding the suitable temperature for growing, as suggested by experimental data [Oda et al., 1980; Loetti et al., 2011], and is consistent with the findings of Marini et al. [2016], who suggested a decrease in active C. pipiens populations under an average temperature increase of 2.0°C. However, Marini et al. [2016] found that C. pipiens were present for a longer period under warmer climate conditions. This is a noteworthy finding because shorting of the active season influences the indication of pandemic. Moreover, Marini et al. [2016] predicted that mosquito populations in future climate conditions would appear earlier than that observed under the current climate conditions because of the rapid population increase owing to increasing temperature during spring. In our simulations, mosquito populations were less likely to be affected by temperature changes during early spring because active adults gradually increased depending on photoperiod. In addition to not changing the population growth rate, the exposure to high temperature during summer lowered the total abundance, leading to an overall reduction of the active period.

The results of our study highlight the importance of incorporating the negative effects of individuals' washout in mosquito population models projecting future dynamics in humid temperate regions. The difference between DCW and DCN model projections were reflected in several parts of mosquito population dynamics based on GCM data. The washout effect led to a smaller active adults' population, a shorter active period, and a rapid shift in population peak timing, which appeared earlier than that in the model without washout (Figure 5 and Table 4). These results suggest that changes in the rainfall pattern would strongly affect mosquito population dynamics under future climate conditions. The effect of water resources is often ignored when modeling mosquito populations in humid regions, but the accuracy of our models greatly improved due to considering water resources. This raises our confidence in future projections.

The underestimation of population peaks under current climate conditions might be because observation errors are not considered or wrong assumptions regarding adults' habitat are made. In addition, underestimation of the occurrence season length in several years (2007, 2012, and 2013) can also be considered to be attributed to the uncertainty of observation. We assumed that the mosquito capture rate was 1% during the entire simulation period for simplicity, but a constant capture rate is unusual in real environments. Experimental confirmation of the range of uncertainty is needed by using outdoor trapping, and this has to be incorporated into the model, since the captured population might considerably change. Moreover, adult mosquitoes can move to a cool place such as tree shades, and thus, population decreases due to high temperature can potentially be mitigated. Increasing surveillance of adults' habitat and meteorological observation would be useful to better estimate mosquito populations.

Although diapause improved the model accuracy when replicating observation data, it might have manipulated mosquito populations because of the assumption that diapause depends only on day length. However, temperature should also be considered when modeling the diapause mechanism because both photoperiod and temperature are related to diapause timing [Sanburg and Larsen, 1972; Oda, 1971]. Oda [1971] mentioned that the effect of photoperiod on the reproductive activity of C. pipiens pallens was weakened under high temperature conditions. In addition, the assumption that all diapause adults never die until diapause breaking should be revised because the diapausing adult population was always larger than the active adult population. Performing a simulation by using the model incorporating the temperature‐dependent diapause would be necessary if the dependency of diapause timing and/or mortality is experimentally clarified because it should affect the future projection. Furthermore, we applied the optimized parameters to current observations to estimate future population dynamics and ignored the possible changes in these parameters. Physiological traits might evolve throughout generations leading to differences in population responses to climate change.

In the present study, the additional mosquito mortality due to washout after heavy rain affected population peak size. Thus, the increase in heavy rain frequency predicted in future climate conditions would reduce mosquito populations inhabiting humid temperate regions. This effect appeared in the historical simulations performed using GCM data, and the difference in population peaks between the DCN and DCW models was larger than that in the validation period (Figures 4a and 4c). This is because of the pattern of continuous rain generated by GCM, where washout is more likely to occur than when using current meteorological data. For the same reason, rainfall patterns such as heavy rain frequency and continuous rainfall need to be included in future GCM data, and the uncertainty in rain‐day frequency and amount of precipitation remains at the regional scale [Chen et al., 2014]. Thus, obtaining precipitation data that reproduce the intensity of rain under future climate conditions is necessary. Furthermore, simulations performed using the PCMP model coupled with several different climate scenarios and/or GCMs would reveal the range of uncertainties of climate change impacts on mosquito population dynamics.

Supporting Information S1

Click here for additional data file.