Glossina palpalis palpalis populations from Equatorial Guinea belong to distinct allopatric clades

Background Luba is one of the four historical foci of Human African Trypanosomiasis (HAT) on Bioko Island, in Equatorial Guinea. Although no human cases have been detected since 1995, T. b. gambiense was recently observed in the vector Glossina palpalis palpalis. The existence of cryptic species within this vector taxon has been previously suggested, although no data are available regarding the evolutionary history of tsetse flies populations in Bioko. Methods A phylogenetic analysis of 60 G. p. palpalis from Luba was performed sequencing three mitochondrial (COI, ND2 and 16S) and one nuclear (rDNA-ITS1) DNA markers. Phylogeny reconstruction was performed by Distance Based, Maximum Likelihood and Bayesian Inference methods. Results The COI and ND2 mitochondrial genes were concatenated and revealed 10 closely related haplotypes with a dominant one found in 61.1% of the flies. The sequence homology of the other 9 haplotypes compared to the former ranged from 99.6 to 99.9%. Phylogenetic analysis clearly clustered all island samples with flies coming from the Western African Clade (WAC), and separated from the flies belonging to the Central Africa Clade (CAC), including samples from Mbini and Kogo, two foci of mainland Equatorial Guinea. Consistent with mitochondrial data, analysis of the microsatellite motif present in the ITS1 sequence exhibited two closely related genotypes, clearly divergent from the genotypes previously identified in Mbini and Kogo. Conclusions We report herein that tsetse flies populations circulating in Equatorial Guinea are composed of two allopatric subspecies, one insular and the other continental. The presence of these two G. p. palpalis cryptic taxa in Equatorial Guinea should be taken into account to accurately manage vector control strategy, in a country where trypanosomiasis transmission is controlled but not definitively eliminated yet.


Background
Following the London Declaration on Neglected Tropical Diseases, Human African Trypanosomiasis (HAT) has been targeted for elimination by 2020 [1]. As Trypanosoma brucei gambiense infection actually causes 98% of the total HAT cases (the remaining are due to T. b. rhodesiense), attention must be focused on this subspecies. Among the elimination strategies, vector control can play an important role [2][3][4], especially in isolated populations, which can be targeted for direct intervention avoiding the reinvasion from neighbouring zones. Islands represent an ideal setting for such strategies as demonstrated by the eradication of Glossina sp. in Unguja [5] and in Principe Islands [2] after a few years of sustained control.
Tsetse flies of the palpalis group (Nemorhina subgenus) are major vectors of T. b. gambiense in West Africa [2]. This group comprises two allopatric subspecies: G. p. palpalis and G. p. gambiensis, which probably derived from an ancestral palpalis species which was isolated in several geographic points when its riverine habitat declined during the last glacial maximum [6,7]. Cumulative evidences support the recognition of G. p. gambiensis and G. p. palpalis as valid specific taxa. For example, using data from the mitochondrial gene cytochrome oxidase 1 (COI), the average genetic distance observed between G. p. palpalis and G. p. gambiensis sequences was 6.6%, which is well above the threshold of 2% divergence for inter-species comparisons [8][9][10][11]. Moreover, experimental crosses between these subspecies yielded sterile males in the offspring [12]. The phylogenetic situation is more complex since recent genetic analyses suggested the existence of at least two distinct cryptic species within G. p. palpalis [10,11]. One circulates in the Western part of Africa (named as Western African Clade or WAC), and the other in the continental part of Equatorial Guinea and the Democratic Republic of Congo (DRC) (Central African Clade or CAC). According to the available data, both types are sympatric in the Fontem focus of Cameroon [11].
In Equatorial Guinea, four historical HAT foci are classically defined: one insular (Luba, on Bioko Island) and three on the mainland (Rio Campo in the north; Mbini in the centre and Kogo in the south) [13]. Due to sustained control measures, less than 10 HAT cases are being detected every year in the three continental foci since 2009 and no cases have been recorded since 1995 in Luba. Vector control activities were never implemented in Bioko Island, and parasite elimination in humans relied on active screening of the population at risk and systematic treatment [14]. Therefore, high densities of G. p. palpalis have been observed in some localities at the south of Luba district and moderate densities in others of the epicentre of the focus [15]. Moreover the presence of T. b. gambiense has also been reported in tsetse flies of Luba despite the absence of human infections, which could be attributed to the existence of reservoirs in the wild fauna, cryptic human infections and/ or low sensitivity of available diagnostic tools [15][16][17]. Because vector control is a key parameter to completely eradicate the parasite [2,4,5] a deep knowledge of the biology of the tsetse fly is a crucial prerequisite. In such a context, the genetic characterization of the Glossina palpalis palpalis, never performed so far in Luba, has become indispensable.
In this study we combine information from both mitochondrial DNA (mtDNA-COI, ND2 and 16S-genes) and nuclear ribosomal (rDNA-ITS1) markers to investigate the phylogeographic origin of G. p. palpalis in the focus of Luba, Bioko Island, using tsetse flies samples captured in a previous epidemiological study [15]. MtDNA has been extensively used in population and evolutionary biology of insects [18][19][20] and metazoa in general [21] due to their particular features: relative ease isolation, simple sequence organization, maternal inheritance, absence of recombination and rapid rate of sequence divergence allowing discrimination of recently diverged lineages [22]. On the other hand, the rDNA internal transcribed spacer 1 (ITS1) is a useful marker for both closely related species and also intraspecific populations of insects [23][24][25].

