Temporal stability of Glossina fuscipes fuscipes populations in Uganda

Background Glossina fuscipes, a riverine species of tsetse, is the major vector of human African trypanosomiasis (HAT) in sub-Saharan Africa. Understanding the population dynamics, and specifically the temporal stability, of G. fuscipes will be important for informing vector control activities. We evaluated genetic changes over time in seven populations of the subspecies G. f. fuscipes distributed across southeastern Uganda, including a zone of contact between two historically isolated lineages. A total of 667 tsetse flies were genotyped at 16 microsatellite loci and at one mitochondrial locus. Results Results of an AMOVA indicated that time of sampling did not explain a significant proportion of the variance in allele frequencies observed across all samples. Estimates of differentiation between samples from a single population ranged from approximately 0 to 0.019, using Jost's DEST. Effective population size estimates using momentum-based and likelihood methods were generally large. We observed significant change in mitochondrial haplotype frequencies in just one population, located along the zone of contact. The change in haplotypes was not accompanied by changes in microsatellite frequencies, raising the possibility of asymmetric mating compatibility in this zone. Conclusion Our results suggest that populations of G. f. fuscipes were stable over the 8-12 generations studied. Future studies should aim to reconcile these data with observed seasonal fluctuations in the apparent density of tsetse.


Introduction
Tsetse flies, Glossina spp (Diptera: Glossinidae) transmit several species of pathogenic trypanosomes causing Human African Trypanosomiasis (HAT) and African Animal Trypanosomiasis (AAT). HAT affects human welfare directly through the chronic and acute forms of the disease caused by Trypanosoma brucei gambiense and T. b. rhodesiense respectively. AAT, on the other hand, stands as a major obstacle to the development of more efficient and sustainable livestock production systems in tsetse-infested areas [1]. A major challenge to controlling HAT is lack of suitable prophylactic drugs and vaccines against trypanosomiasis. Furthermore, chemotherapeutic agents for treatment of HAT are expensive, difficult to administer in remote areas and exhibit poor safety profiles. Consequently, vector control remains a viable alternative for large-scale control of trypanosomiasis.
Understanding tsetse population dynamics is critical for determining which control strategy is most appropriate (e.g., suppression, eradication), for choosing the best method for enacting that strategy (e.g., traps, insecticide-treated cattle, sterile insect technique), and for determining the scale at which vector control activities must be implemented [2]. Determinants of population dynamics include both life history and ecological correlates such as mating system, dispersal ability and population size, which influence the extent to which tsetse populations can recover from refugia following intervention, or re-colonize a cleared zone from neighboring sources. Recently, the use of population genetics has provided insights into tsetse ecology [3], with important ramifications for the implementation of control programs [4]. For example, studies of tsetse in Guinea and Senegal have identified populations that are sufficiently isolated to warrant attempts at complete elimination [5][6][7]. Elsewhere though, studies have documented relatively high levels of gene flow, necessitating integration of barriers into elimination schemes [8] or warranting an area-wide control effort that encompasses the dispersal-linked populations [9,10].
Across Africa, Glossina fuscipes is one of the most important vectors of HAT, transmitting an estimated 90% of all disease cases [11]. Glossina fuscipes is a member of the palpalis group of tsetse, which inhabit low bushes or forests at the margins of rivers, lakes or temporarily-flooded scrub land. In eastern Africa, populations of the subspecies G. f. fuscipes appear to respond to seasonal weather patterns, often disappearing during the bi-annual dry season from sites where they were previously abundant [12]. If populations in refugia are small, then seasonal bottlenecks could result in large temporal changes in gene frequencies.
In order to investigate the impact of seasonal climate changes on population size and to gain further insight into the population dynamics of G. f. fuscipes, we evaluated temporal changes in gene frequencies at one mitochondrial locus and 16 microsatellite loci in multiple Ugandan populations. Our sampling scheme included three populations situated at a zone of contact between two divergent lineages of G. f. fuscipes. These two lineages exhibit distinct mitochondrial DNA (mtDNA) haplotypes and strong differentiation at microsatellite loci, suggesting a long history of isolation, and providing a unique opportunity to monitor their interaction over time [9,10].

