Population Genetic Structuring in Opisthorchis viverrini over Various Spatial Scales in Thailand and Lao PDR

Khon Kaen Province in northeast Thailand is known as a hot spot for opisthorchiasis in Southeast Asia. Preliminary allozyme and mitochondrial DNA haplotype data from within one endemic district in this Province (Ban Phai), indicated substantial genetic variability within Opisthorchis viverrini. Here, we used microsatellite DNA analyses to examine the genetic diversity and population structure of O. viverrini from four geographically close localities in Khon Kaen Province. Genotyping based on 12 microsatellite loci yielded a mean number of alleles per locus that ranged from 2.83 to 3.7 with an expected heterozygosity in Hardy–Weinberg equilibrium of 0.44–0.56. Assessment of population structure by pairwise FST analysis showed inter-population differentiation (P<0.05) which indicates population substructuring between these localities. Unique alleles were found in three of four localities with the highest number observed per locality being three. Our results highlight the existence of genetic diversity and population substructuring in O. viverrini over a small spatial scale which is similar to that found at a larger scale. This provides the basis for the investigation of the role of parasite genetic diversity and differentiation in transmission dynamics and control of O. viverrini.


Introduction
The liver fluke, Opisthorchis viverrini is a food-borne trematode endemic in Southeast Asia, including Thailand, Lao PDR, Vietnam and Cambodia with more than 10 million people infected [1,2,3]. Infection occurs by eating raw or uncooked cyprinid fish containing metacercariae [4,5]. O. viverrini infection is a significant medical problem because of its involvement as a major risk factor causing bile duct cancer (cholangiocarcinoma, CCA) [6]. Liver cancer, predominantly CCA, is the fourth and fifth cause of mortality in males and females, respectively in Thailand [7]. Globally, Khon Kaen Province, Thailand is one of the hot spots of CCA with incidence levels (per 100,000) of 78.4 in males and 33.3 in females [8].
Recently, we reported that O. viverrini does not represent a single species but consists of at least two morphologically similar but genetically distinct (i.e. cryptic) species from Thailand and Lao PDR [9]. We also showed that there were at least six genetically distinct groups that are associated with different major wetlands. Additionally, biological variation between populations of O. viverrini from different wetlands in Thailand and Lao PDR has been detected. For instance, worm recovery as well as the fecundity of O. viverrini from the Songkram River in Thailand was significantly different from other wetland systems (Chi, Mun and Wang Rivers) in Thailand and Lao PDR (Nam Ngum River) [10]. Furthermore, worms belonging to this population were significantly different in body size from populations from the Chi and Nam Ngum River wetlands [10]. The fine scale population genetics of O. viverrini has to date only been studied from a single locality (Ban Phai in Khon Kaen, Thailand), but the results indicated considerable genetic diversity and heterozygote deficiency occurring within a small geographical area [11].
More detailed information on the population genetic structure of O. viverrini is, however, needed to fully determine whether population substructuring and/or differential genetic diversity are associated with geographical differences in distinct wetlands, river systems and flooding patterns [12]. Recently, we characterized, optimized and demonstrated the utility of microsatellite DNA markers for O. viverrini and provided evidence of population subdivision over a large spatial scale with the maximum distant apart of up to 770 km [13]. However, whether such a population pattern occurs over a small spatial scale or not is unknown.
In this study, we examined the genetic diversity and population structure of O. viverrini populations occurring within and between four geographically close localities (small geographical scale population comparisons) less than 60 km apart in Khon Kaen Province, northeast Thailand. Comparisons were also made with data previously reported for populations separated by much greater distances (widely spaced/large scale population comparisons). Finally, potential population mechanism(s) related to the host and environmental factors that may contribute to the current population structure were discussed.

