Landscape and Environmental Factors Influencing Stage Persistence and Abundance of the Bamboo Mosquito, Tripteroides bambusa (Diptera: Culicidae), across an Altitudinal Gradient

The bamboo mosquito, Tripteroides bambusa (Yamada) (Diptera: Culicidae), is a common insect across East Asia. Several studies have looked at the ecology of Tr. bambusa developmental stages separately, but little is known about the factors associated with the persistence (how often) and abundance (how many individuals) of Tr. bambusa stages simultaneously studied across a heterogeneous landscape. Here, we ask what environmental and landscape factors are associated with the persistence and abundance of Tr. bambusa stages across the altitudinal gradient of Mt. Konpira, Nagasaki City, Japan. During a season-long study we counted 8065 (7297 4th instar larvae, 670 pupae and 98 adults) Tr. bambusa mosquitoes. We found that persistence and abundance patterns were not associated among stages, with the exception of large (4th instar) and small (1st to 3rd instars) larvae persistence, which were positively correlated. We also found that relative humidity was associated with the persistence of Tr. bambusa aquatic stages, being positively associated with large and small larvae, but negatively with pupae. Similarly, landscape aspect changed from positive to negative the sign of its association with Tr. bambusa pupae and adults, highlighting that environmental associations change with life stage. Meanwhile, Tr. bambusa abundance patterns were negatively impacted by more variable microenvironments, as measured by the negative impacts of kurtosis and standard deviation (SD) of environmental variables, indicating Tr. bambusa thrives in stable environments, suggesting this mosquito species has a finely grained response to environmental changes.


Introduction
Mosquitoes (Diptera: Culicidae) are among the best-studied insects because of their role as pathogen vectors [1]. The significance of mosquito ecology is increasingly recognized as a key component to successfully manage mosquito nuisance and reduce disease transmission [2]. For example, mathematical modeling [3,4] and field observations [5,6] have robustly suggested that mosquito abundance is a crucial determinant of mosquito-borne disease transmission, a pattern also observed in vector-borne diseases For further details please refer to the inset legend. In the histogram legend, 'Larvae' indicates fourth instar larvae, and 'Small Larvae' first, second, and third instar larvae. The inset map shows the location of Kyushu island in western Japan and Nagasaki city, which is west of Osaka and Tokyo, the two most important cities in Japan.

Mosquito Sampling
Data presented in this study correspond to mosquitoes collected between May 2014 and June 2015. On 17 May 2014, we selected 27 sampling points, or focal trees, by following earlier mosquito research on Mt. Konpira [78]. On that day ovitraps, which were 350 mL Coca-Cola aluminum cans painted with black acrylic paint inside and outside, were placed on each focal tree [23,49]. We used ovitraps in this study given that previous studies have suggested mosquito species richness and abundance in ovitraps and treeholes of a similar volume are comparable at the study site [47,78]. For further details please refer to the inset legend. In the histogram legend, 'Larvae' indicates fourth instar larvae, and 'Small Larvae' first, second, and third instar larvae. The inset map shows the location of Kyushu island in western Japan and Nagasaki city, which is west of Osaka and Tokyo, the two most important cities in Japan.