Tsetse collection and study area
Tsetse flies were collected using biconical traps [13] during the period from March 2008 to January 2010. All sites were sampled in 2008 [10] and then at least one year later in 2009 or 2010. Four sites were also sampled a third time (Table 1). Each fly was stored individually in 80% ethanol.
Tsetse collections were conducted at seven sites spanning central and southeastern Uganda (Figure 1). These sites generally reflected the riverine/woodland habitat preferred by G. fuscipes, but varied somewhat in regard to the immediate environment. Sites at Busime (BU) and Junda (JN) were located in a transition zone between marsh and woodland on the edge of Lake Victoria and Lake Kyoga, respectively. Sites at Bunghazi (BN), Dokolo/Otuboi (OT) and Okame (OK) were situated along permanent streams in a region of mixed agriculture and pastureland. Sampling at Mukongoro (MK) was conducted at the margin of ephemeral wetlands associated with rice cultivation. Sampling at Masindi (MS) was conducted within a region of banana and sugar cane plantations.
The study sites spanned a zone of contact between two divergent groups of G. fuscipes co-occuring in Uganda [9,10]. Sites MK and OT were situated north of the zone of contact and flies here were expected to possess solely northern mtDNA haplotypes. Sites BN, JN, and MS were located at the zone of contact and flies here were expected to possess both northern and southern haplotypes. Sites BU and OK were located south of the zone of contact and flies at these sites were expected to possess exclusively southern mtDNA haplotypes.

DNA Extraction
DNA was extracted from tsetse legs using NucleoSpin 96 Tissue Kits (Clontech, Mountain View, CA) or DNeasy kits (Qiagen, Valencia, CA) following the manufacturer's protocols.

Mitochondrial DNA sequencing
PCR was used to amplify a 570 bp fragment of mtDNA from a random subset of flies from each population using the primers COIF1 (CCT CAA CAC TTT TTA GGT TTA G) and COIIR1 (GGT TCT CTA ATT TCA TCA AGT A) as described by [10]. We amplified COIF1/COIIR1 in a 25 μl PCR reaction containing 1 × buffer (GoTaq colorless, Promega), 0.8 mM each dNTP, 0.4 mM primers, 1.5 mM MgCl 2 and 0.5 U Go Taq polymerase. The amplification involved a denaturation step at 95°C for 8 min, followed by 50 cycles each at 94°C for 30 s, 51°C for 30 s, 72°C for 45 s, with a final extension step at 72°C for 7 min. PCR products were sequenced using an ABI Model 3730 automated sequencer (Applied Biosystems, Foster City, CA, USA).  Table 1. The dotted line indicates the approximate extent of a zone of contact between two historically isolated groups of tsetse.

MS
Electropherograms were visually inspected and sequences were trimmed to remove poor quality data. The resulting sequences (530 bp) were aligned by eye using the computer program Sequencher 4.2.2 (Gene Codes Corporation).

Microsatellite genotyping
We genotyped individual flies at 16 loci. We used 11 of the 13 loci described by [10], excluding D05 and Pgp17 due to possible null allele problems. We also employed five new dinucleotide loci identified in the G. morsitans genome and optimized for use in G. fuscipes: GmmA06, GmmB20, GmmD15, GmmL03, GmmL11 [14]. Amplifications were performed with fluorescently labeled forward primers (6-FAM, HEX and NED) using a touchdown PCR (10 cycles of annealing at progressively lower temperatures from 60°C to 51°C followed by 35 cycles at 50°C) in 12.5 μl reaction volumes employing 1 × buffer, 0.8 mM dNTPs, 2.0 mM MgCl 2 and 0.5 U Go Taq polymerase. PCR products were multiplexed in groups of two or three and then genotyped on the ABI 3730 automated sequencer. Alleles were scored using the program Genemarker v 5.0 (SoftGenetics) with manual editing of the automatically scored peaks.

Marker validation and genetic diversity
Microsatellite loci were evaluated for Hardy Weinberg equilibrium (HWE) and linkage disequilibrium (LD) using Genepop version 4.0 [15]. Markov chain parameters were set at 10,000 dememorizations, 1000 batches, 10,000 iterations per batch for HWE and 100,000 dememorizations,1000 batches,10,000 iterations per batch for LD. We used the method of [16] as implemented in MultiTest v.1.2 to correct for multiple tests. Locus-and population-specific estimates of microsatellite allele frequencies were generated using the program Genalex version 6.3 [17]. We used the program FSTAT version 3.1 [18] to calculate allelic richness and the program Arlequin v. 3.5 [19] to calculate observed (H o ) and expected (H e ) heterozygosity for each population. DnaSP version 5.0 [20] was used to calculate mtDNA haplotype diversity (H d ) and nucleotide diversity (π).

