Increased Adult Aedes aegypti and Culex quinquefasciatus (Diptera: Culicidae) Abundance in a Dengue Transmission Hotspot, Compared to a Coldspot, within Kaohsiung City, Taiwan

The assumption that vector abundance differences might drive spatial and temporal heterogeneities in vector-borne disease transmission is common, though data supporting it is scarce. Here, we present data from two common mosquito species Aedes aegypti (Linnaeus) and Culex quinquefasciatus Say, biweekly sampled as adults, from March 2016 through December 2017, with BG-sentinel traps in two neighboring districts of Kaohsiung City (KC), Taiwan. One district has historically been a dengue transmission hotspot (Sanmin), and the other a coldspot (Nanzih). We collected a total 41,027 mosquitoes, and we found that average mosquito abundance (mean ± S.D.) was higher in Sanmin (Ae. aegypti: 9.03 ± 1.46; Cx. quinquefasciatus: 142.57 ± 14.38) than Nanzih (Ae. aegypti: 6.21 ± 0.47; Cx. quinquefasciatus: 63.37 ± 8.71) during the study period. In both districts, Ae. aegypti and Cx. quinquefasciatus population dynamics were sensitive to changes in temperature, the most platykurtic environmental variable at KC during the study period, a pattern predicted by Schmalhausen’s law, which states that organisms are more sensitive to small changes in environmental variables whose average value is more uncertain than its extremes. Our results also suggest that differences in Ae. aegypti abundance might be responsible for spatial differences in dengue transmission at KC. Our comparative approach, where we also observed a significant increment in the abundance of Cx. quinquefasciatus in the dengue transmission hotspot, suggests this area might be more likely to experience outbreaks of other vector borne diseases and should become a primary focus for vector surveillance and control.


Introduction
The yellow fever mosquito, Aedes aegypti (Linnaeus), and the southern house mosquito, Culex quinquefasciatus Say, are two of the most common urban mosquitoes in tropical and subtropical environments around the world [1][2][3]. Ae. aegypti is the main vector of dengue, chikungunya, and Zika in center in Taiwan and accounts for more than 70% of dengue cases in the whole country (Taiwan Centers for Disease Control). Two highly urbanized (i.e., with over 40% of land cover used for housing and roads) districts of KC-Sanmin and Nanzih-were selected for the study. Sanmin district has historically been a dengue transmission hotspot with a high annual incidence rate (209.5 per 100,000). In contrast, Nanzih has been a dengue transmission coldspot, where dengue incidence rate is about 47.9 per 100,000. The location of the districts and sampling locations within each district are shown in Figure 1, and coordinates and pictures of the places where mosquito traps were located are presented in Supplementary Table S1 [online only]. For each sampling location, we evaluated the percent land cover in a circular buffer of 500 m radius around each trap, the maximum recorded dispersal distance of Ae. aegypti in urban landscapes [64], and a radius were most host seeking Cx. quinquefasciatus occur in relation to their larval habitats of origin [65]. The percent of urbanized land surrounding the traps was (mean ± S.D.): 55.75 ± 5.32 in Sanmin and 49.61 ± 8.49 in Nanzih, values not significantly different according to a Welch's t test (t = 1.501, df = 8.408, p-value = 0.17). Moreover, landscape configuration, assessed through a Principal Components Analysis (PCA) [65,66] of the percent of land used for housing, roads, parks, and other purposes, described in Supplementary Figure S1 [online only], also showed no significant statistical differences between the two sites (t = 2.314, df = 7.195, p-value = 0.05289) when compared with a Welch's t test, a t test that corrects df considering variance differences between the groups compared [67]. Despite the similarities in land cover, the studied districts have different population density per square KM, with Sanmin (17,321 habitants/KM 2 ) more than doubling the population density of Nanzih (7130 habitants/KM 2 ).