Mosquito Sampling
Data presented in this study correspond to mosquitoes collected between May 2014 and June 2015. On 17 May 2014, we selected 27 sampling points, or focal trees, by following earlier mosquito research on Mt. Konpira [78]. On that day ovitraps, which were 350 mL Coca-Cola aluminum cans painted with black acrylic paint inside and outside, were placed on each focal tree [23,49]. We used ovitraps in this study given that previous studies have suggested mosquito species richness and abundance in ovitraps and treeholes of a similar volume are comparable at the study site [47,78]. Ovitraps were uniformly set at 1.2 m from the ground, following Zea Iriarte et al. [78] who found that Tr. bambusa at that height had abundance patterns similar to the ones observed in natural treeholes [47], and ovitraps were fixed with a black cord, which went through a 5 mm diameter hole that served as drainage when liquid contents went over 280 mL [49]. After each ovitrap was fixed it was filled with 280 mL of rain water [49]. Most of the sampling points were located within 5 meters from a mountain path mantained by Nagasaki City Mayor's Office direction of green and recreational areas, "Midori no ka." The altitude of sampling points ranged from 109 to 330 meters ( Figure 1). Ovitraps were set along three transects, and for reference, ovitraps are individually identified in Figure 2. Mosquitoes were collected every two weeks for a whole season. This sampling time was based on the developmental time of Tr. bambusa raised at temperatures fluctuating between 24 and 26 • C, which from egg to adult takes between 20 and 29 days, taking between 4 and 5 days for egg hatching, 13 to 19 days as larvae, and 4 to 5 as pupae before adult emergence [75]. Aquatic stages, i.e., pupae, large (4th instar) and small (1st, 2nd and 3rd instar) larvae, were sampled between 14 June 2014, and 24 June 2015 (28 sampling sessions), provided Tr. bambusa overwinters as larvae [77]. Fourteen sampling sessions, conducted between 18 May and 15 November 2014, covered the adult activity season [14]. Adult sampling was performed within a 2.5 m radius from focal trees with an entomological net (36 cm diam; Model 61-1B; Shiga Insect Co., Tokyo, Japan). Adult sampling at each location consisted of three steps: two-minute sweep, one-minute pause from sweeping, and a second two-minute sweep around the body of the person sampling. The second sweep collects mosquitoes fanned off during the first two minutes of net sweeping. Mosquitoes collected with the net were then sucked with an aspirator and dimethyl ether was put in the aspirator to kill the mosquitoes using a dust blower (Model AD-ECOM; ELECOM Co., Osaka, Japan). Samples were then placed in plastic centrifuge tubes (15 mL; Model ECK-15ML; AS-1 Co., Osaka, Japan), using cotton pads between layers of mosquito samples to minimize anatomical damage. Aquatic stage sampling was done by transferring all contents from each ovitrap to a 150 × 220 × 45 mm clear plastic pan (Mujirushi Ryōhin Co. LTD. Tokyo, Japan), where all pupae were removed, fourth instar Tr. bambusa larvae counted, and the presence of first to third instar larvae (hereafter referred to as "small larvae") checked as positive or negative. Tr. bambusa larvae field identification is facilitated by being the primary, often the only, mosquito species belonging to the Sabethini tribe in Nagasaki [77,78]. The identification of collected adults and pupae was done following the taxonomic key by Tanaka, Mizusawa, and Saugstad [65] using a dissection scope.
Ovitraps were uniformly set at 1.2 m from the ground, following Zea Iriarte et al. [78] who found that Tr. bambusa at that height had abundance patterns similar to the ones observed in natural treeholes [47], and ovitraps were fixed with a black cord, which went through a 5 mm diameter hole that served as drainage when liquid contents went over 280 mL [49]. After each ovitrap was fixed it was filled with 280 mL of rain water [49]. Most of the sampling points were located within 5 meters from a mountain path mantained by Nagasaki City Mayor's Office direction of green and recreational areas, "Midori no ka." The altitude of sampling points ranged from 109 to 330 meters ( Figure 1). Ovitraps were set along three transects, and for reference, ovitraps are individually identified in Fiure 2. Mosquitoes were collected every two weeks for a whole season. This sampling time was based on the developmental time of Tr. bambusa raised at temperatures fluctuating between 24 and 26 °C, which from egg to adult takes between 20 and 29 days, taking between 4 and 5 days for egg hatching, 13 to 19 days as larvae, and 4 to 5 as pupae before adult emergence [75]. Aquatic stages, i.e., pupae, large (4th instar) and small (1st, 2nd and 3rd instar) larvae, were sampled between 14 June 2014, and 24 June 2015 (28 sampling sessions), provided Tr. bambusa overwinters as larvae [77]. Fourteen sampling sessions, conducted between 18 May and 15 November 2014, covered the adult activity season [14]. Adult sampling was performed within a 2.5 m radius from focal trees with an entomological net (36 cm diam; Model 61-1B; Shiga Insect Co., Tokyo, Japan). Adult sampling at each location consisted of three steps: two-minute sweep, one-minute pause from sweeping, and a second two-minute sweep around the body of the person sampling. The second sweep collects mosquitoes fanned off during the first two minutes of net sweeping. Mosquitoes collected with the net were then sucked with an aspirator and dimethyl ether was put in the aspirator to kill the mosquitoes using a dust blower (Model AD-ECOM; ELECOM Co., Osaka, Japan). Samples were then placed in plastic centrifuge tubes (15 mL; Model ECK-15ML; AS-1 Co., Osaka, Japan), using cotton pads between layers of mosquito samples to minimize anatomical damage. Aquatic stage sampling was done by transferring all contents from each ovitrap to a 150 × 220 × 45 mm clear plastic pan (Mujirushi Ryōhin Co. LTD. Tokyo, Japan), where all pupae were removed, fourth instar Tr. bambusa larvae counted, and the presence of first to third instar larvae (hereafter referred to as "small larvae") checked as positive or negative. Tr. bambusa larvae field identification is facilitated by being the primary, often the only, mosquito species belonging to the Sabethini tribe in Nagasaki [77,78]. The identification of collected adults and pupae was done following the taxonomic key by Tanaka, Mizusawa, and Saugstad [65] using a dissection scope.