Study areas and O. viverrini isolates
Khon Kaen Province is geographically centrally located in northeast Thailand and by political division it currently consists of 26 districts. Samples of O. viverrini for this study were obtained from four reservoirs located within three adjacent districts namely Khon Kaen, Ban Phai and Phu Wiang in Khon Kaen Province, Thailand as shown in Figure 1. Sampling localities from Khon Kaen district were from Ban Sa-ard (KBs) and Ban Lerngpleuy (KLp) and the other two were from Ban Phai district (KBp) and Phu Wiang district (KPv). The Ban Phai, Lerngpleuy and Ban Sa-ard localities are connected to the Chi River, whereas the Phu Wiang locality is close to the Nam Phong River upstream from Ubonratana Dam. The Nam Phong River feeds into the Chi River downstream from Ban Sa-ard ( Figure 1). The maximum and minimum geographical distance between localities is between Ban Phai and Phu Wiang (60 km) and Ban Sa-ard and Ban Lerngpluey (10 km), respectively. For comparisons of widely spaced populations, the distance between localities ranged from 225-771 km [13].
Samples of adult O. viverrini were recovered from hamsters experimentally infected with metacercariae obtained from pools of naturally infected cyprinid fish (Cyclocheilichthys armatus). In each sampling locality, approximately 500 fish weighing at least five kilograms were processed by a pepsin digestion method to isolate metacercariae [14]. Of these a random selection of a maximum of 50 metacercariae were fed orally to each hamster and four months post infection the adult worms were recovered from the biliary system of infected hamsters. Five hamsters were used per locality for worm recovery and previous studies have shown that a dose of up to 50 metacercariae per animal for a period of four months infection does not harm the hamsters and causes no morbidity compared with uninfected control hamsters [15].
The worms were identified based on standard morphological methods and washed several times with 0.85% NaCl. A minimum number of hamsters was used to provide a sufficient number of worms for different experiments that we are undertaking. For microsatellite analyses a random sample of 30 individual worms from a mean of 26-34 worms per hamster per locality were selected.

DNA extraction and microsatellite analysis
Individual worms were homogenized on ice in a microcentrifuge tube using a handmade glass pestle. Genomic DNA was extracted by GenomicPrep Cells and a Tissue DNA Isolation kit following manufacturer recommendations (GE Healthcare, NJ, USA). DNA concentration and purity were determined by spectrophotometry (Pharmacia Biotech, Cambridge, UK). Previously isolated and characterized O. viverrini microsatellite loci (12) were used in this study [13]. The forward primer of each pair was modified with fluorescent dye (6-FAM or HEX or NED; PE Applied Biosystems, CA, USA). Microsatellite analyses were performed using a Polymerase Chain Reaction (PCR) containing 1 ng of template DNA, 2 mM Tris-HCl, 10 mM KCl, 2 mM Mg 2+ , 0.2 mM of each nucleotide, 0.2 pmol of each primer, and 0.05 units Taq polymerase (Takara Biomedicals, Tokyo, JP). PCR amplifications were carried out in a BIORAD thermocycler (BIORAD, CA, USA) in a total volume of 25 ml. The cycling conditions included 30 cycles of 1 min at 94uC, 1 min at the optimized annealing temperature, and 3 min at 72uC. PCR products were diluted in HIDI formamide with internal GeneScan size standard, ROX-400HD (PE Applied Biosystems) then loaded on the ABI 3100 DNA sequencer (PE Applied Biosystems). Before analysis, the PCR products were denatured in the thermocycler at 95uC and rapidly cooled on ice. Allele sizes were determined using ABI Prism GeneScan Analysis 3.1 and Genotyper 2.5 (PE Applied Biosystems). PCR reactions were redone in all cases where samples could not be amplified and in each case non amplification was confirmed. To avoid scoring artificial bands resulting in scoring errors, all PCR of samples with electropherograms with many peaks or non-specific products were repeated until unambiguous results were obtained. Furthermore, only clear electropherograms with one or two peaks of the expected size were considered in the analysis.

