Microsatellite DNA Analysis Revealed a Drastic Genetic Change of Plasmodium vivax Population in the Republic of Korea During 2002 and 2003

Background Vivax malaria was successfully eliminated in the Republic of Korea (South Korea) in the late 1970s, but it was found to have re-emerged from 1993. In order to control malaria and evaluate the effectiveness of malaria controls, it is important to develop a spatiotemporal understanding of the genetic structure of the parasite population. Here, we estimated the population structure and temporal dynamics of the transmission of Plasmodium vivax in South Korea by analyzing microsatellite DNA markers of the parasite. Methodology/Principal Findings We analyzed 14 microsatellite DNA loci of the P. vivax genome from 163 South Korean isolates collected from 1994 to 2008. Allelic data were used to analyze linkage disequilibrium (LD), genetic differentiation and population structure, in order to make a detailed estimate of temporal change in the parasite population. The LD analysis showed a gradual decrease in LD levels, while the levels of genetic differentiation between successive years and analysis of the population structure based on the Bayesian approach suggested that a drastic genetic change occurred in the South Korean population during 2002 and 2003. Conclusions/Significance Although relapse and asymptomatic parasite carriage might influence the population structure to some extent, our results suggested the continual introduction of P. vivax into South Korea through other parasite population sources. One possible source, particularly during 2002 and 2003, is North Korea. Molecular epidemiology using microsatellite DNA of the P. vivax population is effective for assessing the population structure and temporal dynamics of parasite transmission; information that can assist in the elimination of vivax malaria in endemic areas.


Introduction
Plasmodium vivax, which is the second most prevalent species of the human malaria parasite, is widely distributed around the world especially in Asia, Melanesia, the Middle East, South and Central America. There are 390 million cases reported annually and, as of 2009, there were 2.85 billion people at risk of infection [1,2]. Treatment and/or prophylactic failures of chloroquine to vivax malaria have been reported from several endemic countries since the late 1980s [3]. Moreover, treatment failures of primaquine to hypnozoites of P. vivax have also been reported, although it should be noted that it is difficult to discriminate between reinfection and parasite tolerance [3]. Severe cases of vivax malaria were also reported in endemic areas [4,5]. Thus, P. vivax deserves more attention than it has previously received [6].
Understanding the genetic characteristics of the malaria parasite population is important for the monitoring of the transmission pattern and evaluation of the effectiveness of malaria control in endemic areas [7][8][9][10]. The transmission dynamics and population structures of P. vivax have recently been reported in some tropical and subtropical areas where the parasites are prevalent year-round or prevalent in a particular season but continuous during the year [11][12][13][14][15][16][17][18][19]; much less is known, however, about these characteristics in temperate areas where vivax malaria is seasonally prevalent but not continuous throughout the year [20][21][22][23][24][25].
In the Republic of Korea (South Korea), which is located in the temperate zone of Asia, indigenous vivax malaria was successfully eliminated by the late 1970s due to an effective program implemented by the National Malaria Eradication Service of the South Korean government with support from the WHO [26][27][28].
There has, however, been a re-emergence since 1993 [29]. The only patients at the beginning of the re-emergence were soldiers from the South Korea and the United States or veterans who had served in the western Demilitarized Zone, which is in the border area between the Democratic People's Republic of Korea (North Korea) and South Korea [30][31][32]. However, the number of infected civilians who lived in or near the area gradually increased [30], which suggested that P. vivax was being transmitted locally between humans and Anopheles mosquitoes. In spite of continuous In a previous study, we assessed the population structure and temporal dynamics of P. vivax transmission in South Korea using the allelic data pertaining to 10 microsatellite DNA loci of 87 isolates (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008), and highlighted the differences between the parasite populations in tropical and temperate regions [33]. We found a gradual increase in the levels of genetic diversity and a decrease in the levels of multilocus linkage disequilibrium (LD) from 1994 to 2008. Choi et al. [34] and Honma et al. [35] reported similar tendencies. Furthermore, we found that two main groups of haplotypes had been transmitting for nearly 10 years. This suggested that in temperate South Korea, the recombination rate of P. vivax was lower than in tropical areas in which genetically identical haplotypes (clones) were seldom seen in two consecutive years. Despite the low recombination rate, other new haplotypes have been observed since 1997. These new haplotypes are genetically distinct from the main two. Taken together, these results suggest that P. vivax is being continually introduced from other population sources.
Our previous study was limited by its relatively small sample size per year (average: 5.8 isolates/year) and some of the sample sizes during the year 2004-2008 were particularly small (2 to 7 isolates, average: 3.6 isolates/year). This was, perhaps, too small to discuss the temporal dynamics of parasite transmission in detail. Thus, in the present study, we increased the number of microsatellite DNA loci (from 10 to 14 loci) as well as the sample size of P. vivax parasites (average: 10.9 isolates/year) from South Korea, and performed population genetic analysis, with a focus on the differences of the parasite populations between successive years. Through this, we aimed to provide a more detailed and precise estimate of the characteristics of the vivax malaria population structure and the temporal dynamics of its transmission. Based on our new findings, we discuss the situation of the P. vivax population in South Korea and provide a possible explanation as to why, in spite of a continuous malaria control program, efforts to eliminate vivax malaria have been unsuccessful.