Sample collection
Fly sampling was carried out in September/October 2007 from five areas known to harbour G. p. palpalis (Avendaño, Drumen, Fortuny, Boloco and Las Palmas). We employed monopyramidal traps [26], which have been successfully applied for vector control and entomological surveys in Equatorial Guinea [27][28][29]. Details about trap distribution are provided elsewhere [15]. Tsetse flies collected were individually stored in absolute ethanol in the field until processed in the laboratory. Species identification was undertaken using the key of Brunhes et al. [30]. Tsetse flies were sent to the National Centre of Tropical Medicine, Institute of Health Carlos III (Madrid, Spain) for subsequent molecular analysis.

Molecular analysis
DNA was extracted from whole flies with SpeedTools Tissue DNA Kit (Biotools, B & M Labs, S.A) following the manufacturer instructions. We analysed three mtDNA (ND2, COI and 16S) and one nuclear (ITS1) markers in our study. COI, 16S and ITS1 sequences were amplified with primers described previously [11,31] whereas new specific primers for ND2 gene from G. p. palpalis were designed (Additional file 1: Table S1).
PCR reactions were performed with 2 μl of each template DNA, 1X buffer (10 mM Tris-HCl, 1.5 mM MgCl 2 , 50 mM KCl, pH 8.3), 100 μM of each dNTP, 0.5 μM of each primer, 1 U of Fast Start Taq DNA polymerase (Roche) and double distilled water (DDW) until reaching 50 μl final volume. The thermal cycling programme started with initial denaturation step of 2 minutes at 95°C, followed by 35 cycles (30 seconds at 95°C, 30 seconds at 55-60°C and 1 minute at 72°C) and a final extension step of 5 minutes at 72°C. Results were visualized in 2% agarose gel stained with ethidium bromide under UV irradiation. After this check, we sent the amplified products to Secugen Sequencing and Molecular Diagnostics (Madrid, Spain) where they were sequenced using Sanger method.

Phylogenetic analysis
MEGA version 5 software [33] was used to calculate the pairwise and average distances expressed as number of nucleotide substitutions per site. Phylogeny reconstruction was performed with Maximum Likelihood (ML), Distance Based Neighbor-Joining (NJ) and Bayesian Inference (BI) methods. Model of evolution was inferred from data using the Find DNA/Protein Model tool of MEGA5. This software implements the Bayesian Information Criterion (BIC) [34] and corrected Akaike Information Criterion (AIC) [35] to measure the goodness-of-fit of the 24 models available. Unless otherwise specified, we applied the model selected by BIC to perform phylogenetic analysis (Table 1). For these three inference methods, all positions containing the missing data were eliminated.
MEGA5 software was also used to construct ML and distance based trees of COI, ND2 and concatenated (COI + ND2) sequences. For ML, the bootstrap method (500 replications) was chosen to test the robustness of the trees [39]. The selected substitution model was applied and rates among sites were treated as Gamma distributed with five discrete categories. We assumed no invariant sites. The ML heuristic method employed was the Nearest-Neighbor-Interchange (NNI).
Distance trees were constructed using the NJ method [40]. Gamma shape parameter was estimated from data. The bootstrap consensus tree was inferred from 2000 replicates.
BI was implemented with BEAST software [41]. This software uses Metropolis-Hasting Markov Chain Monte Carlo algorithm [42]. The default settings were generally used (1x10 7 generations and log parameters sampled each 1000 steps). The first 25% of trees generated was discarded as 'burnin'. Yule process was implemented as tree prior, a simple model appropriate when studying speciation [43,44].