Study Area and Weather Data
Kaohsiung City (KC) is a city with 2.8 million habitants and a density of 942 people/km 2 (Department of Household Registration, Ministry of Interior). KC is the main dengue epidemic center in Taiwan and accounts for more than 70% of dengue cases in the whole country (Taiwan Centers for Disease Control). Two highly urbanized (i.e., with over 40% of land cover used for housing and roads) districts of KC-Sanmin and Nanzih-were selected for the study. Sanmin district has historically been a dengue transmission hotspot with a high annual incidence rate (209.5 per 100,000). In contrast, Nanzih has been a dengue transmission coldspot, where dengue incidence rate is about 47.9 per 100,000. The location of the districts and sampling locations within each district are shown in Figure  1, and coordinates and pictures of the places where mosquito traps were located are presented in Supplementary Table S1 [online only]. For each sampling location, we evaluated the percent land cover in a circular buffer of 500 m radius around each trap, the maximum recorded dispersal distance of Ae. aegypti in urban landscapes [64], and a radius were most host seeking Cx. quinquefasciatus occur in relation to their larval habitats of origin [65]. The percent of urbanized land surrounding the traps was (mean ± S.D.): 55.75 ± 5.32 in Sanmin and 49.61 ± 8.49 in Nanzih, values not significantly different according to a Welch's t test (t = 1.501, df = 8.408, p-value = 0.17). Moreover, landscape configuration, assessed through a Principal Components Analysis (PCA) [65,66] of the percent of land used for housing, roads, parks, and other purposes, described in Supplementary Figure S1 [online only], also showed no significant statistical differences between the two sites (t = 2.314, df = 7.195, p-value = 0.05289) when compared with a Welch's t test, a t test that corrects df considering variance differences between the groups compared [67]. Despite the similarities in land cover, the studied districts have different population density per square KM, with Sanmin (17,321 habitants/KM 2 ) more than doubling the population density of Nanzih (7130 habitants/KM 2 ).  Only one weather station (Kaohsiung Weather Service Office; ID:46744 (120.30 • E, 22.57 • N)) ( Figure 1) exists within the limits of KC, so all the climate variables used in this analysis, including daily average temperature, daily average precipitation, and daily average relative humidity, were summarized from the hourly measurements of the weather station. The study period coincided with Insects 2018, 9, 98 4 of 16 one of the strongest El Niño events on record [66,67], and as consequence of this phenomenon, weather patterns were anomalous in KC. During the last quarter of 2016, precipitation was higher than usual in KC (Supplementary Figure S2A [Online Only]). Meanwhile, mean temperature (Supplementary Figure S2B [Online Only]) was hotter than usual during the whole study period, a pattern also observed for maximum (Supplementary Figure S2C [Online Only]) and minimum temperature (Supplementary Figure S2D [Online Only]). For the analysis, we employed daily data from 1 January 2016 to 31 December 2017, which was processed at biweekly intervals coinciding with mosquito sampling. For each time interval, we computed average values and the standard deviation and kurtosis of the above mentioned weather variables to understand the impacts of environmental variability on mosquito population dynamics [40]. The average distance between the weather station and mosquito traps was 12.13 KM (SD: 3.96 KM, Range: 8.1-16.9 KM). The resulting biweekly weather data can be observed in Figure 2. In Kaohsiung City, winters, from December to February, have low rainfall (Figure 2A), low relative humidity ( Figure 2B), and temperature also drops to a minimum ( Figure 2C). By contrast, during the summer, i.e., June to August, rainfall (Figure 2A), relative humidity ( Figure 2B), and temperature ( Figure 2C) reach their maximum values. With the exception of rainfall, which had a strong variation and which was followed in magnitude by relative humidity, the SD of all temperature variables was similar and low during the study period ( Figure 2D). Meanwhile, the most leptokurtic variable, i.e., with high kurtosis values which represent low variation around the mean [68], where rainfall and maximum temperature, which reached maximum values from November to March, while the most platykurtic variable, more uncertain around the mean value than the extremes of its distribution [69], was minimum temperature ( Figure 2E). Only one weather station (Kaohsiung Weather Service Office; ID:46744 (120.30° E, 22.57° N)) ( Figure 1) exists within the limits of KC, so all the climate variables used in this analysis, including daily average temperature, daily average precipitation, and daily average relative humidity, were summarized from the hourly measurements of the weather station. The study period coincided with one of the strongest El Niño events on record [66,67], and as consequence of this phenomenon, weather patterns were anomalous in KC. During the last quarter of 2016, precipitation was higher than usual in KC (Supplementary Figure  . For the analysis, we employed daily data from 1 January 2016 to 31 December 2017, which was processed at biweekly intervals coinciding with mosquito sampling. For each time interval, we computed average values and the standard deviation and kurtosis of the above mentioned weather variables to understand the impacts of environmental variability on mosquito population dynamics [40]. The average distance between the weather station and mosquito traps was 12.13 KM (SD: 3.96 KM, Range: 8.1-16.9 KM). The resulting biweekly weather data can be observed in Figure 2. In Kaohsiung City, winters, from December to February, have low rainfall (Figure 2A), low relative humidity ( Figure 2B), and temperature also drops to a minimum ( Figure 2C). By contrast, during the summer, i.e., June to August, rainfall ( Figure 2A), relative humidity ( Figure 2B), and temperature ( Figure 2C) reach their maximum values. With the exception of rainfall, which had a strong variation and which was followed in magnitude by relative humidity, the SD of all temperature variables was similar and low during the study period ( Figure 2D). Meanwhile, the most leptokurtic variable, i.e., with high kurtosis values which represent low variation around the mean [68], where rainfall and maximum temperature, which reached maximum values from November to March, while the most platykurtic variable, more uncertain around the mean value than the extremes of its distribution [69], was minimum temperature ( Figure 2E).

Mosquito Sampling
In this study, we collected mosquitoes with a total of 12 BG-Sentinel traps (Biogents, A.G. Regensburg, Germany). Traps were set in the two small area districts of Kaohsiung City: Sanmin (19.79 KM 2 ) and Nanzih (25.83 KM 2 ), as shown in Figure 1. Six traps were placed in each district, from 3 March 2016 to 26 December 2017, a number constrained by our limited budget and manpower but appropriate for the small area districts that we studied. Traps were baited with BG lure, a substance that mimics compounds released by the human skin [70]. Baits were replaced every two months following manufacturer recommendations [63,71]. Traps were operated continuously and mosquito samples were collected by laboratory staff on a bi-weekly basis. All traps were set in sheltered open household spaces, e.g., yards and/or garages, and household owners kindly provided electricity for trap operation. BG traps baited with BG lure were primarily chosen given their ability to capture Aedes (Stegomyia) spp., for example, Ae. aegypti and Aedes albopictus (Skuse) [62,72], but BG traps are known to capture a diverse set of mosquito species, including Culex (Culex) spp. [62,73,74] with medical importance such as Cx. quinquefasciatus in subtropical environments [63,73]. At the landscape level, all traps were located in households with similar environmental characteristics, in urban areas dominated by residential use. Mosquito samples were collected by the laboratory staff and mosquitoes were identified according to morphological characteristics presented in the taxonomic key by JC Lien [75]. Both Aedes spp. and Culex spp. were included in the analysis. Only female mosquitoes were considered in the study due to their importance for disease transmission.

Statistical Analysis
We analyzed average time series for Sanmin (dengue transmission hotspot) and Nanzih (dengue transmission coldspot) districts of KC. We analyzed average time series because a couple of the traps broke during the first year of the study, which we could not replace because of limited funding. When estimating the biweekly mean abundance, we considered all traps available each week, since there was no significant difference on the estimated mean abundance for each mosquito species and study location when these traps were excluded (Appendix A), thus ensuring mean estimates were not biased by including data from the traps that broke at the start of the study. For the time series analysis, we employed autoregressive (AR) models [76]. On seven occasions, it was not possible to collect samples on the exact date scheduled for sampling. Therefore, before the time series analysis, we used the lowess method to interpolate observations for the missing dates [77] in order to fulfill the assumption of data sampled at equal time intervals of AR models [76]. To fit the AR time series models, we followed a protocol described in detail by T Hoshi, Y Higa and LF Chaves [78], where null AR models are fit to mosquito abundance time series after selecting lags by inspecting the autocorrelation function (ACF) and then confirmed by inspecting the partial autocorrelation function (PACF) of the time series. Then, the time series of environmental covariates are pre-whitened using the coefficients of the null AR models. The resulting time series are used to estimate cross correlation functions (CCF), which are employed to select environmental covariates and their lags. The selected variables are then included in a full model for each mosquito time series, which is then simplified by a process of backward elimination based on the minimization of the Akaike Information Criterion (AIC). Prior to fitting each time series model, the mean of each environmental variable is removed in order to ease the interpretation of the time series model intercept as the mean value of the time series [23]. For further details about the methods, please refer to T Hoshi, Y Higa and LF Chaves [78] and LF Chaves, N Imanishi, and T Hoshi [23].

Results
We collected a total of 41,027 mosquitoes, over a total 746 biweekly trap collections. We collected a total of eight Aedes albopictus (Skuse), 2729 Ae. Aegypti, and 38,290 Cx. quinquefasciatus. In each one of the districts, we collected 4 Ae. albopictus during the full study period. In contrast,  Figure 3A,C) than Nanzih (Ae. aegypti: 6.21 ± 0.47 see Figure 3B; Cx. quinquefasciatus: 63.37 ± 8.71 see Figure 3D). Figure 3 also shows that starting in 26 August 2016 we lost two of the traps, reducing the number of traps to five at each site. Overall, the population dynamics of Ae. aegypti at Sanmin ( Figure 3A) and Nanzih ( Figure 3B) did not show major fluctuations, on average, through the study period. Nevertheless, some locations showed "outbreaks" in the sense, that often, sudden extraordinary changes in the abundance of Ae. aegypti were observed, which were larger in Sanmin ( Figure 3A) when compared to Nanzih ( Figure 3B).
Insects 2018, 9, x 6 of 15 quinquefasciatus: 63.37 ± 8.71 see Figure 3D). Figure 3 also shows that starting in 26 August 2016 we lost two of the traps, reducing the number of traps to five at each site. Overall, the population dynamics of Ae. aegypti at Sanmin ( Figure 3A) and Nanzih ( Figure 3B) did not show major fluctuations, on average, through the study period. Nevertheless, some locations showed "outbreaks" in the sense, that often, sudden extraordinary changes in the abundance of Ae. aegypti were observed, which were larger in Sanmin ( Figure 3A) when compared to Nanzih ( Figure 3B). In general, all time series were autocorrelated (Figure 4), and patterns of autocorrelation were mainly suggestive of first order AR processes, i.e., significant autocorrelation at one time lag, for Ae. aegypti in Sanmin ( Figure 4A) and Nanzih ( Figure 4B) and Cx. quinquefasciatus in Sanmin ( Figure 4C) but not in Nanzih ( Figure 4D). Using this information, first order autoregressive models were fitted as null models to all time series but Cx. quinquefasciatus at Nanzih ( Figure 4D). In general, all time series were autocorrelated (Figure 4), and patterns of autocorrelation were mainly suggestive of first order AR processes, i.e., significant autocorrelation at one time lag, for Ae. aegypti in Sanmin ( Figure 4A) and Nanzih ( Figure 4B) and Cx. quinquefasciatus in Sanmin Insects 2018, 9, 98 7 of 16 ( Figure 4C) but not in Nanzih ( Figure 4D). Using this information, first order autoregressive models were fitted as null models to all time series but Cx. quinquefasciatus at Nanzih ( Figure 4D). The PACF of the Ae. aegypti time series from Sanmin furtherly suggested a first order AR process ( Figure 5A). This time series was not associated with the mean ( Figure 5B) and SD ( Figure 5C) of any weather variable, but it was associated positively at lag 1 and negatively at lag 2 with minimum temperature kurtosis ( Figure 5D). By contrast, the PACF of Ae. aegypti in Nanzih did not show any significant autoregression ( Figure 5E) and like the pattern observed in Sanmin the Nanzih time series was not correlated the mean ( Figure 5F) and SD ( Figure 5G) of any weather variable. However, it was associated positively at lag 1 with the minimum, mean and maximum temperature kurtosis ( Figure  5H). Meanwhile, the PACF of the Cx. quinquefasciatus time series from Sanmin confirmed a first order AR process, with a significant correlation at lag 7 that was ignored when fitting the time series models ( Figure 5I). This time series was significantly associated at lag 1 with the mean ( Figure 5J) and SD ( Figure 5K) of temperature (mean, maximum and minimum) and relative humidity and rainfall kurtosis ( Figure 5L). The PACF of Cx. quinquefasciatus at Nanzih, by contrast, showed no significant correlation ( Figure 5M) and was associated with the mean ( Figure 5N) and SD ( Figure 5O) of rainfall with five biweeks of lag, with temperature (minimum, mean and maximum) at lag 0 ( Figure 5N), and one biweek of lag with rainfall kurtosis ( Figure 5P). The PACF of the Ae. aegypti time series from Sanmin furtherly suggested a first order AR process ( Figure 5A). This time series was not associated with the mean ( Figure 5B) and SD ( Figure 5C) of any weather variable, but it was associated positively at lag 1 and negatively at lag 2 with minimum temperature kurtosis ( Figure 5D). By contrast, the PACF of Ae. aegypti in Nanzih did not show any significant autoregression ( Figure 5E) and like the pattern observed in Sanmin the Nanzih time series was not correlated the mean ( Figure 5F) and SD ( Figure 5G) of any weather variable. However, it was associated positively at lag 1 with the minimum, mean and maximum temperature kurtosis ( Figure 5H). Meanwhile, the PACF of the Cx. quinquefasciatus time series from Sanmin confirmed a first order AR process, with a significant correlation at lag 7 that was ignored when fitting the time series models ( Figure 5I). This time series was significantly associated at lag 1 with the mean ( Figure 5J) and SD ( Figure 5K) of temperature (mean, maximum and minimum) and relative humidity and rainfall kurtosis ( Figure 5L). The PACF of Cx. quinquefasciatus at Nanzih, by contrast, showed no significant correlation ( Figure 5M) and was associated with the mean ( Figure 5N) and SD ( Figure 5O) of rainfall with five biweeks of lag, with temperature (minimum, mean and maximum) at lag 0 ( Figure 5N), and one biweek of lag with rainfall kurtosis ( Figure 5P).    Parameter estimates for the best time series models for Ae. aegypti are presented in Table 1. Table 1. Parameter estimates for the best time series explaining abundance changes in Aedes aegypti sampled with BG-sentinel traps at Sanmin (dengue transmission hotsop) and Nanzih (dengue transmission coldspot), Kaohsiung City, Taiwan. In parameters, Intercept is the mean value of the time series and AR = autoregressive term, K = kurtosis, T = temperature, mean = average and min = minimum. Lags are in biweeks.

