OUP user menu

Geographic Distribution and Ecology of Potential Malaria Vectors in the Republic of Korea

(CC)
Desmond H. Foley, Terry A. Klein, Heung Chul Kim, William J. Sames, Richard C. Wilkerson, Leopoldo M. Rueda
DOI: http://dx.doi.org/10.1603/033.046.0336 680-692 First published online: 1 May 2009

Abstract

Environmental geospatial data and adult and larval mosquito collection data for up to 106 sites throughout the Republic of Korea (ROK) were used to develop ecological niche models (ENMs) of the potential geographic distribution for eight anopheline species known to occur there. The areas predicted suitable for the Hyrcanus Group species were the most extensive for Anopheles sinensis Wiedemann, An. kleini Rueda, An. belenrae Rueda, and An. pullus Yamada, intermediate for An. sineroides Yamada, and the most restricted for An. lesteri Baisas and Hu and the non-Hyrcanus Group species An. koreicus Yamada and Watanabe and An. lindesayi Yamada. The relative vectorial importance of these species is unknown, and all, except An. koreicus and An. lindesayi, are predicted to occur widely in the northwest of the ROK where malaria transmission has been sporadic since its resurgence in 1993. Our ENMs suggest that it is unlikely that An. koreicus and An. lindesayi are vectors, but we do not document consistent geographic differentiation that might incriminate any of the other species as vectors. Because all species are predicted to occur in North Korea, we also cannot reject the hypothesis that malaria infected mosquitoes from North Korea may have been the cause of the resurgence of malaria in the ROK. Ecological differentiation of the eight species is inferred from collection locations and 34 environmental layers based on remote sensing and global climatic averages. Interspecific differences were noted, and characterizing mosquito habitats by ground-based and remote sensing methods is proposed.

  • ecological niche modeling
  • distribution
  • malaria
  • ecology
  • Korea

Our understanding of malaria transmission risks and the types of vector control needed for any given location relies in part on what is known, or assumed, about the geographical distribution of each of the potential vector species. The availability of accurate distribution maps for vector species, in combination with demographic and malaria diagnosis data, could assist in determining the extent of the "Mal-area" where malaria transmission is possible (Foley et al. 2008a).

The recent availability of tools for ecological niche modeling (ENM) has made it possible to outline the ecological niche and potential geography of a variety of mosquito vectors (Levine et al. 2004, Moffett et al. 2007, Foley et al. 2008b). Ecological niches are modeled using point-occurrence data as they relate to digital geospatial data layers summarizing environmental variation. Buermann et al. (2008) found that a combination of remote sensing and climatic layers resulted in better ENM performance than using these layers separately. We were interested in using remote sensing and climatic layers to model the distribution of potential malaria vectors in the Republic of Korea (South Korea, ROK). Matching mosquito distribution to foci of malaria transmission could help incriminate those species that are the primary malaria vectors (Peterson 2006a).

Before the 1950s, Plasmodium vivax malaria was endemic and widespread in the ROK (Ree 2000) and was eventually eradicated in 1979 (WHO 1981). In 1993, malaria re-emerged and reached a peak of 4,142 cases in 2000 before falling to 774 cases in 2004 (Lee et al. 2007). Most malaria cases seem to have been contracted near the demilitarized zone (DMZ) that separates North Korea (Democratic Peoples Republic of Korea, DPRK) and the ROK (Lee et al. 2007). One hypothesis to explain the resurgence of malaria in the ROK is that infected mosquitoes traveled across the DMZ from the DPRK (Chai et al. 1994, Kho et al. 1999, Ree 2000). In this case, the vector species responsible for the resurgence of malaria would likely be present on both sides of the DMZ and be able to traverse this environment, possibly feeding off nonhuman hosts. A recent competing hypothesis is that malaria arose spontaneously in both South and North Korea as a result of a combination of troop build-up along the DMZ, increased exposure to local mosquitoes, and the occasional relapse caused by hypnozoites (Hulden and Hulden 2008).