Assessment of genetic diversity
Under our sampling conditions, we assume that undetected haplotypes can exist. In order to assess the number of these unseen haplotypes in our studied area (Luba focus), we used two estimators that calculate haplotype richness: i) the first-order Jackknife richness estimator [45] calculated as followed: Jack1 = S 0 + a1(N-1)/ N, where S 0 = observed number of haplotypes, a1 = number of haplotypes detected in only one fly, N = total number of flies and ii) the Bootstrap richness estimator [46] boot = S 0 + Σ(1-p i ) N , where S 0 = observed number of haplotypes, p i = frequency of ith haplotype and N = total number of flies. Analysis was conducted in R software, with the specific 'vegan' package (http://cran.r-project.org, http:// vegan.r-forge.r-project.org).

Phylogeographic analysis
To infer the haplotype relationships within the data sets, the median-joining network algorithm available in NET-WORK v4.5.1.0 was performed [47], which combines the topology of a minimum spanning tree with a parsimonious search for the missing haplotypes.

Mitochondrial markers
We amplified the mtDNA of COI (622 bp), ND2 (501 bp) and 16S (213 bp) genes for a total of sixty tsetse flies, all coming from Equatorial Guinea. Figure 1 (COI + ND2) and Table 2 show the distribution of mtDNA haplotypes regarding the sampling location. In COI gene six different haplotypes were observed with an overall genetic distance of 0.4%, ranging from 0.2% to 0.8%, in terms of number of base substitution per site. Haplotype 1 predominated over the rest (44/60), haplotype 2 was found in 9 individuals, whereas the other haplotypes were found once or twice (Table 3). For ND2, fifty-four sequences were obtained and four haplotypes were detected with an average distance of 0.4%, ranging from 0.2% to 0.5%. Haplotype 1 was the most Analysis of genetic diversity seemed to be assessed with a reasonable accuracy. The number of estimated haplotypes, n = 13.93 and 11.86 calculated with Jackknife and Bootstrap estimators respectively, was close to the number of observed haplotypes (n = 10) when all the samples are considered as a unique population (Additional file 2: Table S2). This means that the haplotype structure of the population is correctly estimated even with low sample size and underestimation of the haplotype richness.
Since very limited data are available regarding the ND2 and 16S sequences in G. p. palpalis from different geographic origins, we decided to focus the phylogenetic reconstruction in COI gene. To construct more comprehensive ML, distance based NJ and BI trees we included previously published sequences of G. p. palpalis from DRC (EU591840-2), Cameroon (EU591829-31, EU591860 and EU591865), Liberia (EU591857-9), Togo (EU591838-9), Ivory Coast (EU591846-8 and EU591832-3), Burkina Faso (EU591856) and two continental foci of Equatorial Guinea, Kogo and Mbini (EU591825 and FJ767873-6). The source of these sequences is detailed elsewhere [10,11]. G. morsitans (GQ255905) was used to root the tree and G. p. gambiensis (EU591851) sequence was also placed as outgroup. Analysis was based on a total fragment of 622 bp of the COI gene.
ML, distance based NJ and BI trees showed common topologies. All mirrored two major geographic separated clusters from Central (CAC) and Western Africa (WAC) (Figure 2). NJ and BI trees showed more support in this split: 97% bootstrap and 0.9999 posterior probability in BI trees, respectively, than ML inference (88%). Within  For each location is detailed the number of times a given haplotype has been obtained.
WAC, organisation was unclear and again NJ and BI methods yielded a more robust clustering. ML inference method failed to yield a clear discrimination between clades from Western Africa (50% bootstrap) whereas NJ and BI trees split samples from this area in different sister sympatric groups (with 99% bootstrap and 1.0000 posterior probability). Clusters obtained by the NJ and BI methods within the WAC are probably not artifacts since both trees separate exactly the same specimens and in spite of the lower power of the ML tree, the same trend is visible. Regardless of the inference method used the COI haplotypes from the Luba focus grouped within the WAC I cluster in all trees, in clear contrast with the haplotypes from Kogo and Mbini foci (mainland Equatorial Guinea) that belong to the CAC cluster. In accordance to phylogeny trees, network analysis shows the clear separation of the two major geographic clusters (WAC and CAC) ( Figure 3) together with the two subgroups corresponding to DRC and continental Equatorial Guinean haplotypes within the CAC (CAC I and II in Figure 2).
Regarding ND2 gene, five haplotypes were available in Genbank, one from Ivory Coast (EU591895), one from Liberia (EU591884), two from Cameroon (EU591897-8) and one from Equatorial Guinea (EU591905). In accordance with COI data, ML, distance based NJ and BI trees constructed with ND2 and COI + ND2 concatenated genes (Figure 4), exhibited a clear separation between the sequences of mainland Equatorial Guinean haplotypes and those from the Luba focus.
Nuclear ITS1 marker ITS1 size polymorphism was assessed for 59 flies. All PCR products analysed exhibited a unique size around 240 bp, similar to the size reported by Dyer et al. [11] when this fragment was amplified from flies belonging Figure 3 Geographic representation of the COI haplotypes. A) Geographic distribution of the COI lineages; B) Haplotypes network derived from 30 haplotypes of the G. p. palpalis complex. Haplotypes are represented by circles and their frequency is proportional to the area. Network diagrams created using Phylogenetic Network software from Fluxus and using the Median Joining method as described previously [47]. to the WAC. After sequencing, two similar genotypes defined by the (AT) n microsatellite were found. Genotype 1 was characterized by (AT) 9 G(TA) 2 region, whereas genotype 2 exhibited (AT) 10 pure repeats (Additional file 4: Figure S2). As shown in Table 3, genotype 1 was dominant in all the 5 sampled areas of the Luba focus. The two genotypes found were aligned with ITS1 sequences of G. p. palpalis from mainland Equatorial Guinea (Kogo: J767886 and Mbini: J767887-J767888), DRC (J767892 and J767893), Guinea Conakry (EU59 1930), Gambia (EU591931), Burkina Faso (EU591932), Togo (EU591933), Liberia (EU591934) and Ivory Coast (EU5991935) (Additional file 4: Figure S2). The main genotype, found in 44/59 flies and defined by a (AT) 9 G (TA) 2 domain, was not detected in the other published sequences. The second genotype, defined by a (AT) 10 repeat, shared 100% homology with the sequence of a Togo sample.