Temporal genetic differentiation and population stability
For microsatellite data, we used Jost's D EST [21] to quantify genetic differentiation between populations and between temporal samples from the same population. D EST provides a less-biased estimate of differentiation than F ST and related statistics, especially when estimated using highly polymorphic microsatellite loci [22]. Locusspecific calculations of D EST were performed using the web-based program SMOGD [23] and then averaged across loci. For mtDNA data, we used Fisher's exact test and the statistical software SAS version 9.1 to test for differences in haplotype frequencies among temporal samples from the same population. For both microsatellite and mtDNA data, we performed an analysis of molecular variance (AMOVA) as implemented in Arlequin v. 3.5 [19] to characterize the proportion of the variance in microsatellite allele frequencies or haplotype frequencies that was attributable to differences in date of sampling. For this analysis, we used only the two samples from each population that were separated by the longest time interval.
We estimated current effective population sizes based on temporal changes in microsatellite allele at all seven sites. The effective size of a population (N e ) is defined as the size of an ideal population (i.e., one of constant size, discrete generations, and negligible selection and gene flow), which would exhibit the same genetic characteristics as the population at hand [24]. N e , therefore, reflects the rate of change in gene frequencies due to random drift alone [25]. We used two methods: the moment-based approach [26] and a likelihood approach implemented in TM3 [27]. Estimates were generated using the software NeEstimator [28]. For TM3, we employed 100,000 updates and a maximum N e of 20,000.
For all analyses, we assumed that G. fuscipes undergoes approximately 8 generations per year using observations from colony flies (~7.3 generations per year, [29] 8.5 generations per year at 25°C, [30] and those reported in other studies of the palpalis group (G. palpalis gambiensis and G. palpalis palpalis) in Guinea and Equatorial Guinea [5,31]. All populations were evaluated at an interval of at least one year (~8 generations apart). For four populations (BU, OK, BN, MK), we generated estimates at two different sampling intervals (0 to 8 generations, and 0 to 12 generations).
For each temporal sample in all seven populations, we also tested for an excess of heterozygosity relative to observed allelic diversity, which may be indicative of a recent bottleneck [32]. For each temporal sample, tests of heterozygosity excess were performed separately for each microsatellite locus. Significance was assessed across loci using Wilcoxon's test, which is the most appropriate test given the number of microsatellite loci evaluated. All tests were performed using the program BOTTLENECK [33].

Marker validation and diversity
We genotyped a total of 667 tsetse flies at 16 microsatellite loci. We detected 17 values of F IS (out of 288) that exhibited significant departures from HWE at p < 0.05 (Additional file 1: Table S1). Assessed by locus, the number of significant F IS values observed was consistent with chance at an overall value of p < 0.05. Assessed by population, the number of significant F IS values observed was consistent with chance for all populations except the sample representing generation 12 from BN. Following sequential Bonferroni correction, only one locus pair exhibited significant linkage, and only in one population, confirming previous work showing that these loci were unlinked [10,14].
Microsatellite diversity was lowest in the samples from Mukongoro (MK) and highest in the samples from Bunghazi (BN). Allelic richness ranged from 3.0 to 4.4 and expected heterozygosity (H E ) ranged from 0.418 to 0.609, ( Table 1). MtDNA haplotype diversity was relatively low across samples with the number of haplotypes ranging from 1 to 4. As expected, nucleotide diversity was generally higher in populations from the zone of contact which were composed of flies with both northern and southern ancestry. We detected only two haplotypes that had not been previously reported [10]. Both of these haplotypes were recovered in population OT and differed by just one substitution from Hap26 or Hap27 [10].

Temporal variation in genetic diversity
Variation in allele frequencies by population and locus are depicted in Figure 2. Genetic differentiation between samples taken from the same population but at different times was extremely low, and uniformly lower than the differentiation observed between populations. D EST averaged 0.001 for temporally-spaced samples within populations, compared to 0.308 between populations (Additional file 2: Table S2). The largest values of D EST among temporallyspaced samples were observed in Masindi (MS generation 0 vs.13, D EST = 0.019 ± 0.022) and Otuboi (OT generation 0 vs. 11, D EST = 0.013 ± 0.007).
Mitochondrial haplotype frequencies also exhibited little change over time ( Figure 3). We observed a significant change in haplotype frequencies only between the two temporally spaced samples from Junda (JN, p = 0.046). This was attributable to the loss of the two least common haplotypes in the sample representing generation 13.
An analysis of molecular variance using both microsatellite allele frequencies and haplotype frequencies suggested that differences between temporally-spaced samples explained an insignificant amount of the overall genetic variation ( Table 2). Differences among sites, on the other hand, contributed significantly to the overall variation in genetic diversity. The percent variation explained was greater in the case of mtDNA data, compared to microsatellite data.

Effective size
Estimates of N e were generated for the seven populations based on microsatellite allele frequency changes observed among samples collected at different times from the same population. Momentum-based estimates ranged from 216 to infinity, but only the estimate from OT was bounded by a 95% confidence interval that did not include infinity (Table 3). Likelihood estimates ranged from 152 to 19,550 and all estimates were bounded by 95% confidence intervals that included 20,000, the maximum value of N e considered (Table 3).
For populations MK and OK, estimates of N e were similar regardless of whether the calculations were performed using data for generations 0 and 8 or generations 0 and 12. In populations BN and BU, however, estimates of N e derived from the moment method differed by an order of magnitude depending on whether the sample representing 8 generations or 12 generations was included. In population BN, the estimate of N e generated by the Likelihood method was similarly unstable.  Table 4). Samples from Busime (BU) and Mukongoro (MK) exhibited the least evidence for past bottlenecks (all but one p-value >> 0.05), however power to detect a bottleneck may have been low in MK on account of relatively low genetic diversity (Table 1).