Ethics statement
Ethical approval for the use of South Korean P. vivax samples in this study was obtained from the Inje University Busan Paik Hospital Institutional Review Board, Busan, Korea, and performed in accordance with the Ethical Guidelines for Clinical Research issued by the Ministry of Health, Labour and Welfare of Japan on July 31, 2008, and the Ethical Guidelines for Epidemiological Research issued by the Ministries of Health, Labour and Welfare, and of Education, Science, Culture, and Sports of Japan on December 1, 2008. Written or oral informed consent from the patients could not be obtained at each sample collection for the specific purpose of this study because of the longterm prior collection of widely distributed samples. However, no author of the study was involved in gathering patient samples and the donors' personal information was disconnected from the authors. Thus, all of the samples were anonymous and it is doubtful that the analysis of the results would constitute a breach of donor privacy.

Materials
Two hundred and fifty-five P. vivax samples isolated from South Korean soldiers or veterans who had served in the DMZ from 1994 to 2008 were used in the present study. These cases were also diagnosed through microscopic examination of peripheral blood smears after the patients contracted malaria. Blood samples were collected and preserved until use at 230uC.

DNA extraction
Parasite DNA was extracted from frozen whole blood samples by phenol-chloroform extraction after proteinase K digestion [36] or by QIAamp DNA Mini Kit (Qiagen, Valencia, CA, USA) in accordance with the manufacturer's instructions.  Table S1. Four of the 14 microsatellite loci were new markers (MS2, MS3, MS10 and MS16), and were not used in the previous study [33]. PCR primer sets and amplification conditions were consistent with the protocol established by Karunaweera et al. [37]. Fluorescencelabeled PCR products were measured using an Applied Biosystems

