Analysis of the Spatial Distribution of Aedes Albopictus by Indicator Kriging in an Urban Area of Shanghai, China

Background: Aedes albopictus is a well-recognized vector of major arboviral diseases and a primary pest in tropical and temperate regions of China. The current monitoring system, based on sub-district scale, in most cities of China for the spread of Ae. albopictus has not taken spatial distribution for the analysis of density of species. So it is not accurate enough for epidemic investigations, especially in big cities like Shanghai. Methods: In this study, a new surveillance program integrate the actual monitoring locations was used to investigate the temporal and spatial distribution of Ae. albopictus abundance in an urban area of Shanghai, China, from 2018 to 2019 by using the mosquito-oviposition trap (MOT) method. The study area of 14 sub-districts was divided into 133 grids with a side length of approximately 500 m. The vector abundance and spatial structure of Ae. Albopictus were predicted by the indicator Kriging based on eight MOTs in each grid. Meanwhile, the light trap (LT) method was also used for the analysis and compared with the MOT method. Results: A total of 8,192 MOTs were placed in the study area in 2018, and 7,917 (96.6%) were retrieved with a positive rate of 6.45%, while in 2019, 22,715 (97.0%) of 23,408 MOTs were recovered with a positive rate of 5.44%. When using the LT method, 273 (93.5%) and 312 (94.5%) adult female Ae. albopictus were gathered in 2018 and 2019, respectively. The Ae. albopictus populations in the urban area of Shanghai increased slowly from May, reached a peak in July, and declined gradually from September. The MOT positivity index (MPI) showed a signicant positive spatial autocorrelation across the study area, while LT collections indicated a non-signicant spatial autocorrelation. The MPI was suitable for spatial interpolation by using the indicator Kriging, and showed different hotspots in different years. Conclusions: The new surveillance system integrate geographic information can help improve our understanding of the spatial and temporal distribution of Ae. albopictus in urban areas of Shanghai and could provide a practical method for decision-makers to


Background
Aedes albopictus (Diptera: Culicidae), also known as the Asian Tiger Mosquito, has invaded all continents except Antarctica during the last 30-40 years [1,2]. It is an aggressive vector of major arboviral diseases such as dengue, chikungunya, yellow fever, and Zika. In China, Ae. albopictus is a primary nuisance pest and disease vector in the tropical and temperate regions after having adapted to low temperatures [3]. Ae. albopictus is present in regions where Aedes aegypti is absent, including Shanghai [4]. Dengue is one of the most widely transmitted diseases by Ae. albopictus in China [5]. Ae. albopictus was reported as the primary vector of several epidemics in Guangzhou Province (37,354 laboratory-con rmed cases) in 2014 [6] and Zhejiang Province (adjacent to Shanghai) in 2004 [7], 2009 [8] and 2017 [9].
The dengue case reported in 2017 in Shanghai was the rst autochthonous dengue case in the last ve decades in Shanghai [10]. In the following years, three cases have been reported [11]. Since then, Ae. albopictus has been on the top of list for vector control and surveillance in Shanghai. The Public Health Service developed a monitoring system for Ae. albopictus to obtain information regarding its temporal evolution by using the light trap (LT) and mosquito-oviposition trap (MOT) methods in 2010. This system is used for surveillance of Ae. albopictus population and biting rates based on data collected from subdistricts. However, a disadvantage of the surveillance network is that it does not integrate geographic information and the scale is also large.
Geostatistical methods, which integrate the actual locations of samples, have been used to investigate the spatial distribution of mosquitoes [12] and several mosquito-transmitted diseases, including malaria [13,14] and Dengue fever [15,16]. Global and local indicators of spatial association such as Moran's I [17] or LISA (local indicators of spatial association) [18] have been applied to study pests, including mosquitoes [19,20]. These indicators can detect hot spots of mosquito abundance and predict the signi cance of clustering and the effect of disease control [21]. Among the geostatistical methods, Kriging interpolation [22] can predict the vector abundance in unsampled areas. Albieri et al. have used the Kriging interpolation to predict mosquito population distribution at the provincial and municipal scales in northern Italy [23]. Azil et al. have used Kriging to analyze the costs of dengue vector surveillance and control programs in Australia [24]. In addition, Giordano et al. conducted a more e cient larvicide control program for West Nile virus awareness campaigns in Canada by using the Kriging interpolation technique [25].
In China, assessment of the density of Aedes albopictus is usually made based on the positive rate in ovitrap monitoring, and MOT is a standard method for surveillance of the temporal and spatial distributions of container-inhabiting mosquitoes [26]. However, there have been rare reports on the spatial interpolation of positive rates and the indicator Kriging.
The present study was set out to evaluate the temporal and spatial distributions of Ae. albopictus using the MOT method, investigate the autocorrelation of Ae. albopictus abundance and estimate Ae. albopictus abundance at non-sampling locations using the indicator Kriging [27] in an urban area of Shanghai, and identify hotspots and risk areas of high infestation. We conducted a new surveillance program for Ae. albopictus in an urban area of Shanghai from 2018 to 2019. The 14 sub-districts were divided into 133 grids with a side length of approximately 500 meters. Eight MOTs in each grid were applied to predict Ae. albopictus abundance and spatial structure. Meanwhile, the LT method was also used for the analysis and compared with the MOT method.