Discussion
Luba focus is located at the edge between the two main clades of G. p. palpalis, one named WAC (for Western Africa Clade), including flies from Cameroon, Burkina Faso, Ivory Coast, Liberia and Togo; and the other, named CAC (for Central African Clade) represented by flies from the Democratic Republic of Congo and the mainland region of Equatorial Guinea. Phylogenetic analysis, using mtDNA markers, allowed us to cluster the G. p. palpalis population from Luba within the WAC. All the samples of this work were unambiguously included in this group and separated from the CAC regardless the phylogenetic inference method used.
Our results are consistent with the geological history of Bioko, originated by volcanic eruptions in the lower Tertiary period, around 60 million years BP. Bioko is a part of an archipelago which pertains to a large volcanic fracture originating in the South of Lake Chad and extending to Mount Cameroon in the continent [49]. Although politically belonging to Equatorial Guinea, Bioko lies closer to the Cameroon coast (ca. 30 km) rather than to the rest of Equatorial Guinea territories (more than 200 km from mainland and 700 km from Annobon island). During geoclimatic events of Quaternary period, Bioko was linked to mainland given the lower sea level, presumably to Cameroonian coast because of its geological origin. At the end of the last glaciation (around 12,000 years BP), Bioko was isolated by flooding and separated from the continent [50]. It is probable that G. p. palpalis population of Bioko was isolated from Cameroon coast after that event, as suggested for other insect vectors such as Anopheles melas and Simulium yahense [51,52].
Both mtDNA and nuclear markers show a very low genetic intra-variability. MtDNA genes data polymorphism did not exceed 0.8% and ITS1 sequences only yielded two closely related genotypes. For ITS1, the more abundant one is apparently exclusive of Luba focus, suggesting a certain degree of isolation. This is in accordance with previous studies, which found that G. p. gambiensis populations separated from the continent only by 4-5 km of sea show clear evidences of complete isolation. These results are based on wing landmarks, DNA mitochondrial markers and microsatellite dataset [53,54]. Our data should be also completed with an exhaustive coalescent population genetics analysis to support or not the hypothesis of an on-going allopatric speciation in Luba.
Because of the presence of these two sympatric ITS1 genotypes in Bioko Island one could speculate about the past occurrence of at least two separated migration events, as suggested for Anopheles gambiae in Bioko [55]. The ITS1 tandem array sequences are expected to be quickly homogenized by concerted evolution in interbreeding populations, whereas differences are usually observed between non-interbreeding ones [56,57]. However, the low difference observed between our ITS1 genotypes is based in the AT repeats, where a slight heterogeneity can be expected since microsatellite sequence could be evolving faster than the homogenization process. This phenomenon has been previously described in the tandem repeated U2 snRNA gene of primates [58,59].
Other possible explanation for the existence of two distinct ITS1 genotypes is a more recent reintroduction of tsetse flies from continent. Bioko Island is located at 30 km West of the Cameroonian coast, the closest area of the continent. Estimated active dispersal of Glossina sp. in one day is no longer than 1.3 km (in 15-30 min/ day) [60,61]. Additionally, these flies are usually unable to flight for long periods but rather in short bursts lasting between 1 and 2 minutes [62,63]. These observations make extremely unlikely a recent (posterior to the glaciations period) reproductive contact between Luba and other continental foci due to active dispersal. However, human-mediated transportation of flies may not be disregarded. This situation allowed the reinvasion of G. p. palpalis in Principe Island in 1956 despite its eradication in 1914 using mobile traps carried by workers [2]. The distance between this island (ca. 240 km) and mainland is much larger than that of Bioko. It is therefore expected that human movement between Bioko and the coast of Cameroon has been much more frequent since the island was first colonized by the Bubi ethnical group, at the end of the last glaciation [64]. Finally, although the existence of ITS1 hybrid forms was initially suggested in G. p. palpalis of Equatorial Guinea [11], this hypothesis was later rejected by the same authors using microsatellite markers [10]. A set of highly polymorphic nuclear markers should be applied in the future to definitively test this hypothesis in Luba.
The genetic markers used to assess geographical structuring of G. p. palpalis demonstrate the existence of two allopatric taxa in Equatorial Guinea, one in the insular focus, Luba, and the other in two of the mainland foci (Kogo and Mbini). G. p. palpalis populations from Rio Campo, the third mainland focus of the country, have not been studied yet. In Equatorial Guinea, vector control was implemented in Kogo and Mbini from 2002 to 2009, whereas in Luba the active detection and treatment of HAT cases was the only control method employed [14]. Vector control has proved to be efficient at controlling HAT transmission [2,5,13] since tsetse flies show particular features that make them highly susceptible to direct interventions. Firstly, tsetse flies have a very low reproductive rate, given that a female individual has probably only one reproductive mating in its life and deposits only one larva per generation (up to 12 generations with 9-10 days of interval for a lifespan of 2-3 months) [2]. As a result of this low reproductive rate, the population of the vector is usually low comparing to other diptera and small increases in mortality can lead to control [4,65] or even to population extinction [66]. Secondly, the active dispersal of tsetse flies is generally low [60,61] resulting in a more difficult re-colonization of cleared areas. Thirdly, the genetic variability observed within populations tends to be reduced as well [67,68], probably as a consequence of both low reproductive rate and limited dispersal capacity, making more difficult for the selection of new attributes such as resistance to insecticides. There are, however, other behavioural features such as feeding preferences or trap-avoidance that can vary at subspecies and even at population level. For example, within G. p. palpalis, diverse feeding preferences across foci of Cameroon have been observed [69]. Although this observation could be attributed to the opportunistic feeding behaviour of G. p. palpalis [70], it could be also associated with the existence of genetically different G. p. palpalis populations, given the probable isolation of these foci [71]. Indeed, different feeding patterns in two sympatric populations of G. p. palpalis and G. p. gambiensis have also been described in preforest areas of Cote d'Ivoire [72], demonstrating that this feature can differ between closely related subspecies. The feeding behaviour of the vector can be crucial to design an effective vector control campaign since it provides valuable information to understand the epidemiological cycle of the parasite at local level.