Discussion
We assessed changes in genetic composition of seven tsetse populations in southeast Uganda in order to gain insight into the population dynamics of G. f. fuscipes. In general, our results provide evidence for temporal stability of G. f. fuscipes populations over the one to two year period that we examined. With the exception of just one or two populations discussed below, mitochondrial haplotype frequencies and microsatellite allele frequencies exhibited little change over time and effective population sizes were generally large.
Compared to other riverine species of tsetse, estimates of N e for G. f. fuscipes were similar to or larger than estimates for G. palpalis palpalis in Equatorial Guinea [31] and 2 to 3 orders of magnitudes larger than estimates for G. p. gambiensis on islands off the coast of Guinea [5]. Values of N e for G. f. fuscipes populations were also generally larger than estimates for a savannah species, G. pallidipes, in Kenya [34]. The large effective population sizes and overall stability of G. f. fuscipes populations support the hypothesis [35] that seasonal variation in tsetse numbers, in which larva develop in utero, should be relatively small, since they do not depend on surface water or moist media for breeding. Nonetheless, the lack of variation in genetic structure over time is surprising given the reduced abundance of G. f. fuscipes observed during the dry season in Uganda [12]. To reconcile our results with this observation, which may reflect the low efficiency of traps used for monitoring [36,37], we suggest that populations of G. f. fuscipes in dry season refugia remain large, and that seasonal invasion of marginal wet-season habitat (e.g., at Mukongoro, Bunghazi) must occur in waves of tsetse that are large enough to be representative of the refugia population. Large populations of pupa, which develop in  the ground over a period of weeks, may also help to ensure the continuity of tsetse populations and would contribute to reducing the variance in genetic changes over time.
In contrast to the other populations, estimates of N e were low for populations MS and especially OT, where both moment and likelihood methods produced values of only about 200. These values could be indicative of small populations. N e may also be influenced by overlapping generations and temporal variance in reproductive success as well as the forces of selection, mutation and migration. In this study, however, the low values of N e observed in these populations probably reflected small differences in the location of trapping sites used for the two temporal samples. Generation 13 from MS was sampled at a distance of about 4 km from the original site at which generation 0 was sampled. Likewise, generation 11 from OT was sampled at a single site that was 11-20 km from the relatively widely dispersed sites from which generation 0 was sampled. Thus, for these sites, which were the only two sites sampled at different locations across years, fine-scale spatial genetic variation could be responsible for the apparent temporal variation in gene frequencies, thus depressing estimates of N e .  Given that genetic variation in MS and OT samples can probably be attributed to microgeographic variation, the change in genetic composition of the population at JN likely reflects the only significant temporal change observed in this study. Although microsatellite allele frequencies were largely invariate, mtDNA haplotype frequencies here differed significantly between generation 0 and generation 13. Junda (JN), along with sites BN and MS, lies along a narrow zone of contact between two long-diverged and historically-isolated groups of G. f. fuscipes [9,10]. In 2008, populations at all three of these sites harbored both "southern" and "northern" mtDNA haplotypes. Interestingly, in Junda, individuals with the "southern" haplotypes disappeared from the sample after 13 generations. This could be due to a particularly small population of females and stochastic variation in female reproductive success, although in tsetse, the latter is more likely to be true among males than females [5]. Mating success can also be influenced by Wolbachia, a symbiont that may impose mating barriers due to cytoplasmic incompatibility between infected and uninfected tsetse individuals [38], thereby biasing mating in favor of infected females and potentially producing mitochondrial sweeps [39]. Given the change in mtDNA observed at Junda, flies here should be examined for Wolbachia. If present, the zone of contact in Uganda may provide a unique opportunity to monitor symbiont-induced population changes over time.

Additional material
Additional file 1: Table S1. F IS values for the 16 microsatellite loci. Significance was assessed at p < 0.05 (*) and a Bonferroni-corrected value of p < 0.0028 (bold). Low variability precluded calculation of F IS in some populations (n/a).
Additional file 2: Table S2. Pairwise estimates of genetic differentiation (Jost's D EST ) between samples taken from seven populations of G. f. fuscipes. Estimates of differentiation (below diagonal) and associated standard error (above diagonal) between populations of flies sampled at the same site but different times are shaded in grey. 1 Faculty of Science, Gulu University, Uganda. 2 Department of Ecology and Evolutionary Biology, Yale University, New Haven, Connecticut, USA. Ne was calculated using the moments based temporal method of Waples (1989) and the likelihood approach of Berthier et al. (2002). For four populations, Ne was calculated using samples collected at intervals of both 8 and 12 generations from the initial sampling.