Environmental and Landscape Covariates
For the analysis, we used a digital elevation model and ground-and satellite-derived environmental covariates. At each sampling location, we measured canopy openness (Figure 3), i.e., how much light gets to the ground when passing through foliage, including its mean and standard deviation [49]. A ground cover index was estimated using the first component of a principal component analysis [23], the scale for which had negative values for grounds dominated by concrete and positive values for grounds dominated by leaf litter (Figure 3). We also measured air temperature and relative humidity at each sampling location [38], and used mean values from the season long sampling period for the subsequent spatial analyses. When sampling aquatic stages we also recorded ovitrap water temperature [49], of which the season long average values for each ovitrap were used in the statistical analysis. We did not record the water volume from each ovitrap, since this variable has mainly a qualitative effect on this species [77], and instead just recorded whether ovitraps had water or were dry.

Environmental and Landscape Covariates
For the analysis, we used a digital elevation model and ground-and satellite-derived environmental covariates. At each sampling location, we measured canopy openness (Figure 3), i.e., how much light gets to the ground when passing through foliage, including its mean and standard deviation [49]. A ground cover index was estimated using the first component of a principal component analysis [23], the scale for which had negative values for grounds dominated by concrete and positive values for grounds dominated by leaf litter (Figure 3). We also measured air temperature and relative humidity at each sampling location [38], and used mean values from the season long sampling period for the subsequent spatial analyses. When sampling aquatic stages we also recorded ovitrap water temperature [49], of which the season long average values for each ovitrap were used in the statistical analysis. We did not record the water volume from each ovitrap, since this variable has mainly a qualitative effect on this species [77], and instead just recorded whether ovitraps had water or were dry. The satellite data used for this study came from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) on board the NASA Earth Observing System Terra platform, and the Operational Land Imager (OLI) on the joint NASA/USGS Landsat 8 spacecraft. The ASTER instrument captures high spatial resolution surface data in 14 bands, from the visible to the thermal infrared wavelengths, with stereo viewing capability for digital elevation model creation [79]. We used highly resolved (15 m/pixel) ASTER retrieval data from bands 2 (Red) and 3 (Near Infrared, The satellite data used for this study came from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) on board the NASA Earth Observing System Terra platform, and the Operational Land Imager (OLI) on the joint NASA/USGS Landsat 8 spacecraft. The ASTER instrument captures high spatial resolution surface data in 14 bands, from the visible to the thermal infrared wavelengths, with stereo viewing capability for digital elevation model creation [79]. We used highly resolved (15 m/pixel) ASTER retrieval data from bands 2 (Red) and 3 (Near Infrared, NIR) for each sampled location. Landsat 8 OLI data (30 m/pixel) from bands 4 (Red) and 5 (NIR) [80] were Insects 2019, 10, 41 6 of 17 used to supplement the ASTER data. Images were first corrected using the dark object subtraction method [81]. Next, we estimated the normalized difference vegetation index (NDVI), an index of vegetation growth in which high values indicate abundant vegetation and low values the absence of vegetation [82], using the following equation: Given their different spatial resolution, NDVI was estimated separately for ASTER and Landsat 8 data, hereafter, referred to as NDVI-Landsat and NDVI-ASTER. Both Landsat 8 and ASTER data products were retrieved from the online data pool, courtesy of the NASA Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota [83]. The satellite images spanned several months, to ensure the lack of a seasonal bias in the estimates. Details about the images employed to estimate NDVI are presented in supplementary online Table S1. A digital elevation model from the Geospatial Information Authority of Japan [84] was used to estimate a series of landscape parameters that were included in our models for stage-specific persistence and abundance of Tr. bambusa. More specifically, we estimated (units inside parenthesis) elevation (m), slope (degrees), aspect (degrees), flow direction (power of two), roughness (m), and terrain roughness index (m) for each sampling location. Details regarding the estimation of these parameters and their interpretation can be found in Chaves and Moji [49].

