INTRODUCTION
Onchocerciasis, a disease commonly referred to as river blindness, is caused by the filarial worm Onchocerca volvulus. In 1974, the Onchocerciasis Control Program (OCP) was launched as a regional programme in 7 West African countries; this was the first large-scale control programme instituted against onchocerciasis. The initial strategy of the OCP was aerial application of insecticide to kill larvae of the vector in riverine breeding sites. It was thought that to interrupt transmission of the parasite, insecticiding was required for approximately 14 years, the estimated life-span of adult worms in humans. In the late 1980s, the pharmaceutical firm Merck and Co. Inc. donated the microfilaricidal drug, ivermectin (IVM), to the OCP, which included Ghana, and later to the African Program for Onchocerciasis Control (APOC). At the recommended single, yearly dose, IVM suppresses microfilariae (MF) in the skin and eyes, thus reducing the severe pathologies associated with the disease and parasite transmission (Green, Brown and Taylor, 1989). Mass IVM treatment continues in OCP and APOC areas, usually as the sole control measure.
IVM is a broad-spectrum anthelminthic that has had great success against gastrointestinal nematodes. IVM was introduced to the market in 1981 and in 1985 the first case of an IVM-resistant worm was found (Carmichael et al. 1987). Resistance to the avermectin-milbemycin class of drugs is now a major threat to the production of ruminants (Wolstenholme et al. 2004). Development of IVM resistance is influenced by many factors, including, but not limited to, the genetics of the worm, frequency of treatment and use of treatment strategies that restrict refugia (Wolstenholme et al. 2004).
Typically, annual IVM treatment is given to as many eligible members of a community, known to be endemic for O. volvulus, as possible. The population of O. volvulus larvae within the insect vector is small and short lived. The small proportion of the parasite population found in the vector and parasites in untreated members of the community comprise the population in refugia. Targeted treatment coverage is high (76%) and IVM effects are prolonged (in excess of 1 year). Thus, in communities with high annual treatment coverage, refugia may be small and IVM selection pressure high.
A variety of different tests are available to detect anthelminthic resistance in livestock, including post-mortem worm count reduction, fecal egg count reduction test, the egg hatch and larval development assays, the larval motility or paralysis test and, in some cases, PCR. With human onchocerciasis, it is difficult to prove unequivocally whether reduced efficacies (Ali et al. 2002; Awadzi et al. 2004a,b) were caused by drug resistance. The mass treatment of large human populations to break transmission of O. volvulus, and the tendency to increase the frequency of IVM treatment in onchocerciasis from once per year to up to two or more times per year, will increase selection pressure for resistance (Duke et al. 1990, 1992; Awadzi et al. 1999; Gardon et al. 2002). Experimental procedures used to confirm drug resistance in livestock parasites can not be performed on humans or are inappropriate for tissue-dwelling filarial nematodes, which cannot be cultured in animal hosts, and which produce viviparous larvae. In addition, in vitro tests can not readily be performed on O. volvulus as all stages are parasitic and can not be cultured for prolonged periods outside of the human definitive host or Simulium intermediate host. The situation is further complicated because IVM is not markedly macrofilaricidal and the most significant drug effect is its prolonged suppression of reproduction. Thus, alternative measures are required to assess possible resistance selection in O. volvulus. Drug resistance is genetic and will be reflected in changes to genetic parameters. Should resistance develop the allele(s) capable of conferring resistance need to be present in O. volvulus populations before treatment. The early stages of resistance would be manifested as a change in frequency and diversity of alleles linked to resistance. P-glycoprotein is a logical candidate gene because ivermectin is an excellent substrate for P-glycoprotein transport (Didier and Loor, 1996; Pouliot et al. 1997), it has already been implicated in IVM selection in O. volvulus (Ardelli, Guerriero and Prichard, 2005; Eng and Prichard, 2005) and in other IVM-resistant nematodes (Blackhall et al. 1998; Xu et al. 1998; Le Jambre, Lenane and Wardrop, 1999; Sangster et al. 1999; Kerbouef et al. 2003; Wolstenholme et al. 2004) and modifications in this protein could reduce the effective drug concentrations in the parasite (Juliano and Ling, 1976; Lincke et al. 1993). Ardelli et al. (2005) found changes in allelic patterns and a reduction in diversity at many loci in P-glycoprotein in O. volvulus from IVM-treated patients which suggested that IVM is imposing selection on this gene, consistent with a possible development of IVM resistance.
Theoretically, genes under selection can be revealed by unique patterns of linkage disequilibrium (LD) and polymorphism at physically linked loci. However, the effects recombination and mutation may have on the extent and persistence of LD patterns in natural populations is largely unknown. Besides changes in allele polymorphism and genetic diversity, it was of interest to undertake additional genetic analyses to determine whether there were any changes in linkage disequilibrium, the extent of heterozygosity, the proportion of markers in Hardy-Weinberg equilibrium and genotype frequency with increasing IVM treatment. It is also of interest because there are no studies which examine the effects that drug treatment may have on LD in parasitic nematodes. The 28 genetic markers, examined for polymorphism in the Ardelli et al. (2005) study, were analysed to determine whether repeated IVM treatment may lead to disruption of LD in P-glycoprotein, in O. volvulus.
MATERIALS AND METHODS
Origin of Onchocerca volvulus samples
The origins of the O. volvulus samples have been described by Ardelli et al. (2005). A total of 215 adult O. volvulus samples were used in this study. Of these, 79 (31 IVM-treated, 48 non-treated) were collected in the Volta Region of Ghana in 1999 and a further 136 (76 IVM-treated, 60 non-treated) were collected from the same region and villages (untreated: Ho District; treated: Kpando District) in 2002. The treated people had received between 4 and 10 annual IVM treatments. The procedure for isolating the worms has been described by (Ardelli and Prichard, 2004).
Genetic data collection and statistical analysis
Primer pairs, for amplification of 28 regions of a P-glycoprotein gene of O. volvulus (Huang and Prichard, 1999; GenBank Accession no. AY884212), were designed and used to amplify these regions of the gene (Table 1 and Fig. 1). DNA preparation from worms and marker typing (using SSCP or RFLP) followed procedures outlined by Ardelli and Prichard (2004). For each of the 28 regions RFLP and SSCP analyses were conducted and the results were reported by Ardelli et al. (2005). In the present investigation, for each region of the gene and for each population, genotype diversity, Hardy-Weinberg equilibrium (HWE), linkage disequilibrium (LD) and genotype frequency were determined, using the allelic frequencies as outlined in Fig. 1.
Fig. 1. Genetic diversity of the Onchocerca volvulus P-glycoprotein locus and location of 28 markers amplified from the gene. The dark grey boxes indicate transmembrane domains and black boxes indicate the Walker A and Walker B motifs within the nucleotide binding domain. O. volvulus P-glycoprotein mRNA is indicated at the top of the figure and lines, which indicate the regions that were examined, are above the mRNA and numbered 1–28. The frequency of the alleles at each region is shown with the region on the corresponding mRNA sequence indicated by an arrow and a line above the frequency histogram. The frequency histogram for each of the 4 groups is indicated on the figure. The x-axis shows the number of alleles for each region and are numbered consecutively (i.e. region 1 has alleles 1·1, 1·2, 1·3, etc.). The genotype which corresponds to the allele is indicated.
Table 1. Primer pairs for amplification of 28 regions of the Onchocerca volvulus P-glycoprotein locus (Note: Other PCR conditions were annealing time of 30 s, extension temperature of 72·0 °C and extension time of 30 s for 35 cycles.)
For analysis of genomic diversity within treated or non-treated populations, the unbiased estimate of expected genetic diversity (Ĥe) and mean unbiased gene diversity was applied, where hek is the value of he for kth region and r is the total number of regions (n=28) studied; he is an unbiased estimate of genetic diversity for each region and is determined by he=n(1−3x2i)/(n−1), where n is the number of homozygote individuals and x is the corresponding frequency of ith allele at a region in a sample from the population (Nei and Li, 1979; Nei, 1987; Kis-Papo et al. 2003). For each of the 28 regions a value for gene diversity was calculated. The value for gene diversity (see Fig. 2A) for each time-period (i.e. 1999 or 2002) and treatment level (i.e. no IVM, 4-IVM, etc.) is the mean value of gene diversity for all 28 regions (i.e. gene diversity for each region divided by the total number of regions). A parametric paired t-test was used to determine if there was a significant loss in gene diversity with increasing treatment.
Fig. 2. Association of IVM treatment with gene diversity (A), the fraction of markers that depart from HWE expectations (B), the percentage of heterozygotes (C) and fraction of markers in complete linkage disequilibrium (D). For 1999, 48 samples were obtained from non-treated patients and 31 samples were obtained from treated patients (4–6 IVM – n=7; 7 IVM – n=11; 8 IVM – n=13). For 2002, 60 samples were obtained from non-treated patients and 76 samples were obtained from treated patients (4–6 IVM – n=21; 7–8 IVM – n=43; 9–10 IVM – n=13).
HWE expectations for genotype frequencies, and their associated probabilities, were calculated using the program package Helix Tree (http:/www.goldenhelix.com). For each of the 28 regions of the P-glycoprotein locus, the genotype frequency of the alleles was tabulated before and following treatment with IVM. These frequencies were used to calculate the chi-squared distribution, using the formula χ2=Σki=1 (nii−np2i)2/np2i+Σki=1 Σkj=1, j≠1 (nij−2npipij)2/ 2npipj, where the frequency of alleles was expressed as pii=nii/n, where n is the population count, pi the frequency of homozygous alleles, pj the frequency of heterozygous alleles and k represents the degrees of freedom. The Hardy-Weinberg equilibrium coefficient was then defined as pij=2Dij (for i≠j), where D is the Hardy-Weinberg equilibrium coefficient. For each region, HWE was calculated before treatment with IVM and following treatment. A significant departure from HWE (χ2 0·05) following treatment indicated that the region was no longer in HWE.
Probabilities associated with coefficients of LD for all 28 regions and for each population were calculated using the Expectation/Maximization (EM) algorithm as implemented in Helix Tree. Using this program, the D′ contributions were defined as D′=kΣi=1mΣj=1piqj¦D′ij¦, where k and m were 2 markers with alleles l.....k and l.....m, respectively, which have frequencies pi....pk and ql....qm, respectively. The pairwise LD parameter (D′), which ranged from 0·1 to 0·9 was calculated for each of the 28 regions and for each group. This scale of D′ represents the P-values and confidence intervals determined from permutation tests and gives thresholds for the unadjusted P-values at P<0·001 (marginal significance) and P<0·00001 (highly significant).
The software program Haplotyper, which uses the Partition-Ligation algorithm to reconstruct individual genotypes from population genotype data (Niu et al. 2002), was used to estimate P-glycoprotein genotypes. All alleles for each of the 28 regions were used and the program was run for 20 rounds, as suggested for small sample sizes. The posterior probabilities associated with the selected pair of genotypes had to be greater than P=0·9 for the genotype to be used in the analysis of genotype frequencies.
Primers (5′-CAA TCA TGG GGA AGT CCA AAG-3′; 3′-CTC AAA ACC TTC CTT TGC AAT-5′) were used to amplify a region from O. volvulus heat-shock protein-60 (Hsp-60; GenBank Accession no. AF121264). This region of OvHsp-60 was examined in the 4 populations of O. volvulus collected in the Volta Region of Ghana. Gene diversity was determined and genotypes were constructed as described previously for the P-glycoprotein locus. OvHsp-60 was chosen as a control gene as it is not expected to be involved in IVM selection, but is polymorphic.
RESULTS
The 4 populations of O. volvulus sampled from the Volta Region came from untreated and IVM-treated patients (4–10 IVM treatments) collected in 1999 and 2002. Data were collected from 28 regions spanning 10574 base pairs of the O. volvulus P-glycoprotein gene which had a mean of 2·38 allelic polymorphs per region (range 2–10 – Fig. 1). Mean gene diversity (H) systematically decreased with IVM treatment (Fig. 2A) whereas the number of regions that significantly deviated from HWE increased with the number of IVM treatments (Fig. 2B). Gene diversity (as indicated by values of H) significantly decreased (P=0·009) between the 1999 and 2002 samples; samples collected in 1999 were more diverse than those collected in 2002. Interestingly, gene diversity was lower in non-treated samples (H=0·0167) collected in 2002, compared with non-treated samples collected from the same villages in 1999 (H=0·0231). Although correlated, the relationship between departure from HWE and treatment level was not statistically significant (Fig. 2B). Further analysis of any such relationship may require larger sample sizes. The data do suggest that the most significant effects of IVM are observed after at least 6 annual IVM treatments. Although there was an increase in the percentage of heterozygotes with treatment, these were still in relatively low proportions in comparison to the homozygotes (Fig. 2C).
Pairwise LD between loci was determined using ¦D′¦ (Fig. 3A–D). In the present study, the number of regions in LD decreased with treatment level and with time (Fig. 2D). The non-treated populations showed greater LD, in comparison to treated populations. However, LD decreased in both treated and non-treated populations between 1999 and 2002. Although the non-treated samples from 2002 showed more LD than the treated, it was much lower than that for the samples collected in 1999. However, 7 of the 28 regions examined remained in complete LD, despite the number of treatments and the year in which they were collected (Fig. 4A–D).
Fig. 3. Linkage disequilibrium (LD) among the 28 markers on the Onchocerca volvulus P-glycoprotein gene. The pairwise LD parameter (D′) was calculated for each of the 28 regions and these are numbered (i.e. 1–28) on the figure. A scale is shown at the bottom of the figure which represents the P-values and confidence intervals determined from permutation tests, giving thresholds for the unadjusted P-values at P<0·001 (marginal significance, light grey) and P<0·00001 (highly significant, white). Nine shades of grey are depicted on the scale which correspond to D′ ranging from 0·1 to 0·9. White boxes indicate regions in complete linkage disequilibrium (D′=1·0). Black boxes indicate regions in equilibrium (D′=0). More LD was found in non-treated samples collected in the Volta Region in 1999 (n=48) (A), compared with those collected from treated patients in 1999 (n=31) (B), from non-treated patients in 2002 (n=60) (C), or from treated patients in 2002 (n=76) (D). The regions are indicated on the graph to illustrate the pairwise comparisons.
Fig. 4. Schematic representation of Onchocerca volvulus P-glycoprotein mRNA sequence, showing regions predicted to be transmembrane helices (dark grey boxes – TMD) or nucleotide binding domains (black boxes – NBDs) as indicated by the solid horizontal bars. Each TMD is numbered according to its predicted structure. Within each NBD, the locations of the Walker A and Walker B motifs are indicated by uppercase letters. Only those regions which show complete LD are indicated by the boxes containing stripes. (A) No IVM, 1999; (B) IVM, 1999; (C) No IVM, 2002; (D) IVM, 2002.
Genotype frequencies for the 28 regions were derived from genotype data using a Partition-Ligation algorithm (Fig. 5A, B). Assuming random association among the SNPs within each region, a simulation utilizing 20 iterations (as suggested for the software package) based on sample size and observed allele frequencies, estimated 13 (No IVM, 1999), 11 (IVM, 1999), 15 (No IVM, 2002) and 17 (IVM, 2002) genotypes in the O. volvulus populations collected in the Volta Region of Ghana. An increase in the number of rare genotypes was noted in the IVM-treated worms from 2002 (Fig. 5A), although sample size was greater than in 1999.
Fig. 5. Genotype frequency estimation in 2 populations (1999 and 2002) of Onchocerca volvulus collected in the Volta Region of Ghana. Samples from non-treated patients were collected in the Ho District and samples from treated patients were collected from the Kpando district. (A) Genotype frequency from samples collected in 1999. (B) Genotype frequency from samples collected in 2002. Genotypes with the same number are the same for both the 1999 and 2002 data.
Gene diversity was examined in a control gene, OvHsp-60, from O. volvulus (Fig. 6A). Gene diversity in OvHsp-60 was not affected by treatment with IVM. Although gene diversity was lower in the 2002 samples, compared to the 1999 samples, this difference was not significant. Also, the values of H were higher in the control gene than in the OvPgp locus. The genotype frequency for OvHsp-60 was estimated for the 1999 and 2002 samples (Fig. 6B, C). In both the 1999 and 2002 samples, genotype A and H were in the highest frequency.
Fig. 6. Gene diversity and genotype frequency in Onchocerca volvulus heat shock protein-60 corresponding to region 214–590 on the cDNA sequence. (A) Gene diversity in 1999 and 2002 samples; (B) genotype frequency in 1999 samples; (C) genotype frequency in 2002 samples.
DISCUSSION
The results of the present study demonstrate a reduced genetic variation and disruption in LD in O. volvulus populations treated with IVM. This suggests that the widespread use of IVM in Ghana for onchocerciasis control is exerting selection pressure on O. volvulus and as a result is reducing gene diversity, even in worms removed from untreated patients. The recommended dosage of IVM is between 150 and 200 μg/kg and is usually administered as a single annual dose (Ali et al. 2002). At this dosage, dermal MF decrease rapidly with greater than 95% reduction within 1 week of treatment. After 1 year, skin MF densities increase, but the levels are approximately 10–20% of pre-treatment levels. IVM also affects adult O. volvulus fecundity and slows the release of MF from the uterus (Plaisier et al. 1995). It was shown that 3-monthly treatment with IVM for 3 years killed more adult female worms than standard annual treatment (Gardon et al. 2002). Considering the prolonged widespread use of IVM in Ghana and the effects IVM has on reproduction in O. volvulus, it is possible that patients, in the few communities not being treated with IVM, may become selectively infected with larvae that were progeny of worms better able to reproduce, despite IVM treatment. Alternatively, IVM may be selectively eliminating adult worms with specific P-glycoprotein alleles (e.g. susceptible), allowing only resistant worms to produce MF. Either or both of these mechanisms could account for a reduction in gene diversity in non-treated patients. Although the interval between the two collection periods (i.e. 1999 and 2002) is not long, it is sufficient for 3 generations of O. volvulus. If IVM kills certain female worms (Gardon et al. 2002), and affects fecundity in others, this might have an effect on the genotype of O. volvulus transmitted and be reflected in a loss of gene diversity.
The data suggest that the most significant genetic changes caused by IVM are observed after at least 6 annual IVM treatments. Although correlated, the relationship between departure from HWE and treatment level was not statistically significant. However, tests for HWE lack statistical power when HWE expectations are calculated from gene frequencies measured after selection has already occurred (Hartl and Clark, 1999). The dramatic loss of gene diversity observed in non-treated populations suggests that these populations of O. volvulus from the Volta Region have been affected by adjacent selection. Many of the villages that were sampled in this study are located in close proximity to one another. Blackflies are capable of travelling over 500 km, assisted by monsoon winds (Garms, Walsh and Davies, 1979). Thus, fly movements could potentially promote some dispersion of MF from IVM-treated worms into areas where there was no IVM distribution. This may account for the loss in gene diversity in non-treated populations. The results confirm and extend findings (Ardelli et al. 2005; Eng and Prichard, 2005) that IVM is affecting the P-glycoprotein locus, as there was a loss in diversity in this gene with treatment and between sampling periods. Eng and Prichard (2005) examined polymorphism in short sequences of 10 control genes from male O. volvulus collected in 1999, including Hsp-60, and found no evidence of IVM selection on these control genes. Eng and Prichard (2005) found O. volvulus Hsp-60 to be quite polymorphic and therefore a suitable gene to analyse for comparison with P-glycoprotein. This polymorphic control gene did not show a loss of gene diversity following treatment with IVM.
It is not known if IVM treatment of O. volvulus produces a dominant selection model. Measurement of the degree of dominance is generally determined using genetic crosses (Le Jambre, Royal and Martin, 1979; Herlich, Rew and Colglazier, 1981), but this approach is not possible in O. volvulus. IVM resistance in an isolate of the nematode Haemonchus contortus has been examined using genetic crosses (Dobson, Le Jambre and Gill, 1996). This study suggested a single gene or cluster of genes was responsible for a dominant resistance effect. In O. volvulus we found an increase in the percentage of heterozygotes occurred with treatment. However, heterozygotes were still in a relatively low proportion in the samples from the treated people, in comparison to the frequency of homozygotes. There was loss of both homozygotes and heterozygotes, between samples collected in 1999 and 2002, reflecting the decrease in gene diversity.
Intense natural selection is predicted to cause loss of LD between loci physically linked to genes under selection. However, the size of the genomic region and the persistence of LD over time are not stable in natural populations, given that they are subject to recombination and mutation. Parasite transmission rates in the local geographical regions (i.e. Ho and Kpando) were not known. IVM treatment reduces transmission by reducing the density of skin MF available for blackflies (Dadzie, Neira and Hopkins, 2003). However, density-dependent processes in O. volvulus (Basáñez et al. 1996) may amplify selection for drug resistance, inbreeding and reduce LD. While the minimal generation time for O. volvulus is about 1 year, given these conditions, meiotic recombination over a short number of years may not have much impact on LD for a genome segment of only a few kilobases and a loss in genetic diversity due to drug selection would be seen as a decrease in LD.
Several factors will determine the degree of LD, which can affect the chance of identifying an association between genetic markers and variant genes. The number of parasite generations since the original mutation is an important factor because it determines the number of opportunities for meiotic recombination to disrupt LD. The age of the mutations is also of importance as the older the mutation, the more subsequent generations a parasite could have passed through and the more chances of recombination. The transmission rate is also a crucial parameter in LD studies. A higher transmission rate will lead to more chances for meiotic chromosome assortment and crossovers. Passage of O. volvulus larval stages through blackflies is the only opportunity for parasites to enter a different host, and thus increase genetic diversity through sexual recombination. The degree to which sexual recombination occurs in O. volvulus is affected by the chance of a person being bitten by blackflies containing genetically distinct parasite larvae, as recombination will only be observed when infections are genetically heterogeneous. The allelic complexity of a local parasite population has an important impact on the detection of LD. The complexity is influenced by many factors, including inbreeding. Inbreeding would be expected to predominate in regions where infections containing several different genotypes are rare, when genetic diversity is low, or when transmission is low. Inbreeding could create local subpopulation structures that could affect LD. IVM treatment reduces the number of MF in the skin and, as a consequence, fewer parasites would be available for transmission.
The control practices implemented by the OCP have applied similar IVM usage (i.e. single, annual dose of 150–200 μg/kg), with targeted annual coverage in excess of 65% of the eligible population across much of Ghana for several years, which is likely to have imposed strong selection against drug-sensitive phenotypes. Since IVM treatment was only implemented as part of OCP in the late 1980s, O. volvulus populations in Ghana would have been exposed to a maximum of 10 years of IVM treatment for the 1999 samples and 13 years for the 2002 samples, which could include 10 or more generations of O. volvulus. The hallmarks of resistance selection, in particular a reduction in genetic diversity are apparent in O. volvulus removed from patients in the Ho and Kpando districts of the Volta Region of Ghana. Although the patterns of LD varied for each population, 7 markers remained in complete LD for both treatment populations. There is, as of yet, no unequivocally documented resistance to IVM in O. volvulus, although several recent reports have indicated suboptimal responses to IVM in some patients from the Pru and Lower Black Volta river basins of Ghana (Awadzi et al. 2004a,b), in the Khartoum region of Sudan (Ali et al. 2002), and higher levels of O. volvulus recrudescence in OCP areas under IVM treatment than expected (Soumbey-Alley et al. 2004). The spread of IVM resistance in livestock was rapid and was characterized initially by a loss of polymorphism at certain loci. With onchocerciasis, the greatest obstacle to studying the possible emergence of IVM resistance is the lack of tools to detect and recognize it, should it arise. The results of this study demonstrate a dramatic loss of diversity in P-glycoprotein, even in non-treated populations in a region where extensive treatment has occurred over many years, and a disruption in LD. This should serve as an early warning that mass use of IVM is selecting against certain alleles of P-glycoprotein in populations of O. volvulus.
This investigation received financial support from the African Programme for Onchocerciasis Control, the Onchocerciasis Control Program in West Africa and the UNDP/World Bank/WHO Special Programme for Research and Training in Tropical Diseases. B.F.A. was supported by fellowships from the Natural Sciences and Engineering Research Council (NSERC) of Canada and bridging funds from a Fonds Québécois de la recherche sur la nature et les technologies (FQRNT) Centre Grant (Centre for Host-Parasite Interactions). S.B.G. was supported by an NSERC undergraduate research award. We gratefully acknowledge Dr K. Awadzi for collecting Onchocerca volvulus nodules from Ghana in 1999 and Mr M. Osei-Atweneboana for collecting the nodules in Ghana in 2002.