Data analysis
For each locus, the number of alleles, allelic frequencies, and linkage disequilibrium among polymorphic loci using the Markov chain approach [16], and observed and expected heterozygosity were calculated [17]. To avoid genotyping errors (i.e. the presence of null alleles, large allele dropout and scoring errors due to stuttering peaks), the program Micro-Checker version 2.2.3 was used [18] to correct allele frequencies as described by Brookfield [19]. Hardy-Weinberg Equilibrium (HWE) for each locus was examined using the exact test [20]. The fixation index within subpopulations (F IS ) and genetic differentiation between populations (F ST ) based on Weir and Cockerham [21] was determined using F-statistics [22]. The significance of pairwise F ST values was evaluated [21]. The relationship between genetic isolation among localities was assessed by testing for independence between F ST and geographical distances by a Mantel test. All the calculations described above were conducted using GENEPOP Version 3.4 software [20]. Allelic richness and overall estimated F IS of the parasite populations were calculated by using FSTAT [23]. Data analyses were done by comparing small scale geographically defined populations in Khon Kaen Province (the Chi River wetland distinct genetic groups defined by Saijuntha et al. [9]).

Author Summary
Infection with the liver fluke (Opisthorchis viverrini) is a risk factor for cholangiocarcinoma (CCA), which is highly prevalent in Khon Kaen Province, Thailand. Within this Province, there is considerable variation in the epidemiology of opisthorchiasis among districts. Preliminary allozyme and mitochondrial DNA data indicate that genetic variation in O. viverrini occurs even over a small endemic area within the province. Here, we used microsatellite DNA analyses to examine the population genetic structure of O. viverrini from four geographically close localities. Analyses of adult worms based on 12 microsatellite loci revealed varying levels of genetic diversity and population substructuring between the four localities. Worms originating from one locality (Phu Wiang) had significantly higher genetic diversity than the other three localities. Data on the population genetic structure observed in these localities are similar to those found at a larger geographic scale. This provides background data to further investigate the biological and epidemiological significance of these genetic variants.
These closely associated populations were compared with more widely distributed populations (from other wetlands containing distinct genetic groups and/or cryptic species) occurring in Thailand and Lao PDR [13]. By using analysis of multilocus genotypes, genetic ancestry can be inferred regardless of the sampling location of individuals.

Ethics statement
This study was performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Research Council of Thailand. The protocol of animal experimentation was approved by the Institutional Animal Ethics Committee, Khon Kaen University (AEKKU20/2551). All surgery and necropsy was performed under sodium pentobarbital anesthesia, and every effort was made to minimize pain and suffering to the animals.

Genetic diversity
Allele distribution patterns at 12 microsatellite loci varied greatly among samples of O. viverrini (Table 1). For the populations in Khon Kaen Province, 52 alleles were recorded across all individuals (120 worms) at 12 microsatellite loci. The total number of alleles per locus ranged from 2-12. The locality with highest total number of alleles was KBp (44) and the lowest was KBs (34). Localities KLp (40) and KPv (37) had intermediate numbers.
Genotypic disequilibrium analysis indicated that genotypes of most loci were associated randomly (P.0.05) (data not shown). However, linkage disequilibrium was found for Ovms14 and 15, but this was not significant after sequential Bonferroni's adjustment (P.0.01) [24]. Within O. viverrini populations unique alleles are of low frequency (,6%). Three of the four Khon Kaen populations (KPv, KBp and KLp) exhibited unique alleles (3, 2 and 2 respectively) across five loci (Table 1). Widely spaced populations used in the earlier study [13] exhibited a total of 16 unique alleles across eight loci (Table 1). However, three alleles occurring in two Khon Kaen populations were also ''unique'' in three populations from reference 13.
With respect to genetic diversity among Khon Kaen O. viverrini populations as shown in Table 2, the average number of alleles per locus for KPv, KBp and KLp (range 3.33-3.7) were higher than that for KBs (2.83) but the difference was not significant (P.0.05). Expected heterozygosity per locus (H E ) for the four localities  To compare genetic differentiation between localities, F ST statistics were calculated and revealed significant (P,0.05) genetic differentiation between all pairs of localities at all geographic scales (Table 3). Qualitative guidelines suggested by Wright (1978) [22] were adopted, namely, F ST genetic differentiation: 0-0.05 'little'; 0.05-0.15 'moderate'; 0.15-0.25 'great'; and .0.25 indicate 'very great'. The level of genetic differentiation detected here ranged from 0.0002-0.0776 which suggests that O. viverrini in Khon Kaen Province is not panmictic but has low to moderate genetic differentiation among populations.