Statistical Modeling
To study persistence patterns, we counted the total number of times individuals from the different stages (i.e., adults, pupae, and large and small larvae of Tr. bambusa) were present at each location during the sampling period. In contrast, Tr. bambusa abundance patterns were examined using the total number of individuals by stage (adults, pupae, and fourth instar larvae) at each location during the study period. Given the count nature of persistence and abundance [85], we fitted Poisson generalized linear models (P-GLM) and negative binomial generalized linear models (NB-GLM). For each stage, we initially used P-GLM, then, when the P-GLM model diagnostics suggested over-dispersion, we employed the NB-GLM. For both persistence and abundance, we fitted an initial set of 12 models that arose from the combination of variables described in the next lines. All the initial models considered the following covariates: elevation, canopy openness, ground cover index, and the mean, SD, and kurtosis of relative humidity. The mean, SD, and kurtosis of air and ovitrap water temperature were considered in models for adults and aquatic stages, respectively. To account for heterogeneities in vegetation growth at different spatial scales, the mean, SD, and kurtosis of NDVI, estimated with either Landsat 8 (30 m resolution) or ASTER (15 m resolution) data, were considered in each of the models. We considered both SD and kurtosis as both variables give different ideas of variability [86]: while SD measures the dispersal around the mean, larger values indicating an overall larger variability [87], kurtosis indicates if the observations are more variables near (platykurtic, the second and third quartiles are wider than the first and fourth quartiles) or far (leptokurtic, the first and fourth quartiles are wider than the second and third quartiles) from the median of the distribution [35]. The models also included one of the following variables: slope, roughness, or terrain roughness index. Either the aspect or flow direction was incorporated to consider the direction of elevation changes in the landscape. Models for aquatic stages included the natural logarithm of the number of days ovitraps had water as an offset variable, to account for the unequal number of times aquatic stages could be present at a sampling location as modulated by water availability [85].
Each one of the 12 initial models, for both persistence and abundance, was simplified by a process of backward elimination. In this process, covariates are removed one at a time, and the model is selected when an optimum value is found [88]. For model selection, we employed the Akaike Information Criterion (AIC), which weighs the trade-off between model goodness of fit and the number of parameters [88]. Among the models with the same number of parameters, the best model is that which minimizes the AIC. A list of the variables considered in the models is shown in supplementary online Table S2. Model selection for variables explaining Tr. bambusa persistence by stages are presented in: supplementary online Table S3 for adults, supplementary online Table S4 for  pupae, supplementary online Table S5 for fourth instar larvae, and supplementary online Table S6 for small larvae. Model selection for variables explaining Tr. bambusa abundance models by stage are presented in supplementary online Table S7 for adults, supplementary online Table S8 for pupae,  and supplementary online Table S9 for fourth instar larvae. Residuals from the best models were checked for model assumptions. The Moran's I index, a test for which the null hypothesis is spatial independence, was estimated for model residuals to check the assumption of spatial independence [89].

Software
Maps and geographical information systems (GIS) procedures were made using QGIS (version 2.18.10, Open Source Geospatial Foundation, Chicago, USA). All statistical analyses and other plots were made using R (version 3.4. 4, 64 bit, R Foundation for Statistical Computing, Vienna, Austria). We employed the following libraries: moments, for kurtosis estimation; MASS, for negative binomial model fitting; spdep, for Moran's I index; and base, for all other procedures and data plotting.