Author Summary
Vivax malaria is widely prevalent, predominantly in Asia and South America. There were 390 million cases reported in 2009. Worldwide, in the same year, 2.85 billion people were at risk of infection. Plasmodium vivax is not limited to tropical and subtropical regions; it also appears in temperate areas. We previously examined the characteristics of P. vivax in South Korea, a temperate area, temporally, using 10 microsatellite DNA (a short tandem repeat DNA sequence) in the parasite genome and found the population to have low genetic diversity and a low recombination rate in comparison to tropical areas. In the present study, we examine the successive changes of the South Korean populations from 1994 to 2008, and the reasons for the lack of success in eliminating vivax malaria in this country. We found that, in spite of a low recombination rate, outbreeding between different genotypes seemed to increase gradually, and that the genetic composition of the population drastically changed during 2002 and 2003. This suggests that the parasite is continually introduced from other populations, probably from North Korea. The present study also demonstrated the utility of polymorphic DNA markers of P. vivax for the estimation of the transmission situation in endemic areas.
Prism Genetic Analyzer 3130xl with GeneMapper (R) version 4.1 and a 500 ROX size standard (Applied Biosystems, CA, USA).
Amplified different-sized PCR products that used the same primer sets were considered to be individual alleles within a locus, because the variation in size among isolates was consistent with the repeat number in a microsatellite locus [8]. The electropherogram shows peak profiles for the microsatellite loci, based on the fluorescence intensity of the PCR products labeled in this analysis. Multiple alleles per locus were scored if minor peaks were taller than at least one-third of the height of the predominant allele for each locus on electropherograms. Multiple-genotype infections (MGIs) were defined as those in which at least one of the 14 loci contained more than one allele [8].

Multilocus linkage disequilibrium of the P. vivax population in South Korea
Multilocus linkage disequilibrium (LD) was assessed based on the allelic data of 14 microsatellite DNA loci, using the standardized index of association (I A S ) [38,39]. This analysis was performed using the LIAN 3.5 Web interface software [40]. I A S was calculated with the formula I A S = (V D /V e 21)/(l21), with permutation testing of the null hypothesis of the complete linkage equilibrium (I A S = 0), where V D is the observed mismatch variance, V e is the expected mismatch variance and l is the number of examined loci. Significance levels of the observed I A S values were calculated by Monte-Carlo simulation, using 10,000 random data permutations. This statistical method is a variation of that proposed by Smith et al. [38]. The results were standardized by the number of loci in order to enable a comparison of different data sets [33]. This test was applied to the whole data sets in two ways: 1) all haplotype data of the 163 isolates were analyzed, thereby providing complete confidence in the haplotype profile; 2) multilocus genotypes found in multiple isolates were only counted once in the analysis, such that only unique haplotypes were counted, slightly reducing the sample size and removing the possible effect of recent epidemic expansion of particular clones [8].

Genetic differentiation of the P. vivax populations in South Korea between successive years
The extent of population subdivision for all pairs of years in South Korea was estimated based on the allelic data of 14 microsatellite DNA loci, using Weir and Cockerham's Ø estimator for determining F statistics (F ST ) [41]. F ST was calculated using the software FSTAT version 2.9.3.2 and tested for significant difference from 0, based on 1,000 random permutations of the data set [42].

Population structure analysis of the South Korean P. vivax populations by Bayesian approach
The P. vivax population structure in South Korea (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) was estimated based on the allelic data of 14 microsatellite DNA loci using the Bayesian approach model with STRUCTURE version 2.3.4 software [43]. This clustering method assigns samples to the true number of clusters (K) according to allele frequencies of each locus. Structure analysis was performed with 10 runs for each of 10 different K values (1 to 10), with a length of 50,000 burn-in periods followed by 100,000 Markov Chain Monte Carlo replications. The admixture model, with the assumption of both correlated and uncorrelated allele frequencies among the populations, was used in this analysis. Population numbers were inferred by plotting the log probability of the data [Ln P(D)] for each K value. Moreover, in order to estimate the most appropriate number of K (true K), DK, which is the rate of change in the log probability of the data between successive K values, was calculated as described by Evanno et al. [44].