The anopheline fauna of the ROK is relatively well resolved taxonomically (Wilkerson et al. 2003, Li et al. 2005, Rueda 2005) and includes six species of the Hyrcanus Group: Anopheles sinensis s.s. Wiedemann; An. pullus M. Yamada (=An. yatsushiroensis); An. lesteri Baisas and Hu (=An. anthropophagus); An. sineroides S. Yamada; An. kleini Rueda; An. belenrae Rueda; and two non-Hyrcanus Group species An. lindesayi japonicus S. Yamada and An. koreicus S. Yamada and Watanabe. Not all of these species can be identified on morphology, and a polymerase chain reaction (PCR) technique has been developed for species identification (Li et al. 2005). Some information has been collected on the distribution and larval habitat characteristics of molecularly identified anophelines of the ROK (Rueda et al. 2006, Kim et al. 2007, Sames et al. 2008). Previous attempts to incriminate the malaria vector species are flawed because of an imperfect understanding of the anopheline fauna in the ROK. Historically, An. sinensis was considered the primary vector. However, the discovery of additional species and results from field and laboratory parasite studies have suggested that An. kleini and An. pullus are the primary vectors and An. sinensis is a secondary vector (Lee et al. 2007). However, because further mosquito and parasite sampling is needed, we regard here all species as candidate vectors.

Our aim is to present ENM-based distributional predictions for the total anopheline species known from the ROK rather than for any subset of species. A multiple species approach is potentially more informative than a single species approach but has not been commonly undertaken until recently (Barros et al. 2007, Moffett et al. 2007, Rosa-Freitas et al. 2007). The resulting models of potential species distributions can be compared with areas of suspected malaria transmission to better understand the historical and current distribution of malaria risks in the ROK. In particular, we are interested in determining whether geographical differentiation might incriminate any one species as a vector or the species responsible for malaria resurgence. The ecological requirements of each species were studied using data derived from remote sensing and global climatic averages available for each collection location to provide further insight into species distribution.

Materials and Methods

Species Occurrence Data.

The occurrence points comprised 152 unique sites for An. sinensis, 51 for An. pullus, 54 for An. kleini, 49 for An. belenrae, 27 for An. sineroides, 13 for An. koreicus, 16 for An. lesteri, and 10 for An. lindesayi. These were for larval and adult collections from 1998 to 2006, some of which have been published (Rueda et al. 2006, Sames et al. 2008). Species were identified using PCR as described by Li et al. (2005). Clustering of sample points was evident, with most occurring in the northwest of the ROK, near the DMZ. To minimize the effects of spatial autocorrelation, only localities separated by at least 5 km were used, because it was considered that these sites would exhibit sufficient potential variation to be considered spatially independent. Thus, only a total of 106 points were used in the final analysis, i.e., 80 (An. sinensis), 33 (An. pullus), 32 (An. kleini), 36 (An. belenrae), 20 (An. sineroides), 9 (An. koreicus), 12 (An. lesteri), and 9 (An. lindesayi). For the purposes of model evaluation (see below), separate analyses were done for species with low numbers of sample points (<25) and those with high numbers (≥25). All occurrence locality data used in this study are available at http://www.mosquitomap.org/ with the global unique identifiers starting with MMap:KS.

Environmental Data.

