### Analysis of the temporal variation in the larval habitat locations

The temporal variations in mosquito larval habitat locations were analyzed at two levels: (1) the recurrence of larval habitats in the same locations, and (2) occurrence of habitats at nearby locations. The analysis at these two levels can help answer two key questions: (1) how likely it is to find habitats at locations where habitats were previously found, and (2) to what are extent are new habitats spatially associated with habitats that have been found previously (where overlapping habitats are excluded)? These two questions are critical to mosquito habitat management efforts. If larval positive habitats tend to occur in the same or very proximal locations between seasons/years, resource management can be optimized much more easily. For the first level analysis, each of the six point maps was converted into a raster map, in which each habitat location is represented by a pixel and overlapping analysis was selected to determine the changes in habitat pixels between time points. Based on the field observations, the diameter of habitats ranged from 4 m to 15 m. To take into account effects of the varying size of habitats, the pixel sizes range from 4 m to 20 m (allowing for a 5 m GPS error), with a 2 m interval (half of the smallest diameter of habitats). The smallest and largest pixel sizes were chosen based on the observed diameters of the smallest and largest habitats in the field. Even though the pixel sizes were carefully selected, this overlapping approach still has a limitation: the overlapping area is expected to increase with increasing pixel size. To obtain a realistic estimation of the overlapping area, we plotted the percentages of overlapping aquatic habitats as a function of the size of pixels. For the rest of the analysis only the data on anopheline-positive habitats were used, since some stagnant aquatic habitats might not be suitable for mosquito reproduction and our focus is on the habitats that are productive.

For the second level analysis, a Nearest Neighborhood analysis was used to assess the distance between a habitat and its first and second nearest habitat locations. To determine the climatic impact, three types of comparison were performed: (1) the spatial overlapping and nearest neighbor of habitats between the February and the May within the same year; (2) the spatial overlapping and nearest neighbor of habitats in the February between different years and (3) the spatial overlapping and nearest neighbor of habitats in the May between different years. Even though the six datasets were collected separately in six different field surveys, the consistency in the sampling was maintained through the use of the same research group with the same training and sampling methods in the same study area.

### Analysis of the temporal variation in the geographic extent of larval habitats

The spatial patterns of habitats were quantified using two statistics: the dispersion pattern analysis and the compact analysis. The dispersion pattern analysis measures how dispersed habitat locations are around the geographic center of the study area. It takes into account dominant landforms (e.g., a river or cliff) and reveals global trends in a point pattern, such as orientation. The compactness analysis complements the dispersion pattern analysis by quantifying the shape of the range of mosquito larval habitats. The dispersion and compactness of habitat distribution have a direct impact on habitat management: the more dispersed and compact the habitats are in an area, the more management efforts are required for these habitats.

To examine the dispersion pattern of anopheline-positive habitats, the six sets of point data were analyzed using a Standard Deviation Ellipse (SDE) test. SDE calculates an ellipse to describe the dispersion of points, and has been widely used in the geographic analysis of point patterns [21, 22]. In this study, four statistics were derived from the ellipse: the length of the major axis (the maximum dispersion distance), the length of the minor axis (the minimum dispersion distance), the product of the lengths of the major and the minor axes (an approximation of the area of an ellipse), and the orientation of the dispersion. The orientation of an ellipse is measured as the degree by which the major axis is rotated from the geographic east in the counterclockwise direction [23]. The Standard Deviation Ellipse (SDE) was computed using CrimeStat 3.0 [23].

Because habitat geographies often have irregular shapes that can have important implications for habitat management [24], the shape of mosquito larval habitat ranges were explicitly analyzed using convex polygon analysis [25] combined with the Boyce-Clark index [26]. These two methods are widely adopted measures for habitat compactness and shape [24]. The convex polygon analysis identifies the smallest polygon that can enclose all habitat locations. For each set of habitat locations, a convex polygon was calculated and a total of six convex polygons was obtained. The Boyce-Clark index expresses how closely the shape of the habitat distribution (represented by the convex polygon) approximates that of a hypothetical circle. It is calculated by drawing a set of equally spaced radials from a location within a shape to its perimeter and measure. The variation in lengths of the radials is then calculated. In this study, 13 radials with 30 degree space were used (30 degree is selected based on the smallest angle between any two of the neighboring vertexes in each of the six convex polygons). This index ranges from 0 to 200. A value of 0 indicates the shape is close to a perfect circle. A value of 200 indicates the shape is close to a line.

After SDE analysis, convex polygon analysis and Boyce-Clark index were applied to each of the six sets of point data, pair-wise comparisons were carried out to determine the seasonal or inter-annual changes in larval distribution patterns. The three types of comparisons outlined in the first subsection of the statistical analysis were also performed: (1) the difference between the February and the May within the same year (three pairs of datasets); (2) the difference in the February between different years (three pairs of datasets); and (3) the difference in the May between different years (three pairs of datasets).

Analysis of the temporal variation in the niche patterns of larval habitats. Habitat management is often concerned with the ecological niche of species, which can be described by its mean position and breadth along various environmental gradients (often referred to as environmental axes) [

27]. Characterizing the ecological niche of mosquitoes is also critical for mosquito habitat management, as environmental gradients are routinely used as indicators for the distribution of mosquito larval habitats. Our previous study has determined that four environmental gradients, curvature, elevation, distance to streams and wetness index, are related to the mosquito habitats [

12]. This study quantified the niche position and niche breadth along these four environmental gradients using the marginality (M) (Equation 1) and the specialization (S) (Equation 2) measurements [

28]. The marginality (M) is a description of the mean niche position on each selected environmental gradient. As shown in Equation

1, the marginality is calculated by comparing the distance between the mean conditions used by the species and the mean conditions of the study area for that gradient.

where m_{g} is the mean value of an environmental gradient in the study area, m_{e} is the mean value of the environmental gradient in the locations occupied by habitats, and σ_{g} is the standard deviation of the environmental gradient in the study area.

In Equation

1, division by σ

_{g} is needed to remove any bias introduced by the variance of this gradient in the study area. The weight (1.96) ensures that the marginality will tend to be normalized between zero and one. This means if the distribution of a gradient in a study area is normal, 95% of the values lie within 1.96 standard deviations of the mean. A habitat distribution was then considered 'marginal' when its mean condition was far from the mean condition of the study area (its marginality value is close to one). As shown in Equation

2, the specialization (S) measure estimates the variability of habitat conditions by comparing the variance of the habitats on the selected gradient (named niche breadth) with the variance of the study area on this gradient.

where σ_{g} is the standard deviation of an environmental gradient in the study area and σ_{e} is the standard deviation of the environmental gradient in the areas occupied by habitats.

A habitat distribution was considered specialized if it has a narrow niche breadth (its specialization value is larger than one). If a species differs considerably in its marginality and specializations along an environmental gradient, this gradient is less likely to be adequately predictive for the distribution of this species [29, 30].