Genotyping of allelic data
The allelic data of the 14 microsatellite loci were obtained from 163 of the 255 (63.9%) isolates from 1994 to 2008 (average: 10.9 isolates/year) that were used in the study. Allelic data of the 14 loci from the remaining 92 isolates (36.1%) were unavailable or only partially available due to failure in acquiring PCR products from some loci by PCR-based genotyping. Failure was possibly due to there only being a small amount of DNA for PCR amplification or because the quality of the DNA was low after multiple incidences of freeze-thaw. Figure 1 shows the number of isolates found in samples from each year (Fig. 1). In total, 71 unique haplotypes were observed in the 163 isolates (  14) and isolates (from 87 to 163), the number of unique haplotypes was increased from 40 to 71 [33], which enabled us to perform a more precise estimation of the transmission dynamics of the South Korean P. vivax population over the study period. When compared to the previous study, however, the changes in the levels of genetic diversity (heterozygosity) observed in the present study were relatively minor (Table S1, Fig. S1 and S2) [33].
Instances where different sizes of alleles were observed in one locus, were regarded as multiple genotype infections (MGIs)which were observed in some of the 14 microsatellite loci in 160 of the 163 isolates (98.2%). The frequencies of MGIs varied among the 14 loci (0.02 to 0.74; average: 0.22) ( Table 2). We also examined the number of MGI loci per isolate. In the 163 isolates with 14 loci, the highest frequency of MGI loci per isolate was 2 (58 isolates). The frequencies decreased gradually according to the increase in the number of MGI loci (Fig. 2). The highest observed number of MGI loci per isolate was 11 (one isolate). The major alleles in each locus were used for population genetic analysis.
Level of multilocus linkage disequilibrium of the South Korean P. vivax population First, we estimated the level of multilocus LD using allelic data from all of the South Korean isolates (Table 3). When all haplotypes and unique haplotypes were both used in the analysis, significant multilocus LD (P,0.001) was observed in the P. vivax population of South Korea. Second, we calculated temporal changes to the levels of multilocus LD for each year, or 2-3 successive years when the number of isolates per year was less than 10 ( Table 4). Significant multilocus LD (P,0.001) was again observed for all groups of isolates (Table 4). However, when all haplotypes were used in the analysis, the levels of multilocus LD gradually decreased over the study period: e.g., . These tendencies were also observed when unique haplotypes were used (Table 4). Significant multilocus LD (P,0.001) was observed when all haplotypes and unique haplotypes were both used. Thus it can be postulated that there has not been a recent epidemic expansion of any particular clones in the P. vivax population in South Korea [7,8].  (Table 5). Next, because the number of isolates in some years was relatively small (,10 isolates/year), we combined allelic data of the isolates of 2 or 3 successive years in order to increase the sample size of each group (.10 isolates/group), and reevaluated the levels of genetic differentiation. After the allelic data were combined, significant genetic differentiation was observed between the populations of   [43]. Each isolate was represented by a vertical line partitioned into colored segments in proportion to the estimated membership (ancestral population). Different colors represent different genotypes (ancestral population). If one bar has more than two colors, it indicates that the genotype of the isolate is admixed with more than two ancestral populations. Results shown were for K = 2, 3 and 5. The numbers in parentheses after collection year represent the number of isolates for each year. doi:10.1371/journal.pntd.0002522.g001 Population structure analysis of the South Korean P. vivax population by Bayesian approach The results of STRUCTURE analysis using an independent allele frequency model are shown in Figures 1 and 3. The highest value of Ln P(D) was observed at K = 5 (Fig. 4 A). This result suggested that the most appropriate population number (K) of P. vivax in South Korea is 5. However, the differences between Ln P(D) values of each K (3-10) were very small. In order to make the break in slope of the distribution of Ln P(D) salient at true K, we calculated value of DK using the method described by Evanno et al. [44]. The highest value of DK was observed at K = 2, and the second highest was observed at K = 3 (Fig. 4 B). These three K values (K = 2, 3 and 5) were thus adopted ( Fig. 1 and 3). For K = 5, two genotypes (light green, yellow) were predominant from 1994 to 2002. Five genotypes (light green, blue, purple, red, yellow) coexisted from 2003 to 2008 ( Fig. 1 and 3). For K = 2, two genotypes (light green, yellow) coexisted from 1994 to 2008. The light green genotype was predominant from 1994 to 2000, whereas the yellow genotype was predominant from 2001 to 2008 ( Fig. 1 and 3). For K = 3, two genotypes (light green, yellow) coexisted from 1994 to 2001, while three genotypes (blue, light green, yellow) coexisted from 2002 to 2008 ( Fig. 1 and 3). For all K values greater than 2, there was a drastic genetic change in the population structure of P. vivax in South Korea during 2002 and 2003. The results of STRUCTURE analysis, which used a model of correlated allele frequency among population, were almost the same as those of the analysis using the independent allele frequency model.