We obtained GIS raster data layers summarizing elevation (above sea level) and a selection of 30 arc second resolution climate layers for 1950–1990 from the WorldClim dataset (Hijmans et al. 2005, http://www.worldclim.org/). Climate layers included monthly temperature (maximum, minimum, mean), monthly precipitation, and 19 "bioclimatic" layers (representing annual trends, seasonality, and extreme or limiting environmental factors). These layers were: annual mean temperature, mean diurnal range [mean of monthly (max temperature - min temperature) = P2], isothermality [(P2/P7) × 100], temperature seasonality (SD × 100), maximum temperature of the warmest month (P5), minimum temperature of coldest month (P6), temperature annual range [(P5 - P6) = P7], mean temperature of the wettest quarter, mean temperature of the driest quarter, mean temperature of the warmest quarter, mean temperature of the coldest quarter, annual precipitation, precipitation of the wettest month, precipitation of the driest month, precipitation seasonality (coefficient of variation), precipitation of the wettest quarter, precipitation of the driest quarter, precipitation of the warmest quarter, and precipitation of the coldest quarter. Temperature is in degrees Celcius (×10) and precipitation in millimeters.

We also obtained five layers summarizing topography and landform (compound topographic index, slope, aspect, flow direction, and flow accumulation) from the U.S. Geological Survey’s Hydro-1K Elevation Derivative Database (http://edc.usgs.gov/products/elevation/gtopo30/hydro/, native resolution 1 by 1 km). The compound topographic index layer summarizes upward curvature of land surface and consequent water pooling, which is important for mosquito development. Aspect, measured in degrees, can be thought of as slope direction with zero at north and increasing clockwise. Flow accumulation at this resolution translates directly into upstream drainage area in square kilometers.

We included data layers summarizing the monthly photosynthetic capacity or "greenness index" termed the normalized difference vegetation index (NDVI; 0.0091° resolution) (Tucker 1979) from the advanced very high resolution radiometer (AVHRR) satellite (http://glcf.umiacs.umd.edu/data/landcover/) for 1992–1993. Data layers summarizing percent tree cover and 13 classes of land use/land cover (0.0083° resolution) were obtained from the Global Land Cover Facility (http://glcf.umiacs.umd.edu/data/), derived from data for 1992–1993 and 1981–1994, respectively. For the ecological comparisons, we included data on soil taxonomy suborders (0.0333° resolution), obtained from the U.S. Department of Agriculture National Soils Conservation Service (http://soils.usda.gov/use/worldsoils/mapindex/order.html). Udults (Order: Ultisols) and Orthents (Order: Entisols) are the most common soil types in the study area. According to Anon. (1999), Udults are more or less freely drained, relatively humus-poor Ultisols that have an udic (humid) moisture regimen. Udults occur in humid climates, usually receive well distributed rainfall, are often forested in a natural state, and are often cultivated. Orthents are commonly on recent erosional surfaces, which may be because of geological processes or induced by cultivation, mining, or other factors. Orthents occur in any climate and under any vegetation but usually do not occur in areas that have aquic (saturated) conditions, as occurs with a high water table.

All environmental data layers were resampled with bilinear interpolation (Phillips et al. 2006) to a pixel resolution of 0.0083° or 1 by 1 km. Layers were clipped to an area encompassing the area of interest on the Korean peninsula (Fig. 1; 33.0877–39.4875° N, 124.3686-129.885° E), for a total of 508,416 pixels per layer. Because overfitting is more likely in highly dimensional environmental spaces (Peterson et al. 2007a), we used a principal components analysis (PCA; based on correlation matrix; in Minitab 15.1.1.; Minitab, State College, PA) of the 86 environmental layers for ENMs to create 15 new axes summarizing overall variation in fewer, independent dimensions. Distribution models are available at http://www.mosquitomap.org/.

Fig. 1

Predicted lowest presence threshold distribution of An. belenrae, An. kleini, An. koreicus, An. lesteri, An. lindesayi, An. pullus, An. sinensis, and An. sineroides in Korea using the genetic algorithm for rule-set prediction (Garp). Location of known occurrences of each species is shown by black dots.

ENMs.

We used the genetic algorithm for rule-set prediction (Garp; Stockwell and Noble 1992) and a maximum entropy approach (Maxent; Phillips et al. 2006) for ENM development. These algorithms were selected because both have been used in ENM development for mosquitoes, and both have seen intensive analysis of performance (Peterson et al. 2007a, Phillips 2008).

We used a desktop version of Garp (Desktop Garp 1.1.3; Scachetti-Pereira 2003). Garp uses an iterative process of rule selection, evaluation, testing, and incorporation or rejection. A rule is selected and applied to the training data (50-75% of points input into the program). The genetic algorithm in Garp allows rules to "evolve" to maximize predictive accuracy. The change in predictivity between iterations is used to evaluate whether a particular rule should be incorporated into the model, and the algorithm runs until convergence (set at 0.01) or to a maximum of 1,000 iterations. We generated up to 1,000 models for each species and selected the 10 "best" models based on omission and commission error statistics (Anderson et al. 2003): specifically, we used an absolute omission threshold of 10% and the central 50% of the commission index distribution.

Maxent has been described and tested in detail (Phillips et al. 2006, Phillips and Dudik 2008), and we used version 3.2.1 (Phillips 2006). This version has a number of improvements, including a default logistic output format, which can be interpreted as the probability that the species is present, conditional on the environmental conditions (Phillips and Dudik 2008). Maxent is based on the idea that the best explanation for unknown phenomena will maximize the entropy of the probability distribution, subject to the constraint of the environmental conditions observed at sites where species have been detected. We used default parameters for Maxent models, except that a random seed was used.

Model Evaluation.

For Garp output of the four species with 25 or more sample points, extrinsic test data were overlaid onto predicted distributions and observed correct predictions were tallied. The observed degree of coincidence was tested according to random expectations, calculated as the product of the number of testing points and the proportional coverage of the study area predicted to be present. We followed Anderson et al. (2002) and Peterson et al. (2007b) in using cumulative binomial probabilities that the observed success could be the result of random association between predictions and test points as a test statistic. For Maxent output of the four species with 25 or more sample points, model validation statistics were based on 10 replicate random partitions of the localities into the test (25%) and training (75%) data. For each replicate, we calculated the number of test localities omitted from the prediction and applied a binomial test to check statistical significance with regard to the proportion of the study area that was predicted as present.

For the four species with <25 sample sites, for both Garp and Maxent analyses, we used the jackknife (or "leave-one-out") validation approach developed by Pearson et al. (2007), using the program made available by those authors. The more customary receiver operating characteristic analysis was not used in this study because of recent criticism and concerns about bias and artifact (Lobo et al. 2008, Peterson et al. 2008).

Model validation and interpretation requires setting a decision threshold to distinguish areas of "presence" and "absence." We follow the approach of Pearson et al. (2007) by using the lowest presence threshold (LPT). The LPT approach is relatively conservative, identifying pixels that are at least as suitable as the lowest value associated with a species’ presence (Pearson et al. 2007).

Consensus Predictions.

Because some studies (Peterson et al. 2007a) have suggested that Maxent models may be better for interpolation and Garp models better for extrapolation (i.e., transferability), we calculated a consensus of the two model outputs for analyses trained on the full data set of sample points for each species. Following Foley et al. (2008b), we calculated the area of intersection of Garp and Maxent models, but using the LPT as the decision threshold.

Variable Importance and Response Curves.

We used the program PAST 1.27 (http://folk.uio.no/ohammer/past/) to undertake a nonmetric multidimensional scaling analysis of the Euclidian distance of means of the first 15 principal components summarizing the 86 environmental data layers occurring at Korean collection locations for the eight mosquito species. This analysis is appropriate for the reduction and interpretation of large multivariate ecological data sets and was done to assist in understanding differences in the ecological requirements of these species. The algorithm in PAST 1.27 attempts to place the data points in a two-dimensional coordinate system such that the ranked differences are preserved.

We used Maxent to investigate response curves and variable importance using 33 environmental variables (bioclimatic, topography and landform, NDVI, tree cover, land cover, and elevation; Table 1).

View this table:
Table 1

Results

The results from the PCA of environmental data showed that the first principal component (PC1) was most closely related to aspects of temperature, whereas PC2 was related mainly to precipitation and seasonality and PC3 to precipitation, land cover, and NDVI (data not shown). Subsequent PCs emphasized NDVI and temperature and precipitation extremes (PC4), isothermality (PC5), land cover and topography (PC6 and PC12), precipitation and topography (PC7), topography (PC8 and PC9), precipitation (PC10 and PC11), precipitation and NDVI (PC13), slope and NDVI (PC14), and topography (PC15). In all, the first 15 PCs captured 97.4% of the overall variation in the original 86 environmental data layers.

Maxent: Species With High Numbers of Sample Points.

For Maxent predictions of An. sinensis using total point for training, the LPT score was 0.0683. For replicates where points were divided into test and training points, An. sinensis omission error ranged from 0 to 0.1 (mean, 0.02; SD, 0.035), equivalent to 0–2 test from 20 test localities omitted, with 8/10 partitions recording minimum training presence binomial probability of P < 0.05. For Maxent predictions of An. belenrae using total point for training, the LPT score was 0.1452. For replicates where points were divided into test and training points, An. belenrae omission error ranged from 0 to 0.444 (mean, 0.178; SD, 0.167), equivalent to 0-4 test from nine test localities omitted, with 6/10 partitions recording minimum training presence binomial probability of P < 0.05. For Maxent predictions of An. pullus using total point for training, the LPT score was 0.1327. For replicates where points were divided into test and training points, An. pullus omission error ranged from 0 to 0.625 (mean, 0.225; SD, 0.194), equivalent to 0-5 test from eight test localities omitted, with 5/10 partitions recording minimum training presence binomial probability of P < 0.05. For Maxent predictions of An. kleini using total point for training, the LPT score was 0.1117. For replicates where points were divided into test and training points, An. kleini omission error ranged from 0 to 0.429 (mean, 0.214; SD, 0.168), equivalent to 0-4 test from seven test localities omitted, with 4/10 partitions recording minimum training presence binomial probability of P < 0.05. Thus, Maxent distribution models for species with larger sample sizes (25 or more) indicated that An. sinensis had the best support, according to random expectations of omission error and An. kleini the least support.

Maxent: Species With Low Numbers of Sample Points.

For Maxent predictions of An. sineroides using total point for training, the LPT score was 0.2668. After Jackknife validation, An. sineroides omission error occurred in 6 of 20 replicates, with average minimum training presence area = 0.412 and P = 0.005696. For Maxent predictions of An. lindesayi using total point for training, the LPT score was 0.3982. After Jackknife validation, An. lindesayi omission error occurred in four of nine replicates, with average minimum training presence area = 0.215 and P = 0.011035. For Maxent predictions of An. lesteri using total point for training, the LPT score was 0.2412. After Jackknife validation, An. lesteri omission error occurred in 2 of 12 replicates, with average minimum training presence area = 0.034 and P = 0.833333. For Maxent predictions of An. koreicus using total point for training, the LPT score was 0.3992. After Jackknife validation, An. koreicus omission error occurred in 2 of 10 replicates, with average minimum training presence area = 0.163 and P = 0.000015. Thus, Maxent distribution models for species with smaller sample sizes (<25) indicated that all species, with the exception of An. lesteri, gave better omission error than random expectations.

Garp: Species With High Numbers of Sample Points.

For Garp predictions of An. sinensis using total point, the LPT score was threshold 5 (T5). In other words, a collection point was not predicted as present unless at least 5 of 10 of the best models agreed. All thresholds were better than random expectations (P < 0.05). For Garp predictions of An. kleini, the LPT score was T7. All thresholds were better than random expectations (P < 0.05). For Garp predictions of An. belenrae, the LPT score was T8. All thresholds were better than random expectations (P < 0.05). For Garp predictions of An. pullus, the LPT score was T8. All thresholds were better than random expectations (P < 0.05). Thus, Garp distribution models for species with larger sample sizes (25 or more) indicated that all species gave better omission error than random expectations.

Garp: Species With Low Numbers of Sample Points.

For Garp predictions of An. sineroides, the LPT score T10. After Jackknife validation, An. sineroides omission error occurred in 17 of 20 replicates, with average minimum training presence area = 0.143799 and P = 0.427249. For Garp predictions of An. lindesayi, the LPT score T10. After Jackknife validation, An. lindesayi omission error occurred in 10 of 10 replicates, with average minimum training presence area = 0.0002 and P = 1.0. For Garp predictions of An. lesteri, LPT was T10. After Jackknife validation, An. lesteri omission error occurred in 7 of 11 replicates, with average minimum training presence area = 0.016503 and P = 0.000007. For Garp predictions of An. koreicus, the LPT score was T10. After Jackknife validation, An. koreicus omission error occurred in 8 of 10 replicates, with average minimum training presence area = 0.002 and P = 0.000009. Thus, Garp distribution models for species with smaller sample sizes (<25) indicated that only An. lesteri and An. koreicus gave better omission error than random expectations.

Species models based on all data points (Figs. 1 and 2) identified distribution patterns that differed between species and algorithms. Garp predictions tended to be more extensive than Maxent predictions, as has been found in other studies (Peterson et al. 2007a). The predicted potential distribution of species more commonly collected (An. belenrae, An. kleini, An. pullus, and An. sinensis) showed the most extensive distribution. The predicted potential distribution of species least commonly collected (An. koreicus, An. lesteri, An. lindesayi, and An. sineroides) showed the least extensive distribution. Garp and Maxent models differed most noticeably for An. koreicus and An. lindesayi, which also had the lowest numbers of sample points. The Maxent model for An. koreicus identified three distinct areas of occurrence in the ROK compared with one diffuse area in the Garp model. The Maxent model for An. lindesayi predicts that this species occurs more to the south in the ROK than the Garp model (which did not significantly differ from random expectation). The Maxent model for An. sineroides predicts that this species occurs along the east coast and has a different southern distribution in the ROK compared with the Garp prediction (which did not significantly differ from random expectation). The intersection of the models using the LPT threshold is shown in Fig. 3 for an area north of Seoul, where resurgent malaria transmission has occurred. All species were predicted to occur in the northwest of the ROK and in North Korea.

Fig. 2

Predicted lowest presence threshold distribution of An. belenrae, An. kleini, An. koreicus, An. lesteri, An. lindesayi, An. pullus, An. sinensis, and An. sineroides in Korea using a maximum entropy (Maxent) approach. Location of known occurrences of each species is shown by black dots.

Fig. 3

Consensus by intersection of predicted lowest presence threshold distribution models of Garp and Maxent distributions of An. belenrae, An. kleini, An. koreicus, An. lesteri, An. lindesayi, An. pullus, An. sinensis, and An. sineroides. The area shown encompasses the zone of recent malaria transmission in Korea. Pixels where mosquito species are predicted present are colored gray except those (in black) that fall within areas (outlined) where civilian cases of vivax malaria were recorded by Kho et al. (1999) for 1994–1997. Thus, the larger the amount of black pixels the greater the evidence for that species being a vector.

Variable Importance and Response Curves.

Output distribution models using bioclimatic, topography and landform, NDVI, tree cover, land cover, and elevation variables (see Table 1) were similar to but less general than models using PCA-derived variables. Analysis of the jackknife of regularized training gain showed that the NDVI for June gave the most useful unique information (i.e., decreased the gain the most when it was omitted) for most species and the precipitation of the wettest month resulted in the highest gains when used in isolation (i.e., the most useful information by itself) for four of the eight species.

Figure 4 shows a minimal spanning tree of the nonmetric multidimensional scaling analysis of the means of the first 15 principal components for collection locations of the eight mosquito species. Average environmental conditions for the non-Hyrcanus Group species An. koreicus and An. lindesayi were different to Hyrcanus Group species. However, this analysis showed that the environment of An. lesteri was the most dissimilar of the Hyrcanus Group species.

Fig. 4

Minimal spanning tree of a nonmetric multidimensional scaling analysis of the means of the first 15 principal components summarizing 86 environmental data layers occurring at Korean collection locations of An. belenrae, An. kleini, An. koreicus, An. lesteri, An. lindesayi, An. pullus, An. sinensis, and An. sineroides

Summary data for the eight mosquito species for 1-km2 grids that coincided with collection locations is given in Table 1 for the above 33 environmental layers plus soils. Because of the potential for bias because of some species (e.g., An. sinensis) being more extensively distributed in the warmer southern regions of the ROK than other species, comparison of environmental conditions at locations in Table 1 were limited to those occurring north of latitude 36.5° N (all species were recorded from this region). Based on these data, the following environmental conditions were characteristic of the Hyrcanus Group species (compared with non-Hyrcanus Group species): higher annual mean temperature, lower mean diurnal and annual temperature range, lower annual precipitation (particularly low precipitation in the warmest and wettest quarter), lower precipitation seasonality, lower elevation, lower slope, southerly rather than south westerly aspect, lower flow accumulation (i.e., less upstream drainage), lower tree cover (except for An. sineroides), occurring in cropland and wooded grassland rather than deciduous broadleaf forest, and on Udults soil suborder types (except for An. sineroides) rather than on Orthents.

The following lists relative tendencies of environmental conditions for Korean species within the Hyrcanus Group: An. sinensis, higher NDVI in most months and not as cold in the coldest part of the year; An. belenrae, high flow accumulation (i.e., greater upstream drainage) and predominantly cropland; An. kleini, southwesterly aspect; An. pullus, coldest in the coldest part of the year, high topoindex (wetness), high flow accumulation (i.e., greater upstream drainage), low slope, and bimodal aspect; An. sineroides, low temperature seasonality and range (not as cold in the coldest part of the year and not as warm in the warmest part of the year), low precipitation seasonality (not as wet in the wettest part of the year and not as dry in the driest part of the year), low topoindex (wetness), high slope, low flow accumulation (i.e., less upstream drainage), high elevation, high tree cover, and predominantly Orthents soil type; An. lesteri, high precipitation and temperature seasonality (drier and wetter, and warmer), low mean diurnal temperature range, low isothermality, low NDVI for all months, low slope, southwesterly aspect, low tree cover (typically cropland), and low elevation.

A one-way analysis of variance (ANOVA) of data for each environmental layer was performed using Minitab for the Hyrcanus Group species. This analysis indicated a significant difference in means for isothermality (F = 3.60, df = 5, P = 0.004) and mean temperature of the wettest quarter (F = 2.48, df = 5, P = 0.034), with a nearly significant result for mean temperature of the warmest quarter (F = 2.26, df = 5, P = 0.051). In all cases An. lesteri was significantly different from the other species (P = 0.05), based on the Fisher least significant difference (LSD) method. An. lesteri has low scores for mean diurnal temperature range but high scores for temperature annual range, resulting in low scores for isothermality. This suggests that this species prefers warmer and more stable temperatures during the summer than other species within the Hyrcanus Group. Only one significant differences for the 15 PC layers was recorded (one-way ANOVA, F = 2.31, df = 5, P = 0.047) for PC5, which mainly reflects the differences between An. lesteri and the other Hyrcanus Group species in terms of isothermality.

Discussion

The species identity of malaria vectors in the ROK is unknown, so we developed ENMs for all anopheline species within the ROK to shed light on their distributions, ecology, and possible relation to malaria risk in that country. Prediction successes of distribution models were better than random except for Garp models for An. sineroides and An. lindesayi and Maxent models for An. lesteri. An. belenrae, An. pullus, and An. sinensis have been variously recorded at sites in the DPRK by Rueda and Gao (2008) (georeferencing error noted for the Pyeongseong City collection). The Kaesong City sites were positive for these three species in our predictions. Although our Garp prediction for An. lindesayi agrees with many of the sites listed in Sames et al. (2008), it did not predict any occurrence on Jeju Island. Lee (1994) reported this species on Jeju Island (Korea), but this location may be marginal, because subsequent attempts to collect this species there failed (Kim et al. 2005, Sames et al. 2008). The spatial resolution of this analysis (1 km2) may may be too coarse for smaller habitats, and the number of locations for this species available for modeling was low. Inclusion of additional records, especially at lower latitudes, should improve predictions.

One hypothesis to explain the resurgence of malaria in the ROK is that infected mosquitoes traveled across the DMZ from the DPRK (Chai et al. 1994, Kho et al. 1999, Ree 2000). Our results indicate that all eight mosquito species considered here occur on both sides of the DMZ and are candidates for a north to south border crossing. Thus, we did not document geographic differentiation that could reject the hypothesis that malaria-infected mosquitoes from North Korea may have been the cause of the introduction and resurgence of malaria in the ROK.

Did we document geographic differentiation that might incriminate any of the other species as a current or past vector? Historically, malaria was widespread in the ROK; for instance, a 1965 map from the National Malaria Eradication Service in Chai (1999) showed malaria widespread but with focal areas of greater incidence in northern Gyeongsangbuk and the north and southeast of Gyeonggi provinces. All species except An. lesteri, An. koreicus, and An. lindesayi are predicted to be extensive in both of these areas. Murdoch and Lueders (1953) reported a 0.83% malaria parasite infection rate for An. sinensis s.l. collected in the Wonju area in August 1952. At that time, morphologically identified An. sinensis could include other members of the Hyrcanus Group, but not An. sineroides, An. lindesayi, and An. koreicus. We found that all species, except An. koreicus, An. lesteri, and An. lindesayi, are predicted to occur extensively in that region. Regarding more recent malaria transmission, Fig. 3 shows areas where civilian cases of vivax malaria were recorded in 1994–1997 (from Kho et al. 1999). All species are predicted to be present where transmission is thought to have occurred, but the least co-occurrence of species and malaria was recorded for An. koreicus and An. lindesayi, with An. sineroides distribution failing to co-occur with cases in the south. These species are also infrequently collected in adult and larval surveillance (Kim et al. 2007). Thus, distribution modeling suggests that An. koreicus and An. lindesayi, at least, are unlikely to have played a significant role in past and present malaria transmission in the ROK. However, it is possible that more than one anopheline species has been involved in malaria transmission in the ROK.

Before eradication in the 1970s, malaria transmission occurred widely throughout the ROK, possibly because of local conditions, such as agricultural practices, favoring vector abundance, and survival, or because of the extent of human/mosquito contact and lack of rapid medical treatment that favors transmission. The fact that resurgent malaria has not returned to preeradication limits suggests that conditions in the ROK have changed. For example, Korean forests were badly degraded through the first half of the 20th century, but between 1967 and 1995, forested areas increased to 6.3 million ha as a result of large-scale reforestation (Food and Agriculture Organization of the United Nations Forestry county profiles, http://www.fao.org/forestry/23809/en/kor/). The seasonal incidence of malaria in 1913 versus 1993–1996 reported in Chai (1999) showed a shift, from June and August peaks to only August, which could indicate either change in the incubation period of the parasite, the involvement of different vector species, or a change in the phenology of the vector species. Recent evidence suggests that ≈60% of the malaria cases are latent with the onset of symptoms 6–18 mo after infection (T. A. Klein, personal communication).

Considering ecological requirements, Hyrcanus Group species seemed to differ from non-Hyrcanus Group species (Fig. 4), and An. lesteri was the most divergent among the former. The frequency of occurrence of An. sinensis increases to a maximum for locations with a south-facing aspect, whereas aspect is predominately south southeast and south southwest for An. pullus and in the southwest and northwest quadrants for An. lindesayi (Fig. 5). No trend with time of year could be discerned (data not shown). Could this difference be a function of species using different land cover categories, which in turn occur more frequently at different aspects? We graphed the aspect of all grid cells (0.0083° resolution) for Korea according to land cover category (Fig. 6). Forest areas (categories 1-5) were more frequently northeast to southeast facing, whereas crops and grassland areas (7–11) were more frequently west to southwest facing. The forest finding agrees with Johnson et al. (2002), who showed that site quality for oaks in West Virginia (a similar latitude to the ROK) was maximal at 81° azimuth and least favorable at 261°; increasing slope reduced site quality caused by associated increases in solar radiation. The cool and moist microenvironment in the northeast quadrant resulted in thicker leaf litter layers, which have various effects on soil development processes (Johnson et al. 2002). However, despite An. lindesayi and An. koreicus being most frequent in forest areas, they are more likely to be found in areas facing southwest, and despite An. sinensis being most frequently found in crops and grassland areas, it is most often found in land facing south. The interaction of environmental layers may be important for mosquito ecology. For example, slope and aspect would be expected to interact with tree cover that may relate to a preference for shaded habitats. As mentioned above, land cover in the ROK is dynamic. Therefore, landscape and tree cover estimates from remote sensing data drawn from the 1980–1990s may not apply to mosquito collections conducted after 2000.

Fig. 5

Frequency distribution of aspect (degrees clockwise from north) for Korean collection locations for three species (An. sinensis, An. pullus, and An. lindesayi).

Fig. 6

Frequency distribution of aspect (degrees clockwise from north) for grid cells (1 km2) in the study area of Korea, according to landscape category (1 = evergreen needle leaf forest, 4 = deciduous broadleaf forest, 5 = mixed forest, 6 = woodland, 7 = wooded grassland, 8 = closed shrubland, 9 = open shrubland, 10 = grassland, 11 = cropland, 12 = bare ground, and 13 = urban and built). Note that forest areas (categories 1–5) are more frequently east to southeast facing, whereas crops and grassland areas (7–11) are more frequently west to southwest facing. Categories 2–3 are not present.

According to Kim et al. (2007), An. kleini numbers peak during late spring and early summer, whereas An. sinensis peaks during late summer and fall. The more southerly extension of An. sinensis into tropical regions of Indochina supports the idea that this species is adapted to higher temperatures. Much of the difference in average temperatures for An. sinensis and An. kleini is caused by temperatures during the colder months. Overwintering adaptations such as diapause may not be well developed in An. sinensis. In contrast, the non-Hyrcanus Group species were found in locations with the lowest (below freezing) winter temperatures and may have well-developed overwintering adaptations. An. lindesayi have been found overwintering in the ROK as larvae, even in icy water (H.-C. Kim, personal communication).

Kim et al. (2007) found larvae of An. sinensis mainly from water seepage springs and marshes, An. pullus from rice fields and stream margins, An. belenrae from pond and lake habitats, and An. lindesayi from pools of fast running streams. Some concordance between results from larval habitat studies and remote sensing data can be seen. For example, we found that An. pullus occurred mainly in croplands with low slope, high flow accumulation, and high wetness (topoindex), as you would expect of rice field habitats, and An. lindesayi occurred in areas of high flow accumulation and high slope, as you would expect of fast running streams. Combining data on larval microenvironment and high-resolution remotely sensed data may yield powerful insights into ecological differentiation.

Maxent response variables identified NDVI and climate data as important for prediction models, thereby supporting the contention of Buermann et al. (2008) that a combination of remote sensing and climatic layers is more important for ENMs than either is in isolation. Our distribution models are based largely on yearly climate averages (e.g., WorldClim), and as such, portray general potential distributional tendencies. In addition, the realism of ENMs will depend in part on the biological relevance for mosquito biology of environmental layers such as NDVI and aspect. Absence data are less informative for assessing model accuracy because ENM shows potential distributions rather than actual distributions; other factors such as historical, physical, climatological, and biotic constraints may limit distribution (Soberón and Peterson 2005, Peterson 2006b). Improvements in model accuracy could be made by using more widely distributed sample points, absence as well as presence information, analyses by month rather than year, and by using more up-to-date, higher resolution environmental layers. For example, more detailed information on rice irrigation, tree cover, and landscape classes may improve predictions.

Although these models propose to identify the potential limits of species distributions, they do not define population densities and survival, two critical components of vectorial capacity. For example, An. kleini and An. pullus, although relatively widely distributed, in agreement with model predictions, accounted for >50% of all anophelines collected at military installations near the DMZ. South of Seoul, these two species account for <5% of the total population of mosquitoes in both larval and adult collections. The distribution of malaria transmission in the ROK has changed in the last century, but it is unknown what role changes in species distribution, if any, have played on changes in malaria distribution. Weather anomalies and disasters, including war, could affect the distribution and abundance of mosquitoes because of, for instance, changes in landscape and/or farming practices. Thus, mosquito species distribution models are only the first step in the process of modeling the entomological risk factors involved in malaria transmission.

Acknowledgements

We thank J.-P. Chretién for advice and support and A. T. Peterson and M. Papeş for advice on ENM. Funding for this project was provided by the U.S. Department of Defense Global Emerging Infections Surveillance and Response System, Silver Spring, MD (Project G00018_08_WR).

Footnotes

  • This research was performed under a Memorandum of Understanding between the Walter Reed Army Institute of Research and the Smithsonian Institution, with institutional support provided by both organizations. The published material reflects the views of the authors and should not be construed to represent those of the Department of the Army or the Department of Defense.

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0/), which permits non-commercial reuse, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions{at}oup.com

References Cited

View Abstract