Parameters (lag)
Sanmin The biweekly average number (Intercept in Table 1) of Ae. aegypti mosquitoes was 50% higher in Sanmin (around 9) when compared with Nanzih (around 6) and in both locations the dynamics were sensitive to changes in minimum and mean temperature kurtosis, which were the more platykurtic environmental variables at KC during the study period ( Figure 2E). Parameter estimates for the best time series models for Cx. quinquefasciatus are presented in Table 2. Cx. quinquefasciatus average abundance at Sanmin (around 143) more than doubled what we found at Nanzih (around 63). In both districts increasing temperature had a negative impact on Cx. quinquefasciatus abundance, with a larger coefficient in Sanmin (Table 2), but only at Nanzih did the SD of Rainfall had a positive impact (Table 2) on Cx. quinquefasciatus abundance. Table 2.
Parameter estimates for the best time series explaining abundance changes in Culex quinquefasciatus sampled with BG-sentinel traps at Sanmin (dengue transmission hotsop) and Nanzih (dengue transmission coldspot), Kaohsiung City, Taiwan. In parameters Intercept is the mean value of the time series and SD = Standard deviation, T = temperature, max = maximum, R = Rainfall. Lags are in biweeks.

Discussion
Our data shows that, during the study period (March 2016-December 2017), a dengue transmission hotspot consistently had more mosquitoes than a neighboring transmission coldspot at KC, Taiwan. This fact supports the hypothesis that differences in transmission between these two districts might be related to differences in vector abundance. Moreover, the fact that both Ae. aegypti and Cx. quinquefasciatus were more abundant in the hotspot than the coldspot highlights the need to focus vector surveillance and control in hotspot areas. Specific to dengue transmission ecology in KC, our data showed the presence of Ae. aegypti near households, an observation furtherly supporting that Ae. aegypti might play the most important role in dengue virus transmission in southern Taiwan [79,80].
Our longitudinal mosquito data also highlights the importance of high order environmental variability, measured by kurtosis, on forcing mosquito abundance fluctuations, a phenomenon observed elsewhere for Ae. aegypti [32,39], Cx. quinquefasciatus [81][82][83] and other mosquito species [23,38,71,84], with patterns following the prediction of Schmalhausen's law, the biological principle stating that organisms are sensitive not only to average environmental conditions but also to the environmental variability per se [25]. In our study, both Ae. aegypti and Cx. quinquefasciatus were sensitive to temperature, the most platykurtic environmental variable during the study period, yet Ae. aegypti fluctuations were associated with changes in temperature kurtosis, while Cx. quinquefasciatus tracked changes in maximum temperature. These differences might be related to different impacts of temperature on the life history of these two mosquito species. For example, Ae. aegypti abundance changes might reflect an accelerated larval development to avoid death by extreme hot temperatures, and/or by larval habitat desiccation [32,39,85]. Meanwhile, abundance changes in Cx. quinquefasciatus could be related to decreased oviposition [86,87] and reduced aquatic survival at extreme high temperatures [88]. The different feeding behaviors of Ae. aegypti (day time biter) and Cx. quinquefasciatus (night time biter) might also be associated with the environmental conditions at their active time [10,89].
Although our longitudinal study was long enough to show the importance of platykurtic environmental factors on mosquito abundance fluctuations, it is an open question if impacts of the El Niño Southern Oscillation (ENSO) on dengue transmission [20,[90][91][92] could be related to changes in vector abundance. For example, in leishmaniasis we have shown that vector abundance and transmission increase during the cold ENSO phase in Panamá [93,94]. One limitation of this study is that given the lack of baseline longitudinal studies it is impossible to assess if ENSO impacted the population dynamics of both Ae. aegypti and Cx. quinquefasciatus, especially considering this study was done during one of the strongest ENSO events on record [66,67], which altered weather patterns at KC (Supplementary Figure S2 [Online Only]), and also altered weather and phenological patterns of mosquitoes elsewhere in east Asia, as documented for Tripteroides bambusa (Yamada) in Nagasaki, Japan, where the comparison of records from the 1970's and 1980's with data from 2016 and 2017 showed an increase in winter temperatures and that pupal recruitment season has lengthened its duration [95].
One potential limitation of this study was the use of weather data from KC's weather station, which was around 12 km apart from the study sites. Our previous experience has been that over spatial scales similar to the one of this study [68,95], and even larger [71] temperature fluctuations are more synchronous than mosquito abundance, and very similar to observations from weather stations located at a distance similar to the one between KC's weather station and our study locations [95]. Based on that, we believe data from KC's weather station was appropriate to evaluate the association between weather fluctuations and mosquito abundance at the two study sites. Nevertheless, it is an open question whether more finely grained weather data, from weather sensors attached to the traps, could help to improve the understanding of weather impacts on vector abundance at finer spatial scales.
Finally, our findings, in conjunction with observations about the high correlation between entomological indicators and dengue transmission [54], support the need to strengthen systematic vector surveillance at sites that consistently behave as dengue transmission hotspots and the importance of mosquito abundance for transmission, since we significantly found more mosquitoes in Sanmin (the dengue transmission hotspot) than in Nanzih (the dengue transmission coldspot), though these patterns might also be affected by other environmental, demographic, and social aspects important for vector-borne disease transmission [24,25]. In that sense, the development of intelligent mosquito traps that automatically count and identify mosquitoes is a very promising tool that could increase the reliability of early warning systems for vector-borne disease transmission, which primarily rely on weather data [96,97] given delays associated with the processing of entomological information [49,98].

Conclusions
Our data suggest that spatially heterogeneous Ae. aegypti abundance might be a driver of spatial differences in dengue transmission at KC. Our comparative approach, where we also observed a significant increment in the abundance of Cx. quinquefasciatus in the dengue transmission hotspot, suggests this area might be more likely to experience outbreaks of other vector borne diseases and should become a primary focus for vector ecology research, surveillance and control.