Comparison with populations at larger geographic scales
Data obtained from the present study concerning genetic diversity and population differentiation were compared with a previous report which analyzed five widely spaced populations of O. viverrini in Thailand and Lao PDR [13]. This was possible because the same 12 microsatellite loci were used in the analysis. As shown in Tables 2 and 3, the measurements of allelic diversity, expected heterozygosity, allelic richness and F ST from this study were not different from widely spaced populations. The overall F ST for the populations in Khon Kaen Province was 0.038 [95% confidence interval (CI) 0.002-0.105], and that of the more widely spaced populations in the previous study was 0.043 (CI = 0.016-0.075) [13]. Within Khon Kaen Province, (Chi River wetland genetic group) in this study, 67% and 33% of the populations had 'little' and 'moderate' genetic differentiation, respectively (Table 3). Of the five widely spaced populations studied [13], 60% had 'little' and 30% had 'moderate' genetic differentiation. 'Great' genetic differentiation (F ST = 0.165) occurred in only a single comparison, between LP and NP. Comparisons between Khon Kaen and widely spaced populations yielded 60% and 40% of populations with little and moderate differentiation, respectively. A Mantel regression test was done to determine the correlation between the F ST and geographical distance between populations. No correlation between genetic and geographic distance was found among populations from Khon Kaen as well as between Khon Kaen and widely spaced populations, whether distances were calculated in straight lines or along river courses (Table 3). At the broader geographic scale, no correlation between genetic and geographic distance was found when distances were measured along river courses, but was found when distances were measured in straight lines (Table 3). That is, isolation-by-distance was only indicated at the broader geographic scale and only when straight-line distances between populations were used.