Discussion
This is the second report on our 15-year-long longitudinal study on P. vivax population genetics in South Korea using highly polymorphic neutral markers. In the present study, we assessed characteristics of the population structure and temporal dynamics of parasite transmission in a more detailed and precise manner than was employed in our previous study, using the allelic data of 14 microsatellite DNA loci of the 163 isolates in South Korea (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008).
The results of FSTAT analysis were basically consistent with the results of the STRUCTURE analysis based on the Bayesian approach ( Fig. 1 and 3) [43,44].
It is noteworthy that some negative F ST values were observed in this study. According to Á rnason and Pálsson, a negative F ST value indicates great differences between two random individuals from the same population, rather than between two random individuals from different populations [45].   We estimated the population number, population structure and temporal dynamics of P. vivax transmission in South Korea using the STRUCTURE version 2.3.4 software [43,44]. However, the estimation of the population number (K) using this method was problematic. When the log probability of the data [Ln P(D)] was used for estimating the most appropriate population number (K), 5 (K = 5) showed the highest value of Ln P(D) (Fig. 4 A). However, when DK was used, which is the rate of change in log probability, 2 (K = 2) showed the highest value of DK (Fig. 4 B). As a result, it was not possible to estimate true K in this analysis.
In our previous study using eBURST analysis [46], we found that two predominant haplotypes were observed from 1994 to 2005 and that new haplotypes had also been appearing since 1997 [33]. In a basic sense, our present study supports these previous findings. When the population number was assumed to be 2 (K = 2), two different genotypes (light green, yellow) coexisted over the study period, but the frequencies of the two genotypes changed around 2001: the light green genotype was predominant until 2000, while yellow has been predominant since 2001 ( Fig. 1 and  3). In contrast, when the population number was assumed to be greater than or equal to 3 (K$3), two different genotypes coexisted until 2001, and thereafter, new genotypes appeared and were predominant from 2002. In the present study, when K was $3, the two genotypes that were predominantly observed until 2001 corresponded to groups 1 and 2 of the P. vivax populations in South Korea of the previous study [33].
Honma et al. also conducted STRUCTURE analysis of the P. vivax population in South Korea (1997)(1998)(1999)(2000)2007) using allelic data of 13 microsatellite DNA markers [35]. The Honma et al. study used 9 out of the 13 microsatellite DNA markers that were used in the present study (MS3, 5,6,8,9,10,15,16,20).. They observed the highest value of Ln P(D) when K was 5 (K = 5) and also the highest value of DK when K was 3 (K = 3). For these reasons, these two K values were adopted. They found a drastic change in the population structure of P. vivax between the groups of 1997-2000 and 2007 at both K values. The latter group showed an increase in genetic diversity. However, the year in which the parasite population changed could not be estimated because of their sampling limitation. In the present study, we estimated that the parasite population structure in South Korea underwent a drastic change during 2002 and 2003. This change of the population structure was observed at K = 3-10.
Some of the previous studies reported that the level of genetic diversity in the South Korean P. vivax population has been increasing in recent years [33,34]. However, the multilocus LD levels of the population were shown to have been gradually decreasing:  (Table 4). These results were consistent with the results of the STRUCTURE analysis in that the variation of the alleles in the population suddenly increased from 2002. These changes might enhance outbreeding between different genotypes in the population.  In the case of P. falciparum populations, Tanabe et al. observed a strong negative correlation between ''within-population genetic diversity'' and ''geographic distance from sub-Saharan Africa over Africa, Asia, and Oceania'' [47]. In contrast, other studies suggested that the levels of genetic diversity of the parasite populations are higher in high-transmission areas and lower in low-transmission areas [7,8]. Some exceptions, however, have been reported [48]. In contrast, a negative correlation as was observed with P. falciparum, has thus far never been reported with P. vivax. Likewise, there have been no reports of a correlation between the levels of genetic diversity and the levels of malaria endemicity in the P. vivax population. Relatively high levels of genetic diversity are reported even in low transmission areas, such as Sri Lanka [37]. This feature could be attributed to relapse, owing to hypnozoites in the liver of the vivax malaria patient.
In South Korea, the number of vivax malaria cases fluctuated between 864 and 2,227 cases per year from 2004 to 2011. This suggested dramatic changes in endemicity in the endemic area [1]. Moreover, previous studies suggested that the recombination rate of the P. vivax population in South Korea was very low [33]. Therefore, the increase in genetic diversity and the decrease of LD levels of the P. vivax population in South Korea cannot be attributed to increased outbreeding between different genotypes. Thus, we would prefer to attribute it to the introduction of a genetically different P. vivax population (or populations) from other sources.
The majority of malaria patients were South Korean soldiers or veterans who served in or near the DMZ, and some civilians who lived in the same area. The width of the DMZ is 4 kilometers or less [49,50]. Given that the Anopheles sinensis, the main vector  species of vivax malaria in South Korea has been demonstrated to fly up to 12 km in one night [49], the geographical origin of the other population sources is probably North Korea. Furthermore, the reported number of vivax malaria cases in North Korea was many times higher than that in South Korea [1]. Why did the parasite population structure in South Korea change drastically  during 2002 and 2003 and what happened to the parasite population in that time? One possible scenario is that the genetic diversity of the parasite population in North Korea increased and the consequent outbreeding between different genotypes could have also occurred frequently in the endemic area in 2001, because the number of P. vivax cases in the North Korea was the highest (300,000 cases) in that year [1]. Then, this genetically diverse population could have been introduced into South Korea by Anopheles mosquitoes flying over the DMZ in 2001. Thereafter, those new genotypes started to appear in South Korean patients from 2002 due to the 8-13-month incubation period of P. vivax on the Korean peninsula [26]. That the effect of high transmission in North Korea was not reflected in higher transmission in South Korea in the following year is probably due to the use of chloroquine and primaquine for chemoprophylaxis by South Korean soldiers [50]. Although the proportion of civilian patients in South Korea increased, in 2001 and 2002, about 60% of the patients were soldiers or veterans [51]. Therefore, the number of vivax cases in South Korea did not increase during 2002 and 2003 but the genetic profile of the population changed drastically, owing to the results of genetic drift such as the founder effect (introduction by infected Anopheles mosquitoes) and/or genetic bottleneck (by chemoprophylaxis).
Chemoprophylaxis was provided to South Korean soldiers before 2001, but not to all soldiers [50]. For example, chloroquine prophylaxis was emphasized for the soldiers who served in Paju, Yeoncheon, and Cheolwon counties, whereas it was not provided to the soldiers who served in Pocheon and Hwacheon counties in 1999. These differences might influence the population structure of P. vivax in South Korea. Moreover, relapse and asymptomatic parasite carriage might play an important role in the population change. However, such factors have been implicated in the parasite population since 1993, with the reemergence of vivax malaria in South Korea. Thus, relapse and asymptomatic parasite carriage can be thought to have had little influence on the population change. The possibility exists, although it is mere speculation, that both the drastic change in the number of North Korean vivax malaria cases and the use of chemoprophylaxis by South Korean soldiers have influenced the structure of the parasite population in South Korea.
In conclusion, molecular genetic epidemiology using polymorphic DNA markers of the P. vivax genome is a very useful tool for assessing the population structure and temporal dynamics of the transmission of the parasites, the knowledge of which may lead to the effective control of vivax malaria in endemic areas.