Selection of the study area
Shanghai is situated at 31°12′N north latitude and 121°30′E east longitude in the east part of the alluvial plain of the Yangtze Delta, adjacent to the Yangtze River estuarine and the East China Sea. It has four distinct seasons and abundant precipitation with a subtropical monsoon climate. The mean annual temperature in Shanghai is approximately 17 ℃, and the mean annual precipitation is more than 1100 mm, with 53% occurring between June and September. The study area is situated in the center of Shanghai, China, with a total area of 37.37 km 2 , measuring 6.15 km from east to west and 11.93 km from south to north (Fig. 1). In 2018, the study region included 14 sub-districts with a resident population of 1,057,700 [28].

Meteorological data
The monthly total precipitation and monthly mean maximum and minimum temperatures were calculated based on data from the China Meteorological Administration [29]. Weather variables were recorded at Xujiahui, which is located approximately 2 km from the study area.

Entomological survey
The abundance of Ae. albopictus was analyzed by using the MOT (Tian®, Kaiqi Co. Ltd, Shanghai, China) and LT methods. The MOT trap [30] consists of a transparent cylindrical plastic jar (100 mm high, 70 mm diameter, 66 mm internal diameter) with a concave bottom (20 mm inward) and a black top cover with three conical openings of 100 mm diameter. When used as a collection container, a white lter paper of 70 mm diameter, used as an egg deposition substrate, was placed inside the bottom of the MOT, and 25 ml of dechlorination water was poured into the jar to keep the paper moist but not submerged. MOTs were placed outdoors on grasslands, kept away from direct sunlight, rain, and wind at ground level by a skilled technician, and maintained unchanged until the end of the study.
The MOTs were placed once a month in 2018 between April and November. To enhance the data in peak period, in 2019, the frequency was increased to once a week in week 20, week 23, week 25, weeks 27-39 and weeks [41][42][43][44][45][46]. MOTs that were removed, emptied, or interfered for any reasons were excluded from further analysis. After four days, each MOT was collected to the laboratory. Species identi cation was performed with a stereomicroscope, and the MOT positivity index (MPI) was calculated as: MPI = number of the Aedes-positive MOTs / total number of MOTs retrieved × 100%. In 2018, the center of the study area was under construction and could not be placed with MOTs, so we chose the rest area and divided it into 128 grids with a side length of approximately 500 m as surveillance units. We found a residential area with vegetative coverage as the monitoring point in each unit and put eight MOTs around a center point (Fig. 1). The MPIs of the eight MOTs were used to represent the Ae. albopictus density at that point. In 2019, we added ve surveillance units in the middle of the study area, which lead to a total of 133 surveillance units in 2019.
The LTs baited with carbon dioxide. On the third Wednesday of each week from May to November, the Centers for Disease Control and Prevention set two LTs in every sub-district throughout Shanghai as part of a city-wide mosquito surveillance program. LTs were usually collected from 4 pm to 10 pm. Contents of LTs were sent to the laboratory for species identi cation, and only female mosquitos from the traps were collected for data analyses.
The locations of MOTs and LTs were georeferenced using GPS CHCNAV X360H (WGS 1984 coordinate system) and later projected to the Shanghai local coordinate system. The georeferenced positions of MOTs and LTs in 2018 and 2019 are presented in Fig. 1.