Discussion
The population genetic data on O. viverrini from the four Khon Kaen localities considered in this study showed that there was considerable variation in allelic diversity, heterozygosity and allelic richness. Particularly, allelic richness for worms from KPv was significantly higher than for worms from the other three localities. The cause of this genetic diversity may be due to the transmission  dynamics of the parasite's life cycle as a consequence of selection against specific genotypes of parasite by different species of host (snails, fish or humans). Of these hosts, snails in particular show high levels of genetic diversity [25] and the potential for coevolution between parasite and host species, i.e. the Bithynia snail intermediate host, which in turn may play a role in the observed genetic, biological and/or morphologically variation in this parasite.
An alternate explanation of relatively low genetic diversity in three localities (KBs, KBp and KLp), as opposed to high diversity in KPv, may be due to the history of the parasite control program by chemotherapy. All of these areas are endemic for opisthorchiasis and there are records of praziquantel treatments [26,27,28], although no details on the frequency and coverage are available. It is possible, as for Schistosoma mansoni, that parasite genetic diversity is reduced after praziquantel treatment [29]. In this study the mean number of alleles per locus ranged from 2.83 to 3.7 which is lower than that found in other trematode parasites such as schistosomes [30,31].
O. viverrini sensu lato contains at least three species each genetically distinct with two morphologically similar (hence cryptic species) in river wetlands in Thailand and Lao PDR and one morphologically and biologically distinct isolate in Sakon Nakhon and Nakhon Phanom in Songkram River wetland, Thailand [9,10]. Additionally there are five distinct genetic groups of isolates that correspond to different wetlands which may in turn be different species within the complex as they have fixed genetic differences to an extent that define such species in many other parasite taxa [32].
Interestingly, even with the sample size of 30 adult worms per Khon Kaen locality analyzed in this study, a private or unique allele was observed in five of the 12 microsatellite loci we examined, albeit at low frequency (,6%). It is possible that with a larger sample size more unique alleles can be detected and could be used as markers to differentiate between populations of O. viverrini. In the case of other parasites, such as S. mansoni, a similar range of frequency of unique (private) alleles (1.1-4.1%) was observed [30].
Only a single locus (Ovms10) showed significant negative F IS values for all samples. This could be the result of possible extensive migration rates of the second intermediate host causing cross fertilization between distinct populations of O. viverrini.
Most pairs of populations presented in Table 3 are ''significantly'' differentiated from each other, indicating lack of panmixia. Actual values of F ST suggest mostly low to moderate differentiation between populations. Gene flow between the closely spaced populations analysed in Khon Kaen Province and the populations of O. viverrini from other wetlands in Thailand and Lao PDR is unlikely via river flow and flooding patterns of snails and fish hosts as each wetland is distinct. Additionally, the Nam Ngum River, which contains a genetically very distinct O. viverrini cryptic species, enters the Mekong River 578 km upstream from where the Mun River meets the Mekong River. Whether, human and/or food (in this case fish) cross border movement may provide an avenue for parasite gene flow in Thailand and Lao PDR requires further investigations.
The observed levels of genetic differentiation between the four spatially close populations of O. viverrini (10 to 60 km separation) within the Chi River wetland in Khon Kaen and more distant populations (up to 770 km separation) including Thailand and Lao PDR river wetland indicates the existence of intra-specific population structuring [13]. A significant correlation between genetic differentiation and geographical distance observed only for the distant populations (although only when distances were measured in straight lines) suggests that geographical separation is an important factor for population structure. However, connectivity along a river course between localities, and other factors such as drug treatment and the level of elevation, may influence the variation in genetic differentiation as hypothesized in schistosomes [29,31]. Thus, data from this and the previous study [13] suggest that O. viverrini is not a panmictic population but rather is differentiated genetically into different gene pools, indicating the existence of intraspecific population structures that may be associated with different cryptic species within defined wetlands that make up the O. viverrini species complex. Although the average pairwise F ST within a single species for the Khon Kaen populations in this study was not different from that between cryptic species in the widely spaced populations [13], the F ST estimated from the overall loci of the cryptic species populations (0.043) was higher than that for the closely spaced Khon Kaen populations (0.038). Substantial genetic differentiation (F ST = 0.165) was observed herein only between the cryptic species at a large geographical scale. This provides further independent evidence to support the hypothesis of the existence of cryptic species of O. viverrini within Thailand and Lao PDR [9] and is similar to the situation found for S. japonicum [31].
Several factors such as the intensity of O. viverrini infection, the life style and behavior, as well as genetic polymorphism in human genes and the frequency of praziquantel treatment may influence the transmission dynamics of O. viverrini and thus the incidence of CCA [4,33,34,35]. Khon Kaen Province is known to have a generally high prevalence of opisthorchiasis and incidence of CCA [36], which includes the same small geographical area examined here. Different levels of genetic diversity of O. viverrini have been found between localities in close proximity and even within the Lawa Lake, Ban Bhai district [11,37]. These results suggest that O. viverrini from different reservoirs and streams, and hence different ecological and geographical environments, have different levels of gene/allele frequency and/or number of genotypes, as has been shown for schistosomes [38,39,40,41,42].
Although expulsion chemotherapy to collect adult worms from humans is possible for O. viverrini [43,44], it is challenging under field conditions and has its limitations. For instance, worm recovery is unpredictable, most infected individuals have light infections and hence a low worm burden. In this study, caution in interpretation of the result is needed since only adult worms from experimental animals were used and these may not fully reflect the situation in the field due to potential laboratory host selection. Further study on patterns of genetic diversity and population genetic structure of O. viverrini using eggs from feces, cercariae from snails and metacercariae from fish should be compared to the adult worms to determine whether there is any affect of host selection or not as has been shown for schistosomes [45].
In conclusion we have shown that substantial genetic diversity and population genetic differentiation exists between four geographically close localities in the northeast of Thailand. Significantly higher allelic richness was found in worms from the KPv locality compared to worms from the other three localities. The overall genetic diversity and population structure observed at a relatively small spatial scale with a maximum population separation of 60 km was largely similar to that found at a much larger scale where the populations analyzed were separated by a distance of up to 770 km. The level of genetic differentiation was, however, significantly correlated with the distance between populations.