Results
Tr. bambusa stages had different persistence patterns. Figure 1 shows mosquito persistence patterns across the study site. Small and fourth instar larvae were present at the study site through all the study periods ( Figure 1). Across all sampling locations, the average persistence (±SD) for fourth instar larvae was 23.7 ± 4.0, and for small larvae was 22.7 ± 3.8. As shown in Figure 1, small larvae were always present at three sites, and fourth instar larvae at five sites. A factor likely modulating persistence of the aquatic stages was the proportion of time an ovitrap was dry. Figure 1 also shows that Tr. bambusa was present across all land cover use types. In contrast to the patterns observed for the larvae, pupae and adults were less persistent during the study period. The mean ± SD persistence for pupae at the study site was 7.8 ± 2.3, while for adults it was 1.7 ± 1.0. The maximum number of times pupae and adults persisted at a given location were 11 and 4 times, respectively. Figure 4 shows Tr. bambusa persistence associations by stage. No clear association was observed among the persistence patterns of adults and pupae ( Figure 4A), adults and fourth instar larvae ( Figure 4B), adults and small larvae ( Figure 4C), pupae and fourth instar larvae ( Figure 4D), or pupae and small larvae ( Figure 4E). Estimated Pearson correlations among those stages were not statistically significant (p > 0.05). In contrast, the Pearson correlation between fourth instar larvae and small larvae ( Figure 4F) persistence was positive (r = 0.68) and statistically significant (t = 4.691, df = 25, p < 0.05).
The statistical analysis of the association between Tr. bambusa mosquito persistence and environmental variables showed that adult persistence (Table 1) was negatively associated with landscape aspect. The direction of the slope (landscape aspect) indicated that adult mosquitoes were more common on the northern transect, where the slope was overall directed towards the northeast ( Figure 5A). Meanwhile, pupae persistence (Table 1) was positively and significantly (p < 0.05) associated with the aspect, meaning pupae were more persistent where the slope had a southwest orientation ( Figure 5A). Pupae persistence was also positively associated with ovitrap temperature (p > 0.05), indicating that pupae persistence was more common in the northern and western transects ( Figure 5B). In contrast, pupae persistence (Table 1) was negatively associated with both the mean (p > 0.05) and SD (p < 0.05) of relative humidity, which implies pupae were likely more common at the less humid sampling locations ( Figure 5C,D). The fourth instar larvae persistence (Table 1) was positively associated with mean relative humidity (p < 0.05), indicating this stage was most common at the sites with high relative humidity ( Figure 5C). Small larvae persistence (Table 1) was also positively associated (p < 0.05) with high relative humidity sites ( Figure 5C), but high vegetation growth (NDVI Adult persistence as function of (A) pupae persistence; (B) fourth instar larvae persistence; and (C) small larvae persistence. Pupae persistence as function of (D) fourth instar larvae persistence; (E) small larvae persistence. Fourth instar larvae persistence as function of (F) small larvae persistence. In all plots, numbers and colors represent, respectively, the sampling locations and transects described in Figure 2.  Adult persistence as function of (A) pupae persistence; (B) fourth instar larvae persistence; and (C) small larvae persistence. Pupae persistence as function of (D) fourth instar larvae persistence; (E) small larvae persistence. Fourth instar larvae persistence as function of (F) small larvae persistence. In all plots, numbers and colors represent, respectively, the sampling locations and transects described in Figure 2.    Adult abundance as function of (A) pupae abundance and (B) fourth instar larvae abundance. Pupae abundance as function of (C) fourth instar larvae abundance. In all plots, numbers and colors represent, respectively, the sampling locations and transects described in Figure 2.
The statistical analysis of Tr. bambusa stage-specific abundance patterns showed that adults ( Table 2) increased their abundance as a function of the average ( Figure 5B) and SD ( Figure 5C) of relative humidity. Adult abundance decreased with increasing variability, both SD ( Figure 7A) and kurtosis ( Figure 7B), in vegetation growth at a scale of 30 m (NDVI-Landsat), but also with canopy openness ( Figure 7C). Adult abundance also decreased with kurtosis relative humidity ( Figure 7D) and aspect ( Figure 5A). Meanwhile, Tr. bambusa pupae abundance (Table 3) was negatively associated    Adult abundance as function of (A) pupae abundance and (B) fourth instar larvae abundance. Pupae abundance as function of (C) fourth instar larvae abundance. In all plots, numbers and colors represent, respectively, the sampling locations and transects described in Figure 2.
The statistical analysis of Tr. bambusa stage-specific abundance patterns showed that adults ( Table 2) increased their abundance as a function of the average ( Figure 5B) and SD ( Figure 5C) of relative humidity. Adult abundance decreased with increasing variability, both SD ( Figure 7A) and kurtosis ( Figure 7B), in vegetation growth at a scale of 30 m (NDVI-Landsat), but also with canopy openness ( Figure 7C). Adult abundance also decreased with kurtosis relative humidity ( Figure 7D) and aspect ( Figure 5A). Meanwhile, Tr. bambusa pupae abundance (Table 3) was negatively associated Figure 6. Abundance association among Tripteroides bambusa stages at Mt Konpira, Nagasaki, Japan. Adult abundance as function of (A) pupae abundance and (B) fourth instar larvae abundance. Pupae abundance as function of (C) fourth instar larvae abundance. In all plots, numbers and colors represent, respectively, the sampling locations and transects described in Figure 2.
The statistical analysis of Tr. bambusa stage-specific abundance patterns showed that adults ( Table 2) increased their abundance as a function of the average ( Figure 5B) and SD ( Figure 5C) of relative humidity. Adult abundance decreased with increasing variability, both SD ( Figure 7A) and kurtosis ( Figure 7B), in vegetation growth at a scale of 30 m (NDVI-Landsat), but also with canopy openness ( Figure 7C). Adult abundance also decreased with kurtosis relative humidity ( Figure 7D) and aspect ( Figure 5A). Meanwhile, Tr. bambusa pupae abundance (Table 3) was negatively associated with the kurtosis of NDVI-Landsat ( Figure 7B), mean ( Figure 5B) and SD ( Figure 5C) of relative humidity, and the kurtosis of ovitrap water temperature ( Figure 7E). On the other hand, pupae abundance was positively correlated with aspect ( Figure 5A) and roughness ( Figure 7F). The abundance of Tr. bambusa fourth instar larvae (Table 4) was negatively associated with the NDVI-Landsat SD ( Figure 7A) and kurtosis ( Figure 7B), elevation ( Figure 7G), relative humidity kurtosis ( Figure 7D), ovitrap water temperature mean ( Figure 5D) and SD ( Figure 7H), landscape roughness ( Figure 7F), and aspect ( Figure 5A). with the kurtosis of NDVI-Landsat ( Figure 7B), mean ( Figure 5B) and SD ( Figure 5C) of relative humidity, and the kurtosis of ovitrap water temperature ( Figure 7E). On the other hand, pupae abundance was positively correlated with aspect ( Figure 5A) and roughness ( Figure 7F). The abundance of Tr. bambusa fourth instar larvae (Table 4) was negatively associated with the NDVI-Landsat SD ( Figure 7A) and kurtosis ( Figure 7B), elevation ( Figure 7G), relative humidity kurtosis ( Figure 7D), ovitrap water temperature mean ( Figure 5D) and SD ( Figure 7H), landscape roughness ( Figure 7F), and aspect ( Figure 5A).  The models presented in Tables 1-4 had Moran's I indices supporting spatial independence. The models for pupae (Table 3) and fourth instar larvae (Table 4) abundance had statistically significant (p < 0.05) overdispersion parameters, indicating that using a negative binomial model was an appropriate statistical modeling choice.   The models presented in Tables 1-4 had Moran's I indices supporting spatial independence. The models for pupae (Table 3) and fourth instar larvae (Table 4) abundance had statistically significant (p < 0.05) overdispersion parameters, indicating that using a negative binomial model was an appropriate statistical modeling choice.

Discussion
Understanding the factors shaping population abundance and persistence are a major goal of ecology [10,[90][91][92]. In the applied context of insects with medical and economic importance, the need to understand these factors becomes even more relevant, since insect population management can help to reduce disease transmission [2,8], reduce nuisance [93] or increase food production [94]. Concerning the bamboo mosquito, our results show that factors associated with its persistence and abundance can change through its different life stages, in the case of persistence with some covariates changing the sign of their association through consecutive life stages. This result is interesting as it highlights that for insects with complex life cycles [95] the information provided by one life stage is not necessarily useful to predict spatial abundance patterns across all its life stages. Surprisingly, the spatial persistence of larvae and pupae were not correlated, despite the limited dispersal ability of mosquito aquatic stages in treeholes and the relative homogeneity of ovitraps as larval habitats. The lack of correlation in the spatial persistence of larvae and pupae indicates that heterogeneities in mosquito emergence are not solely related to habitat quality, as recognized for mosquitoes with medical importance [96,97]. This complex relationship was shown by the association between landscape aspect and Tr. bambusa persistence, where an increase in aspect had a negative impact on adult persistence. However, the same factor, landscape aspect, had a positive impact on pupae persistence. Similar patterns were also observed for Tr. bambusa stage abundance, while aspect had a negative impact on adult, which have a greater dispersal ability by their ability to fly [1], and fourth instar larvae abundance, the impact was positive for pupae, something that might be related to the density-dependent effects on persistence [14,24]. This pattern of variable association by stage was not unique to landscape variables, a similar pattern of changing associations between Tr. bambusa life stage and environmental variables was observed for relative humidity, which had a positive impact on larvae but negative on pupae persistence. Moreover, relative humidity had a positive impact on adult, but negative on pupae abundance.
A second pattern observed in Tr. bambusa spatial persistence and abundance was that different life stages were associated with different covariates. For example, small larvae were the only stage whose persistence was associated with NDVI from ASTER, while canopy openness was negatively associated with adult abundance, roughness positively associated with pupae and elevation negatively with fourth instar larvae. These unique stage-specific patterns of association have important implications for modeling the ecological niche of organisms with complex life cycles [98]. On the one hand they illustrate that ecological niche predictions based on specific stages of insects with complex life cycles are prone to be biased by ignoring information about all the life cycle stages, not to mention other caveats related to species interactions [38] and evolutionary changes, related to the response to changing environmental variables, over short periods of time [99]. On the other hand, inferences about ecological niches based on all life stages might improve predictions by ecological niche models as they are commonly implemented, for example, as has been done for medically important mosquito species [100].
In addition, Tr. bambusa persistence and abundance patterns were not only related to average values of environmental variables, but also to their patterns of variability, something to be considered when modeling the ecological niche of organisms, which tend to be based on mean environmental variables [98]. The relevance of environmental variability for Tr. bambusa and other organisms' persistence and abundance patterns is predicted by Schmalhausen's law. Both persistence and abundance patterns are predicted by Schmalhausen's law, the biological principle that states organisms are sensitive not only to average environmental conditions but, also, to their associated patterns of variability [54,101]. In this study, we considered patterns of environmental variability by measuring the SD and kurtosis of ecological variables. We found that overall, an increasing environmental variability had negative impacts on both Tr. bambusa persistence and abundance. This pattern is the opposite of what we have previously observed for medically important species co-occurring with Tr. bambusa at Mt. Konpira, where the globally invasive species Aedes albopictus and Aedes japonicus (Theobald) and the increasingly common Aedes flavopictus Yamada tended to be positively associated with increasing levels of environmental variability [24,38,49]. In contrast to what has been observed for Ae. japonicus [38,49], where increasing environmental variability was positively associated with its abundance, Tr. bambusa tended to be negatively associated with increasing environmental variability. This result might be the key to understanding, why despite having a propagule pressure similar to that of major globally invasive mosquito species [71,72], Tr. bambusa has not become widely distributed.

Conclusions
Our data showed that Tr. bambusa thrives under stable environmental conditions, suggesting the species has a fine environmental grain in response to changing environments [102]. This result suggests that evolutionary changes in response to climate change in this species could be minimal if its fitness set is convex or have directional changes if its fitness set is concave [102]. Then, these two scenarios, respectively, imply an expectation to see no genetic changes (convex fitness set) in this species or a reduction of local genetic diversity (concave fitness set) and the emergence of a genetic structure along its geographic distribution. Finally, given the ease to study Tr. bambusa ecology in the field, this species is an excellent model to study genetic changes in response to climate change and ecological traits used to predict the invasive potential of mosquitoes given its close association with globally invasive mosquito species over their native range.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2075-4450/10/2/41/s1, Table S1: Codes of Landsat 8 and Aster images employed to estimate NDVI at the sampling locations, Table S2: Abbreviations for the variables considered in the models; Table S3: Model selection for the best model explaining the persistence of Tripteroides bambusa adults at Mt Konpira, Nagasaki, Japan. The best model is bolded; Table  S4: Model selection for the best model explaining the persistence of Tripteroides bambusa pupae at Mt Konpira, Nagasaki, Japan. The best model is bolded; Table S5: Model selection for the best model explaining the persistence of Tripteroides bambusa fourth instar larvae at Mt Konpira, Nagasaki, Japan. The best model is bolded; Table S6: Model selection for the best model explaining the persistence of Tripteroides bambusa small larvae at Mt Konpira, Nagasaki, Japan. The best model is bolded; Table S7: Model selection for the best model explaining Tripteroides bambusa adult abundance at Mt Konpira, Nagasaki, Japan. The best model is bolded; Table S8: Model selection for the best model explaining Tripteroides bambusa pupae abundance at Mt Konpira, Nagasaki, Japan. The best model is bolded; Table S9: Model selection for the best model explaining Tripteroides bambusa fourth instar larvae abundance at Mt Konpira, Nagasaki, Japan. The best model is bolded.