Cluster analysis
The monitoring data were assigned into MOT2018, MOT2019, LT2018 and LT2019 groups and analyzed. The total collections of each LT and the mean MPI of each unit were calculated. Geostatistical analyses were conducted using ArcGIS 10.3 Spatial Statistics Tools (ESRI, China), and data in Microsoft Excel 2019 were imported into ArcGIS. The area of the minimum enclosing rectangle around the input features was used for cluster analysis based on Euclidean distance.
The average nearest neighbor (the mean distance from each feature to its nearest neighboring feature) analysis was then performed to calculate the nearest neighbor index. The index compared the mean distance to the null hypothesis states in which traps were randomly distributed, with index < 1 indicating a trend of clustering, while index > 1 indicating dispersion or competition.
We evaluated whether the mosquito abundances were spatially autocorrelated by calculating the incremental spatial autocorrelation (Global Moran's I at multiple distances). The Global Moran's I was tested using the permutation procedure based on feature locations and attribute values (total collections of each LT and mean MPI of each unit) against the null hypothesis (the absence of spatial autocorrelation). This analysis identi es the spatial patterns in the study area but did not indicate where such clusters occur, which could be determined by Local Moran's I. The local Moran's I analysis was performed to assess the presence of hot spots with statistically signi cant clusters, cold spots, and spatial outliers.

Geostatistical analysis
The indicator Kriging [27] was used for geostatistical analysis. The Kriging and co-kriging techniques were performed using the ArcGIS 10.3 Geostatistical Analyst extension, including exploratory statistical analysis, variogram modeling, and producing a probability or standard error of indicator. The Kriging interpolation method was used to quantify the spatial structure of the data and predict species abundance at unsampled locations on the basis of the mean MPI of each unit. MPI was made into a binary (0 or 1) variable, where 1 means MPI value above the threshold, and 0 means MPI value below the threshold. The resulting interpolation map shows the probabilities of exceeding the threshold. The co-kriging uses information on several variables. We used several thresholds for MPI and then the binary data on each threshold (covariable) to predict the threshold of primary interest by co-kriging. The prediction surface with the lowest error output and standard error surface was applied and subtracted from the study area. A leave-one-out cross-validation method was used to determine whether the Kriging interpolation provided reliable estimates of the indicator at unsampled locations. The criteria used for accurate prediction in the cross-validation were requested to be the following: Root mean square standardized (RMSS) approximately 1, mean standardized (MS) approximately 0, and root mean square (RMS) approximately the average standard error (ASE).

Mosquito collection
A total of 8,192 MOTs were placed in the study area, and 7,917 of them (96.6%) were retrieved with a positive rate of 6.

Spatial distribution of MOT and LT traps
We used the average nearest neighbor analysis to estimate whether the spatial autocorrelation of Ae. albopictus abundance was in uenced by the locations of MOTs and LTs. The average nearest neighbor analysis tested whether the traps were randomly distributed within the study area, and produced the mean, minimum, and maximum distances between traps ( Table 1). The results showed that the MOTs exhibited a dispersed distribution pattern in 2018 and 2019, while the LTs showed a random distribution pattern in 2018 and 2019.
The random or dispersed distribution pattern of MOTs and LTs indicated that the signi cant clustering of Ae. albopictus abundance was not due to the location of the traps in the study area.

Cluster analysis of Aedes albopictus abundance
We performed the incremental spatial autocorrelation analysis to evaluate the spatial autocorrelation of Ae. albopictus abundance across the study area. Only the MOT method demonstrated statistically signi cant positive spatial autocorrelation based on the mean MPI of each unit, which peaked at 609 m in 2018, and 542 m in 2019. The LT method did not show statistically signi cant spatial autocorrelation of Ae. albopictus abundance based on the total collections of each LT (Table 2).
Local Moran's I was then used to determine the locations of hot spots, cold spots, and spatial outliers. Results showed 13, 14, 2, and 1 location of clusters or outliers for MOT2018, MOT2019, LT2018, and LT2019, respectively (Fig. 6).

Prediction of Aedes albopictus abundance at non-sampling locations
We found that MPI of each unit was suitable for Kriging interpolation, and probability maps were created for the abundance.
There was no spatial autocorrelation for LT2018 and LT2019, and spatial interpolation was not permitted. We constructed the indicator Kriging interpolation by transferring the mean MPI of each unit into a binary with a threshold value of 5. Semivariograms created with indicator Kriging showed spatial dependence (range) within approximately 2000 m and 900 m for MOT2018 and MOT2019 (Figs. 7 and 8), respectively, beyond which the semivariance remained constant. The best-tting model for 2018 was the stable model with a nugget effect value (the semivariance value at zero distance) of 0.096, and a partial sill value (the constant semivariance value beyond the range) of 0.160, and the spherical model t the data of 2019 best, with a nugget effect value of 0, and a partial sill value of 0.259.
Prediction maps (Fig. 9) associated with those of standard errors (Fig. 10) based on MOT2018 and MOT2019 data showed that the highest mosquito abundance and strong spatial clustering were located in the south and north regions of the study area in 2018, and in the south and central areas in 2019. The prediction of standard errors quanti ed the degree of uncertainty for each location on the surface. According to this analysis, the prediction error was the lowest around where MOTs were set in the study area. Overall, the leave-one-out cross-validation statistics (Table 3) with the value of RMSS approaching 1 showed that the predicted models were reliable for map production.
Furthermore, Co-kriging improved the prediction of indicator Kriging. By using the Co-kriging method, inclusion of multiple thresholds of MPI generated several indicator variables for the same dataset. The use of covariables in the co-kriging analysis resulted in lower root mean square errors (RMS) for prediction (Table 3)

Discussion
According to the "Implementation Program of National Vector Surveillance" launched by China's Center for Disease Control and Prevention, MOT and LT should be routinely applied for Ae. albopictus monitoring in China. In this project, each subdistrict as a monitoring unit in Shanghai has a sample size of 2 LTs and 50 MOTs in a park and a residential area, but the spatial scale is limited for evaluating Ae. albopictus distribution [31,32]. To obtain more accurate seasonal and spatial distribution density of species, we developed this scheme without additional fund demand by dividing the sub-district into grids. We also evaluated the feasibility and accuracy of MOTs in combination with geostatistical analysis as a practical tool for monitoring temporal and spatial distribution of Ae. Albopictus in an urban area of Shanghai.
The results showed that MPI peaked in July in both 2018 and 2019, while the LT collections peaked in July in 2018 and in August in 2019. However, the indices of LT remained at a high level from July to September in both years. Consistent with the study by Gao et al., there was a signi cant correlation between monthly sampling yields [33]. The Ae. albopictus populations in the urban area of Shanghai slowly increased from May, peaked from July to September, and declined after September, which is coincident with the seasonal high temperatures and precipitation, and also consistent with previous reports [34,35].
Moreover, although there was little difference in the quantity of Ae. albopictus captured by the LT method between 2018 and 2019, the MPI in 2018 was higher than that in 2019. This can be explained by more precipitations in the 2019 mosquito season, which bought more water-lled habitats for mosquito reproduction. As Ae. albopictus has a strong oviposition preference for physical and chemical stimuli over a range of distances [36], the increased number of water-lled containers would distract gravid Ae. albopictus females from laying eggs in the MOTs [37,38]. Therefore, our results indicated that both MOT and LT could be applied for surveillance of the seasonal uctuation of Ae. albopictus, but MPI and LT might be different between years.
In this study, we did not observe spatial autocorrelation and found less Local Moran's I clusters for the LT collections of different periods, possibly due to the small number or the low density of LTs during those periods. Since the spatial analyses used distances to establish neighbors, low amounts of neighbors might result in lower statistical signi cance. Adult Ae. albopictus were dispersed at a maximum distance of 600 to 800 m [39,40]. Duncombe et al [41] suggested that the mosquito traps be placed at a distance of less than 1200 m from each other, while in this study the maximum distance to the nearest neighbor was over 1500 m ( Table 1). The spacing of LTs was too large to detect spatial autocorrelation in the sample area [42]. In our study, there were 2 LTs arranged in each sub-district with different areas, which contributed to the different spatial density of traps. More LTs are needed and should be well distributed for Ae. albopictus surveillance in Shanghai.
In this study, there was signi cantly positive spatial autocorrelation of MOT2018 and MOT2019 across the study area with a maximum distance of peak Global Moran's I around 600 m ( Table 2). The semivariograms showed that mosquito collections from MOT more than 2000 m in 2018 and 900 m in 2019, did not show spatial autocorrelation (Figs. 7 and 8). The MOT collections were suitable for spatial interpolation using the indicator Kriging technique. We chose MPI of 5 as the threshold, which was also used as the threshold for Ae. albopictus sampling in Shanghai, China [33]. The indicator Kriging maps revealed different oviposition hot spots in 2018 versus 2019, which may be attributed to the following reasons: rst, the mosquito abundance and seasonal distribution vary from year to year due to changes in temperature, precipitation, and humidity [43,44]; second, the result of MOT2018 may provide a reference for strict vector control in high-density areas in 2019, causing a lowered MPI density in these sites compared to other regions in 2019. The Co-Kriging is recommended to improve the accuracy of spatial interpolation and the prediction of Ae. albopictus abundance [45]. The cross-validation results indicated that the estimated MPI at unsampled locations was reasonably acceptable, despite some limits due to the uneven distribution of the data. Interpolation error of the Kriging method occurs most frequently on borders of study areas due to the "edge effects" [46] associated with the fewer neighbors. The new Ae. albopictus monitoring system is a reliable tool to determine the spatial variation within MPI data at the grid-scale.
As far as we know, this is the rst study to apply MOT together with geostatistical methods to develop a routine surveillance program in China. We found that MOT and LT could be used to monitor the relative size of mosquito populations [47] and that there was a signi cant correlation between monthly MPI and LT in our study. The use of LT collections can help set a threshold based on the correlations between LT and MOT. Besides, the indicator Kriging method is more suitable for MOT data with the advantage of applying binominal variables and accommodating zero values of the data. By using the indicator Kriging and an appropriate threshold, the possibility that MOT considerably underestimates Ae. albopictus population in areas of high Ae. albopictus density can be avoided [33]. In addition, this study was conducted to evaluate the spatial distribution of Ae. albopictus by indicator Kriging and co-Kriging based on the yearly data, which can provide an overview of high-risk areas. However, hotspots and impact levels both vary on different scales for different reasons [48].
This study was limited by the lack of combination of the breeding sites and habitats of Ae. Albopictus, such as the Breteau index. Future studies should include the preferred larva breeding sites and adult habitats, as well as NDVI (Normalized Difference Vegetation Index) and the human population in spatial modeling of abundance, which will increase the accuracy and comprehensiveness of the model.

Conclusions
In conclusion, the new surveillance system with MOT based on grids is affordable and can predict areas of high and low vector density in Shanghai. This new surveillance system can improve our understanding of the spatial and temporal distribution of Ae. Albopictus in urban regions of Shanghai and is a practical method for decision-makers to use for vector control and management of mosquitoes. In future studies, the exploration of larger scale use of this monitoring program is needed, and more data are needed to validate this new method. The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.