INTRODUCTION
Fish-borne trematodiasis is becoming increasingly recognized as a serious public health problem in Southeast Asia. Three species within the family Opisthorchidae most commonly implicated as aetiological agents are Opisthorchis viverrini, Opisthorchis felineus and Clonorchis sinensis. These parasites infect at least 40 million people within a disease risk region estimated to include 700 million people (WHO, 1995; Keiser and Utzinger, Reference Keiser and Utzinger2005).
Of these three species of opisthorchid trematode, only O. viverrini is classified as a type 1 carcinogen because of its role as an initiator of chronic inflammation and the subsequent development of cholangiocarcinoma (CCA) (IARC, 1994; Sithithaworn et al. Reference Sithithaworn, Haswell-Elkins, Mairiang, Satarug, Mairiang, Vatanasapt and Elkins1994; Honjo et al. Reference Honjo, Srivatanakul, Sriplung, Kikukawa, Hanai, Uchida, Todoroki, Jedpiyawongse, Kittiwatanachot, Sripa, Deerasamee and Miwa2005). This liver fluke is most commonly found in Southeast Asia (WHO, 1995). An estimated 9 million people are infected in Thailand with the region of highest endemicity being in the north and northeast (Sithithaworn and Haswell-Elkins, Reference Sithithaworn and Haswell-Elkins2003). There is considerable variation, however, in parasite prevalence and disease presentation in different geographical areas (Sithithaworn and Haswell-Elkins, Reference Sithithaworn and Haswell-Elkins2003; Sriamporn et al. Reference Sriamporn, Pisani, Pipitgool, Suwanrungruang, Kamsa-ard and Parkin2004), the latter of which may be associated with genetic differences among parasites (Saijuntha et al. Reference Saijuntha, Sithithaworn, Wongkham, Laha, Pipitgool, Tesana, Chilton, Petney and Andrews2007).
Currently, there is inadequate information on the level of genetic variation within and among populations of O. viverrini throughout its geographical distribution. Previous studies (Ando et al. Reference Ando, Sithithaworn, Nuchjungreed, Tesana, Srisawangwong, Limviroj and Chinzei2001; Saijuntha et al. Reference Saijuntha, Sithithaworn, Wongkham, Laha, Pipitgool, Petney, Chilton and Andrews2006a, Reference Saijuntha, Sithithaworn, Wongkham, Laha, Pipitgool, Petney and Andrewsb; Reference Saijuntha, Sithithaworn, Wongkham, Laha, Pipitgool, Tesana, Chilton, Petney and Andrews2007; Sithithaworn et al. Reference Sithithaworn, Nuchjungreed, Srisawangwong, Ando, Petney, Chilton and Andrews2007) have reported genetic variation among isolates of O. viverrini from Thailand and the Lao PDR. These studies have been based on sequence analyses of the mitochondrial (mt) DNA cytochrome c oxidase subunit 1 gene (cox1) (Ando et al. Reference Ando, Sithithaworn, Nuchjungreed, Tesana, Srisawangwong, Limviroj and Chinzei2001), RAPD analyses (Sithithaworn et al. Reference Sithithaworn, Nuchjungreed, Srisawangwong, Ando, Petney, Chilton and Andrews2007) and multilocus enzyme electrophoresis analyses (Saijuntha et al. Reference Saijuntha, Sithithaworn, Wongkham, Laha, Pipitgool, Petney, Chilton and Andrews2006a, Reference Saijuntha, Sithithaworn, Wongkham, Laha, Pipitgool, Petney and Andrewsb, Reference Saijuntha, Sithithaworn, Wongkham, Laha, Pipitgool, Tesana, Chilton, Petney and Andrews2007). The findings of these studies have suggested that there is population substructuring among O. viverrini populations from different geographical areas. Furthermore, Saijuntha et al. (Reference Saijuntha, Sithithaworn, Wongkham, Laha, Pipitgool, Tesana, Chilton, Petney and Andrews2007) demonstrated that O. viverrini represents species complex of at least 2 species. Nonetheless, more DNA-based markers need to be established to determine the magnitude of genetic variation among populations of this liver fluke throughout different parts of its distributional range.
Mitochondrial DNA sequence data have been used widely to examine the population genetic structures of animals, including parasitic platyhelminths (e.g. McManus and Bowles, Reference McManus and Bowles1996; Le et al. Reference Le, Blair and McManus2000), because of the apparent maternal mode of inheritance and high mutation rates. Such sequences have proven useful for the analysis of within-species variation in parasitic flatworms, such as Fasciola hepatica (see Semyenova et al. Reference Semyenova, Morozova, Chrisanfova, Gorokhov, Arkhipov, Moskvin, Movsessyan and Ryskov2006), Paragonimus westermani (Park et al. Reference Park, Im and Yong2003), Clonorchis sinensis (see Lee and Huh, Reference Lee and Huh2004; Park, Reference Park2007) and Echinococus granulosus (see Bowles et al. Reference Bowles, Blair and McManus1992, Reference Bowles, Blair and McManus1994; Maillard et al. Reference Maillard, Benchikh-Elfegoun, Knapp, Bart, Koskei, Gottstein and Piarroux2007). For example, cox1 and the NADH dehydrogenase subunit 1 gene (nad1) have been shown to be useful population genetic markers for Schistosoma japonicum (see Bowles et al. Reference Bowles, Hope, Tiu, Liu and McManus1993; Sorensen et al. Reference Sorensen, Drew, Brindley, Bogh, Gasser, Qian, Chiping and McManus1998). Bowles et al. (Reference Bowles, Hope, Tiu, Liu and McManus1993) examined the degree of genetic divergence in cox1 between S. japonicum isolates from China and the Philippines using sequence data of cox1, whereas Sorensen et al. (Reference Sorensen, Drew, Brindley, Bogh, Gasser, Qian, Chiping and McManus1998) compared the level of genetic variation among populations of S. japonicum from 6 different geographical regions in mainland China. Morgan and Blair (Reference Morgan and Blair1998) have also shown the utility of the nad1 gene sequences for species identification and for comparing the rates of genetic divergence strains of echinostomes. Thus, cox1 and nad1 are of potential use to examine genetic variation among geographical isolates of O. viverrini.
In the present study, we investigated the level of genetic variation among isolates of O. viverrini from 14 different geographical areas in Thailand and the Lao PDR using cox1 and nad1 sequence data. In addition, we determined the phylogenetic relationships of O. viverrini to other trematodes based on analyses of sequence data of the cox1 gene.
MATERIALS AND METHODS
Collection of adult O. viverrini
Metacercariae of O. viverrini were obtained from infected cyprinid fish from 14 different geographical localities within Thailand and the Lao PDR (Fig. 1 and Table 1). Fish were mixed and digested with 0·3% pepsin solution. Metacercariae were identified using a light microscope and 50–100 metacercariae were used to infect hamsters. After 4–6 months, adult worms were collected from the bile duct, identified using morphological characters, then washed extensively in physiological saline, placed in pools (⩾20 worms) into microcentrifuge tubes, frozen live and stored at −80°C until used for the molecular analyses.
Localities in Thailand and the Lao PDR where Opisthorchis viverrini were collected (see Table 1 for the details and abbreviations of each locality).
Details of the source of the 14 isolates of Opisthorchis viverrini used in this study
Molecular methods
Pools of adult worms (⩾20 worms) were crushed using a tissue grinder. Lysis buffer and proteinase K (200 μg/ml) were added to each homogenized sample. Genomic DNA (gDNA) was extracted using the phenol/chloroform method described by Sambrook et al. (Reference Sambrook, Fritsch and Maniatis1989). Part of the nad1 gene was amplified by polymerase chain reaction (PCR) using the forward primer MNDI-A (5′-TAC GCA GGT GGT TTG GTT GG-3′) and the reverse primer MNDI-B (5′-CCC AAA GCT CAC ATC CTT GT-3′). These primers were designed based on a nad1 sequence for O. viverrini (GenBank Accession number DQ119551). Part of the cox1 gene was amplified by PCR using the forward primer MCOI-A (5′-TTT TTT GGG CAT CCT GAG GTT TA-3′) and the reverse primer MCOI-B (5′-TAA AGA AAG AAC ATA ATG AAA ATG AGC-3′) (Ando et al. Reference Ando, Sithithaworn, Nuchjungreed, Tesana, Srisawangwong, Limviroj and Chinzei2001). PCR reactions were conducted in a total volume of 25 μl containing 25 mm MgCl2, 10×buffer, 2·5 mm dNTPs, 25 pmol of each primer, 0·15 U of Go Taq DNA polymerase (Promega) and 100–300 ng of gDNA. No gDNA (i.e. negative) controls were included in each run. PCR reactions were performed using a MyCycler Thermal cycler (Bio-Rad Laboratories, Hercules, CA, USA). The cycling conditions used were: an initial denaturing cycle of 5 min at 94°C, then 30 cycles of 94°C for 1 min, 52°C (for nad1) or 50°C (for cox1) for 1 min and 72°C for 2 min. Amplicons were compared electrophoretically on 1·5% agarose gels (Sambrook et al. Reference Sambrook, Fritsch and Maniatis1989) stained with ethidium bromide and visualized using UV illumination. Specific bands were excised from agarose gels, purified using the UltraClean PCR kit (Mo-Bio Laboratories) and subjected to automated DNA sequencing using the same primers (in separate reactions) as employed for the primary PCR.
Data analysis
The nucleotide sequences of the nad1 and cox1 for O. viverrini obtained during this study have been deposited in the GenBank databases (Accession numbers EU022337-EU022350 and EU022351-EU022364, respectively). Nucleotide sequences were translated into amino acid sequence using Translation Tool 1.3 (http://www.bioinformatics.vg/bioinformatics_tools/tranlatetool.shtml), employing the codon table of invertebrate mitochondrial code. BioEdit version 5.0.6 was used for both nucleotide and amino acid sequence alignments. In addition, included in the nad1 alignment were sequence data for 5 isolates of O. viverrini retrieved from GenBank: Isolate OVN3 (Accession number DQ119551) from Khon Kaen, Thailand; OVBD (Accession number DQ882172) from Binh Dinh, Vietnam; isolate OVDL3 (Accession number DQ882174) from Dak Lak, Vietnam; isolate OVPY3 (Accession number DQ882173) from Phu Yen, Vietnam; and isolate OVL (Accession number DQ882175) from Vientiane, Laos. Two species of Taenia (Accession numbers NC_004826 and AF338826) and Ascaris suum (Accession number NC_001327) were included as outgroups. Also included in the cox1 alignment were publicly available sequence data for other trematodes, including O. felineus (Accession numbers DQ469316 and DQ469317), Clonorchis sinensis (Accession numbers AF181889 and AF188122), 2 species of Metagonimus (Accession numbers AF096231 and AF096230), 8 species of Schistosoma (Accession numbers NC_002544, NC_002529, AF295106, AF101196, NC_008074, DQ157223, AJ519521 and AY157201), 3 species of Echinostoma (Accession numbers AF025829, AF025824 and AF025823), 2 species of Fasciola (Accession numbers AB207180 and AF216697), Fasciolopsis buski (Accession number EF027094) and 7 species of Paragonimus (Accession numbers NC_002354, AY618842, AF538944, AY618838, AY618834, AY618806 and AF159599). Sequences of cox1 for the nematode A. suum (Accession number NC_001327) and 5 species of taeniid cestode (Accession numbers AF297617, AB018440, NC_004826, AB107245 and AB271234) were included as outgroups for the phylogenetic analyses. Phylogenetic trees were constructed using the PHYLIP program version 3.6.6. A distance matrix was calculated using the Kimura two-parameter model and the tree constructed using the neighbor-joining method (Saitou and Nei, Reference Saitou and Nei1987). The relative support for the clades in the NJ analyses was determined using 1000 bootstrap replicates.
RESULTS
The partial nad1 nucleotide sequences derived from all 14 geographical isolates were compared over an alignment length of 668 bp. There were only 2 alignment positions (nos. 434 and 470 Table 2) at which nucleotide variation was detected among the isolates included in the present study. However, comparison with previously published sequences revealed sequence variation at 9 alignment positions (Table 2). Four different haplotypes (designated as OVN1-OVN4) were detected among the 14 isolates of the present study. Haplotype OVN1 was detected in 6 isolates of O. viverrini, 2 from Thailand (CP and KPv) and 4 from Laos (NG, TH, VT and SV). Haplotype OVN2 was detected in 4 isolates from Thailand (KS, BR, NP and LP). Haplotype OVN3 was detected in 2 isolates from Thailand (KLp and SK), whereas haplotype OVN4 was found in another 2 Thai isolates (KBp and MS) (Table 2). One of the publicly available sequences (Accession number DQ119551) was identical in sequence to that of haplotype OVN3. The 4 other O. viverrini sequences available publicly represented a unique haplotype compared with the sequences of the isolates examined in the present study (Table 2). The nad1 nucleotide sequences of 8 different haplotypes were conceptually translated to peptide sequences (221 amino acids in length). An alignment of these amino acid sequences revealed differences at 4 positions (Table 2). The 8 nucleotide sequences translated to 5 different amino acid sequences.
Nucleotide and amino acid sequence differences in the nad1 and cox1 among haplotypes representing different isolates of Opisthorchis viverrini
1 Genotypes OVN1=isolated from CP, KPv, NG, TH, VT and SV; OVN2=KS, BR, NP and LP; OVN3=KLp, SK and DQ119551; OVN4=KBp and MS; OVBD=DQ882172; OVDL3=DQ882174; OVL=DQ882175; OVPY3=DQ882173; OVC1=KBp, KPv, CP, SK, NP, BR, LP, NG, TH and SV; OVC2=KLp, KS, MS and VT.
2 Position of DNA sequence differences.
3 Position of amino acid sequence difference.
A phylogenetic analysis of nucleotide sequence data revealed that haplotype OVN1 was genetically more similar to OVN4, whereas haplotype OVN2 was more similar to OVN3. These 4 groups were more closely related to each other than 3 isolates from Vietnam and 1 isolate from Laos (Fig. 2).
Phylogenetic tree depicting relationships between Opisthorchis viverrini isolates inferred from a neighbor-joining method of analysis of the nad1 nucleotide sequence data using two species of Taenia and Ascaris suum as outgroups. Bootstrap values (>50%) are indicated above branches.
The cox1 sequences from all 14 isolates of O. viverrini were compared over an alignment length of 394 bp. Only 2 sequence variants were detected (designated as haplotypes OVC1 and OVC2) and these differed at 2 alignment positions (Table 2). The mutational change at alignment position 370 resulted in a change in amino acid sequence (Table 2). Haplotype OVC1 was detected in 7 isolates from Thailand (designated KBp, KPv, CP, SK, NP, BR and LP) and 3 isolates from Laos (NG, TH and SV), whereas haplotype OVC2 was detected in 3 isolates from Thailand (designated KLp, KS and MS) and 1 isolate from Laos (VT).
A phylogenetic analysis of the 2 cox1 haplotypes of O. viverrini and other species of trematode (Fig. 3) revealed strong statistical support (i.e. bootstrap value of 91%) for a clade comprising sequences from O. viverrini, O. felineus and C. sinensis (members of the family Opisthorchidae). There was also some support (i.e. bootstrap value of 61%) for a clade that included 2 species of the genus Metagonimus with the opisthorchid nematodes. There were 3 other clades with varying levels of statistical support. The second clade contained all 7 species of Paragonimus (bootstrap value of 96%) that differed in sequence by 9–19%. The third clade, with moderate bootstrap support (71%) comprised 3 genera Echinostoma, Fasciola and Fasciolopsis, with sequence divergence among species ranging from 6 to 25%. The fourth clade with limited bootstrap support (65%) comprised the 8 species of Schistosoma, which differed in sequence by 7–23%.
Phylogenetic relationships of Opisthorchis viverrini with the other parasitic trematodes based on a neighbor-joining analysis of cox1 nucleotide sequence data. Ascaris suum (Nematoda) was used as the outgroup in the analyses. Bootstrap values (>50%) are indicated above branches.
DISCUSSION
Although genetic variation within O. viverrini sequences of the mtDNA gene cox1 has been reported previously (Ando et al. Reference Ando, Sithithaworn, Nuchjungreed, Tesana, Srisawangwong, Limviroj and Chinzei2001), this is a first study that has used sequence data of the nad1 gene to investigate genetic variation within this parasite species. In the present study we detected limited genetic variation among 14 geographical isolates of O. viverrini for the parts of the cox1 (394 bp) and nad1 (668 bp) genes. For both cox1 and nad1, only 2 nucleotide positions (0·5% and 0·3% respectively) were variable among O. viverrini isolates. However, in both genes, there were changes in amino acid sequences among isolates. A comparison of the nad1 sequences with those for 4 O. viverrini isolates deposited in the GenBank database revealed additional variable positions in the sequence alignment. Thus, at the nucleotide level there are at least 8 nad1 haplotypes for O. viverrini, representing 4 different amino acid sequence types.
Surprisingly, isolates separated by small geographical distances, did not necessarily have the same nad1 haplotype. For example, for the 3 isolates from Khon Kaen province, KPv had the haplotype OVN1, KLp the haplotype OVN3 and KBp the haplotype OVN4. This finding does not agree with the results of a multilocus enzyme electrophoretic study of isolates from the same geographical areas using 32 enzyme loci (Saijuntha et al. Reference Saijuntha, Sithithaworn, Wongkham, Laha, Pipitgool, Tesana, Chilton, Petney and Andrews2007). These authors found that O. viverrini isolates separated by small geographical distances and within the same defined wetland formed a distinct cluster. This difference may be a reflection of the different number of genes examined. The addition of sequence data for more, independent genes may indeed provide more reliable results and hence biological interpretations regarding the level of genetic variation detected among isolates O. viverrini examined in this study and their phylogenetic relationships to other trematodes.
A previous study (Ando et al. Reference Ando, Sithithaworn, Nuchjungreed, Tesana, Srisawangwong, Limviroj and Chinzei2001) using cox1 sequences compared 5 O. viverrini isolates from Ubon Ratana, Lerngpleuy, Ban Phai, Mahasarakham, and Chatturat. Each isolate had a different cox1 sequence. The cox1 sequences for these isolates could not be incorporated in the present study because the data were not available in any public database. In the present study of 14 isolates, 2 different haplotypes were detected, and these differed at 2 nucleotide positions. A possible explanation for the differences between studies may relate to the reproductive turn-over rate in O. viverrini populations, which occurs seasonally, as shown in the different rate of infection in cyprinid fish in each month (Sithithaworn et al. Reference Sithithaworn, Pipitgool, Srisawangwong, Elkins and Haswell-Elkins1997). Alternatively, preferential selection could result in lower numbers of haplotypes of O. viverrini being detected between the two studies.
The results of our study show that sequence data of cox1 can be used to establish evolutionary relationships of parasitic trematodes at the genus and species levels. Campos et al. (Reference Campos, Cummings, Reyes and Laclette1998) used the sequence data of the 18S nuclear ribosomal DNA gene to infer the phylogenetic relationships of taxa within the phylum Platyhelminthes. They found a clade comprising Fasciola, Fasciolopsis and Echinostoma, and that Opisthorchis was not closely related to the other parasitic trematodes included in our study. The phylogenetic tree derived using sequence data of the cox1 showed similar relationships. Furthermore, each cluster of different genera of parasitic trematodes was separated into clades that were correlated with the type of second intermediate host they utilize and their mode of transmission to the definitive host. Of the three clades of food-borne trematodes, the first (the genera Opisthorchis, Clonorchis and Metagonimus; order Opisthorchiida) use cyprinid fish as their second intermediate host, the second (Paragonimus spp.; Order Plagiorchiida) parasitize crabs, and the third (the genera Fasciola, Fascolopsis and Echinostoma; order Echinostomida) lack a second intermediate host (i.e. Fasciola, and Fascolopsis) and the definitive host consumes the metacercariae attached to water plants, or they use snails, mussels and frogs as second intermediate hosts (i.e. Echinostoma spp.). The fourth clade represents the waterborne trematodes (Schistosoma spp.; order Strigeatida) which lack a second intermediate host in their life cycle. The definitive host is infected by the ceracariae that directly penetrate the skin. This mode of transmission may be due to co-adaptation between clades of parasites and specific species of intermediate hosts. For example, it has been suggested by Saijuntha et al. (Reference Saijuntha, Sithithaworn, Wongkham, Laha, Pipitgool, Tesana, Chilton, Petney and Andrews2007) that O. viverrini may co-evolve with their first intermediate hosts, snails of the genus Bithynia. They found congruence between the phenogram of the genetic variation in O. viverrini from different geographical areas and the phenogram of their gastropod intermediate hosts based on analyses using multiple enzyme loci.
In conclusion, the number of different nad1 haplotypes of O. viverrini was greater than that previously reported for this parasite based on cox1 sequences. Hence, nad1 may prove useful to examine in greater detail the magnitude of population variation within O. viverrini, whereas the cox1 gene is more useful for examining of the phylogenetic relationships of food-borne trematodes.
This research was supported by a Wellcome Trust Collaborative Research Initiative Grant, the Thailand Research Fund through the Royal Golden Jubilee Ph.D. Program (grant no. PHD/0221/2544 to W. S.) and the Faculty of Medicine, Khon Kaen University Overseas Visiting Professor Program. Dr Petney's travel was supported by Deutsche Forschungsgemeinschaft grant (PE 1611/1-1).
