Multilocus Analysis of Divergence and Introgression in Sympatric and Allopatric Sibling Species of the Lutzomyia longipalpis Complex in Brazil

Background Lutzomyia longipalpis, the main vector of visceral leishmaniasis in Latin America, is a complex of sibling species. In Brazil, a number of very closely related sibling species have been revealed by the analyses of copulation songs, sex pheromones and molecular markers. However, the level of divergence and gene flow between the sibling species remains unclear. Brazilian populations of this vector can be divided in two main groups: one producing Burst-type songs and the Cembrene-1 pheromone and a second more diverse group producing various Pulse song subtypes and different pheromones. Methodology/Principal Findings We analyzed 21 nuclear loci in two pairs of Brazilian populations: two sympatric populations from the Sobral locality (1S and 2S) in northeastern Brazil and two allopatric populations from the Lapinha and Pancas localities in southeastern Brazil. Pancas and Sobral 2S are populations of the Burst/Cembrene-1 species while Lapinha and Sobral 1S are two putative incipient species producing the same pheromone and similar Pulse song subtypes. The multilocus analysis strongly suggests the occurrence of gene flow during the divergence between the sibling species, with different levels of introgression between loci. Moreover, this differential introgression is asymmetrical, with estimated gene flow being higher in the direction of the Burst/Cembrene-1 species. Conclusions/Significance The results indicate that introgressive hybridization has been a crucial phenomenon in shaping the genome of the L. longipalpis complex. This has possible epidemiological implications and is particularly interesting considering the potential for increased introgression caused by man-made environmental changes and the current trend of leishmaniasis urbanization in Brazil.


Introduction
Speciation events are the result of a complex array of interesting and dynamic biological processes that remain only partially understood [1]. The intricate balance among mutation, recombination, gene flow, genetic drift and natural selection can either unify or differentiate genetic variation, and this consequently may or may not promote the appearance of new species [2][3][4][5]. The evolutionary history described by individual genes may be quite variable, and the pattern of relationships among closely related species can be discordant [6][7][8]. In principle, the use of multiple loci should give a more complete picture of the history of divergence of species complexes, and comparisons across genes can reveal whether all loci fit a simple model of dichotomic phylogeny. Use of multiple loci can also reveal whether retention of ancestral polymorphisms, introgression or selective pressures can explain incongruities in the evolutionary histories of species complexes [5,8,9].
Gene flow during a speciation process can be evidenced when some loci show little or no differentiation, while others show a large level. In the last decade, many studies have been carried out in species that have recently diverged, and it appears that divergence and speciation may often occur in the presence of gene flow [10][11][12][13][14][15][16][17][18][19][20][21]. This is also true for insect disease vectors, where the process of divergence with gene flow is particularly interesting beyond the standard evolutionary viewpoint. Genes involved in vectorial capacity, insecticide resistance and adaptation to different environmental conditions could be introgressing between sibling species with important epidemiological consequences [20,[22][23][24][25][26].
The most important neotropical vector of Leishmania infantum, the causative agent of American visceral leishmaniasis, is the sand fly Lutzomyia longipalpis (Diptera: Psychodidae), a complex of sibling species with a large distribution area ranging from northern Argentina and Uruguay to Mexico [27][28][29][30][31][32][33][34][35][36]. In Brazil, despite some incongruence among genetic markers [32,34], an integrative analysis using a combination of biochemical, behavioral and molecular traits strongly supports a number of sibling species having different levels of differentiation [33].
The Brazilian populations of L. longipalpis s.l. can be divided in two main groups. The first group is genetically homogeneous and widely spread, and probably represents a single species. Males of this species produce Burst-type copulation songs and the Cembrene-1 pheromone. The second group is very heterogeneous and likely represents a number of putative incipient species with more restricted distributions. These sibling species produce different subtypes of Pulse-type copulation song in combination with different sex pheromones (Germacrene, Himachalene, Cembrene-1 and Cembrene-2) [33].
The coexistence of sibling species in an overlapping geographic area is one of the best pieces of evidence for the existence of a species complex [37,38]; in at least three localities in Brazil, siblings of L. longipalpis s.l. are present in sympatry [33]. This is true for the Brazilian municipality of Sobral (Ceará State, Northeast Brazil), where two L. longipalpis sibling species were observed. In this locality, males of these two species can be distinguished by the number of pale spots on the abdomen (one or two pairs of spots: ''Sobral 1S'' and ''Sobral 2S,'' respectively). Crossing experiments show that these two siblings have strong reproductive isolation, which is consistent with the fact that their males produce different pheromones and copulation songs [29,35,39,40]. In addition, molecular markers such as microsatellites and nuclear genes clearly indicate that Sobral 1S and Sobral 2S represent two sympatric species in Sobral [41][42][43][44]. However, some of these molecular markers also show signs of introgression [41,43], and this could explain the differences in interpretation among early studies regarding the status of the Brazilian populations [32,34].
In the current study, we conducted a multilocus approach using 21 nuclear loci to estimate and compare levels of divergence and gene flow between the sympatric siblings from Sobral and two allopatric species from the localities of Lapinha (Minas Gerais State) and Pancas (Espírito Santo State) in Southeast Brazil. Pancas and Sobral 2S are probably populations of the same species, whose males produce Burst-type songs and the Cembrene-1 pheromone, while Lapinha and Sobral 1S are two putative incipient species that share the same pheromone (Germacrene) and produce, respectively, Pulse song subtypes P2 and P3 [33]. In the absence of interbreeding, sympatric populations of two species should not be more similar to each other genetically than to allopatric populations of these two species. This kind of comparison using a larger number of unlinked loci allows to distinguish among common ancestry and the effects of introgression. Moreover, it brings new insights of how this is shaping the genome of the L. longipalpis species complex.

Methods
Sand flies were collected from Sobral (Ceará state) (3u419S, 40u209W) in northeastern Brazil, Pancas (Espírito Santo state) (19u139S, 40u519W) and Lapinha (Lagoa Santa, Minas Gerais state) (19u389S, 43u539W) in southeastern Brazil ( Figure 1). A permit for sand fly collection in Brazil was obtained from the Brazilian Ministry of Environment (SISBIO #26066-1). Sand flies were captured using CDC light-traps near human habitation with permission from local homeowners. In addition, the collections were usually supported by the local vector surveillance authorities from local State Health Departments. Samples were identified according to [45], and only males were used to avoid misidentification, as females are difficult to distinguish from other closely related species.
The remaining 18 loci (Table S2) are new markers randomly selected from a screen of cDNA sequences performed in our lab. The selected loci are responsible for different functions (Table S1), show a high similarity at the protein level when compared to Anopheles gambiae and/or Drosophila melanogaster and contain one intron in the latter species under the assumption that they might also be present in the L. longipalpis studied fragment, although that was not always the case.
DNA samples used for PCR were prepared as previously described [24]. PCR was carried out using proofreading Pfu DNA Polymerase (Biotools) for 35 cycles: 95uC for 30 s, 55uC for 30 s, and 72uC for 30 s. PCR products were purified using the Wizard PCR Preps kit (Promega) or Micro Spin S-400 HR (GE Healthcare) and were cloned into the pMOS Blunt-ended Vector (GE Healthcare). The plasmid DNA was prepared using the alkaline lysis method in 96 well micro-plates [47] and filtered in multiscreen filter plates. Sequencing was carried out using the ABI Prism Big Dye Terminator Cycle Sequencing Ready Reaction V3.1 kit (Applied Biosystems) in an ABI 3730 DNA Sequencer by

Author Summary
The sand fly Lutzomyia longipalpis, the most important vector of visceral leishmaniasis in the Americas, is a complex of cryptic species distributed from Argentina to Mexico. There is evidence for the existence of a number of closely related sibling species of this complex in Brazil that differ in their male mating songs, sex pheromones and molecular markers. We compared the levels of divergence and gene-flow (introgression) in sympatric and allopatric pairs of sibling species of this complex using several molecular markers. Our results suggest that the L. longipalpis species complex in Brazil has an intricate evolutionary history with a crucial role for introgression, which varies between loci and is asymmetrical between the sibling species. This introgression potentially has important epidemiological implications because it could affect genes that influence the relative role that the different sibling species play as vectors. In addition, our results are particularly relevant considering the potential for increased introgression caused by man-made changes to the environment and the current trend of urbanization of visceral leishmaniasis in Brazil.
We analyzed our data set under the isolation-with-migration model available with the IM software [12] to discriminate between the relative effects of divergence and gene flow. This multilocus analysis was performed using non-recombining blocks (NRB) of the 21 different nuclear loci. Sites within indels or ambiguous alignment were removed (Table S3). To select the NRB, we used IMgc software [56] for the Sobral 1S and 2S sequences and manually selected the blocks and sequences of Pancas and Lapinha according to the NRB from Sobral to ensure consistency. Some putative recombinant sequences were excluded from the data set. The IM model considers parameters for splitting time (t), bidirectional gene flow after splitting (m 1 and m 2 ) and current effective population sizes (h 1 , h 2 and h A ). To fit the model to our data, the IM software used a Bayesian framework that gave estimates for the posterior probability densities of the model parameters using a Markov chain Monte Carlo simulation [12]. We carried out two sets of analyses, one comparison involving the sympatric siblings from Sobral (Sobral 1S vs. Sobral 2S) and one for the two allopatric siblings (Pancas vs. Lapinha). For each analysis, we ran simulations assuming the Infinite Site model (IS) [57] of sequence evolution recommended for nuclear loci. In addition, for each comparison (sympatric or allopatric) we ran simulations that estimated the population migration rates and simulations that estimated the per locus migration rates. We conducted preliminary runs to determine appropriate priors for subsequent runs. After that we performed three to five independent runs with different starting points to ensure that the values converged on similar estimates, with ten or twenty independent chains under Metropolis coupling. Each chain was initiated with a burn-in period of 200,000 updates, and the total run length of each analysis was between 30 million and 40 million updates, using the following input parameters for simulations: m = 25 and 100, h = 5 and 10, t = 5. Convergence was assessed by estimating the effective sample size, which was always over 50. The six demographic parameters were calculated from the values of the bin with the highest count; these values corresponded to the average calculated from three to five independent runs. For credibility intervals for each parameter, we recorded the 90% posterior density interval, the shortest span that includes 90% of the probability density of a parameter. To convert the parameter estimates into demographic values, we used Drosophila melanogaster synonymous and non-synonymous substitution rates for nuclear genes, 1.56610 28 and 1.91610 29 per site per year, respectively [58]. From the 21 loci used in this study, we calculated a geometric mean of mutation rate per loci per year, m = 1.77610 26 , that was used to rescale the IM parameter estimates. This value was used to convert the parameter estimate t to the number of years since population splitting (t), assuming that L. longipalpis s.l. has about ten generations per year, and to convert the estimates of the population mutation rates (h 1 , h 2 and h A ) into estimates of the effective population sizes (N 1 , N 2 and N A ).

Intra-population nucleotide variation
Our data set included 21 loci. The markers are responsible for different functions and potentially spread throughout the genome of L. longipalpis, based on their positions in A. gambiae and D. melanogaster (Table S1, see also Methods for more details). Table 1 shows a summary of the multilocus DNA polymorphisms for the four L. longipalpis s.l. populations sampled in eastern Brazil ( Figure 1). In general, for each molecular marker, the levels of variation within populations were similar. However, we observed considerable variation in the levels among loci. The CG9297, rpL36, tfIIAL and obp19a loci were the most polymorphic, and sesB, eno, tropC and CG9769 were the least polymorphic of the markers analyzed. Indels were observed in thirteen markers, and all were located in introns. Table 2 shows the Tajima's D statistics [59]. Although not significantly different from zero after Bonferroni's correction, values were negative in most cases, possibly indicating population expansion. Additional tests, such as the Fu's F s [60] and Ramos-Onsins and Rozas R 2 [61] (considered more powerful to detect population growth [62]) also support the expansion hypothesis, mainly in Sobral (Table S4)

Differentiation and introgression between sympatric and allopatric pairs of sibling species
We carried out a comparison of the differentiation between the two pairs of sibling species: sympatric (Sobral 1S vs. Sobral 2S) and allopatric (Lapinha vs. Pancas). As mentioned before, Pancas and Sobral 2S are populations of the Burst/Cembrene species while Lapinha and Sobral 1S are two putative incipient species producing the same pheromone (Germacrene) and Pulse song subtypes P2 and P3 [33]. for all pairwise comparisons. A tree of these four populations based on the mean F ST values clearly shows that the closer genetic similarity of the two Burst-type populations (Pancas and Sobral 2S) at one side and the two Pulse-type populations (Lapinha and Sobral 1S) in the other (Figure 2). In addition, it also shows that although Lapinha and Sobral 1S produce different subtypes of Pulse songs and probably represent incipient species, the overall level of genetic divergence between them is only slightly higher than between the two Burst-type populations.
The sympatric pair shared polymorphisms in almost all loci and only exhibited fixed differences in para (Table 3). On the other hand, the allopatric pair of species only shared polymorphisms in about half of the 21 loci and a total of 17 fixed differences were observed in six loci. The lower level of divergence between sympatric species than allopatric species was also observed after inspection of the observed F ST values (Table 3 and Figure 3). The mean pairwise F ST was more than twice as high for the allopatric pair (0.45360.245) than for the sympatric pair (0.21360.209), and this difference was significant (t = 23.413; p,0.01). While all 21 F ST values were significant in the case of the allopatric pair, five of the values for the sympatric siblings failed to reach significance. Out of 21 loci, only two (sec22 and up) showed higher F ST values between the sympatric pair than between the allopatric populations.
The data in Table 3 and Figure 3 strongly suggest the occurrence of introgression between the sympatric populations. In addition, the estimated F ST statistics were highly variable between loci in both comparisons. Part of the variation was likely caused by differences in the rate of evolution of the different loci. However, as shown in Figure 3, differences in F ST values between the allopatric and sympatric comparisons tend to be lower for loci that have high F ST values in the sympatric pair. In fact, Figure 4A shows that there is a highly significant negative correlation (r = 20.896; p,0.001) between the ranked (from the smallest to the highest) normalized differences in F ST values between the sympatric and allopatric pairs at each loci and the respective rank of F ST values in the sympatric pair at the same loci (see figure legend for more details). No correlation is observed in the case of the allopatric populations (r = 0.231; p = 0.3137) ( Figure 4B). These results suggest that selection is acting as a filter on gene flow and, as a result, introgression is differential, affecting some loci to a much greater extent than others. This shows that, while one can predict that a molecular marker with a high F ST value in a sympatric comparison will most likely show a similar value between two allopatric sibling species, the same is not true for many of the high F ST values observed between allopatric populations. A gene that differentiates sympatric species is likely a good general marker for the complex. The same is not necessarily true for markers that differentiate allopatric species. Figure 5 illustrates this further by comparing the haplotype networks of a few markers (fcop, rpL17A, per, up, sec22 and para) in sympatric and allopatric populations. The first two, rpL17A and fcop ( Figure 5A), represent examples of markers which showed high F ST values in allopatry but low values in sympatry ( Figure 3). The allopatric networks of these two loci show much better separation between the haplotypes, especially in the case of rpL17A, which shows no separation at all in sympatry. On the other hand, in the case of the four loci with high F ST values in sympatry (per, up, sec22 and para), the networks are somewhat similar in the allopatric and sympatric comparisons ( Figure 5B). Therefore, these markers seem to be more reliable for the study of the divergence between populations of the L. longipalpis complex as they avoid the problems caused by strong introgression.
The results presented above show a complex pattern of divergence and gene flow between the Brazilian species of the L. longipalpis complex. In particular, we found evidence of differential introgression among loci in the sympatric pair of species. To explore in more detail the gene flow between sympatric and allopatric relatives, we performed a multilocus analysis using the IM software.

Isolation with migration analysis
Because the IM software requires the absence of recombination within loci, a non-recombining block (NRB) was extracted from each gene alignment, and some putative recombinant sequences were excluded from the data set (Table S3, see Methods). To evaluate the effects of the decrease in fragment length and/or the number of sequences per loci, we conducted polymorphism and    (Table 4) were estimated to infer and compare the population history of the two sympatric and the two allopatric sibling species of the L. longipalpis complex. Figure 6 shows the posterior probability distributions for six parameters with single narrow peaks and bounds that fell within the prior distributions. In all cases, each of the replicates yielded a posterior distribution with identical and well-defined modes. The marginal posterior distribution for h of the sympatric species showed slightly different distributions ( Figure 6A). The current effective population sizes estimated from these values indicate a three-fold increase for Sobral 1S and a two-fold increase for Sobral 2S since their splitting from the ancestral population. On the other hand, the h parameters for Pancas and Lapinha indicate that these two allopatric relatives and the ancestral population have similar effective population sizes ( Figure 6A). The marginal posterior distributions for t suggest that, as expected, the sympatric pair of species split from the ancestral population at approximately the same time as the allopatric pair ( Figure 6B). Therefore, assuming mutation rates similar to Drosophila (see Material and Methods), the splitting event that separated the Burst-type (Sobral 2S and Pancas) and Pulse-type song (Sobral 1S and Lapinha) lineages occurred approximately 0.5 MYA (Table 4).
Shared variation may be the result of gene flow or the incomplete sorting of ancestral polymorphism. The distinction between both possibilities is the main goal of several statistical tests [63] like the one implemented by IM through coalescent simulations [12]. Based on multilocus estimates of population migration rates showing nonzero peaks, our results strongly suggested that after separation, the two lineages have undergone gene flow in both directions for the sympatric and for the allopatric    Figure 6C). Although the effective calculated number of gene migrants per generation indicates a bidirectional migration in both cases, gene flow is five to ten times higher in the sympatric species. In addition, the observed level of introgression is highly asymmetrical, with Sobral 2S receiving about seven times more migrants than Sobral 1S. The estimates of the mean time of migration events [13] indicate that most of them occurred between 0.13 MYA and 0.20 MYA, suggesting lower current levels of gene flow.
We also estimated the per locus migration rates. Figure 7 shows the marginal posterior distributions for each gene. Migration at several loci was observed in both pairs of species. However, although sympatric and allopatric comparisons showed a similar number of loci that exhibited signals of gene flow, in general the allopatric comparison showed much lower levels of gene flow. In addition, most loci showed unidirectional gene flow and introgression from Sobral 1S to Sobral 2S. Weak divergence and extensive gene flow may result in a flat and diffuse posterior distribution [64], and this most likely explains the results observed for a number of loci in the case of the Sobral sympatric siblings.

Discussion
The analysis of the divergence and gene flow at 21 nuclear loci in sympatric and allopatric Brazilian populations of L. longipalpis s.l. shows an intricate evolutionary history. The results suggest that, at least for the Brazilian species, introgressive hybridization has played a crucial role in the speciation process in this complex. In closely related species, this phenomenon seems to be common, and multilocus analyses have been useful in these cases [4,5]. Overall, our results strongly suggest that introgression has played an active role in shaping the genomes of species in the L. longipalpis complex, as previously shown, for example, in Anopheles, Heliconius and Drosophila closely related species (e.g. [8,22,[65][66][67].
The occurrence of sympatric species of the L. longipalpis complex in Sobral was strongly supported by our analysis, and this finding agrees with previous work [32]. Moreover, the high number of shared polymorphisms, some full shared haplotypes in a number of loci and the lack of fixed polymorphisms (except those previously reported for para [44] were the outstanding genetic traits observed in those two sibling species. The extremely variable levels of divergence among loci were the first clue that an evolutionary model of isolation with migration might fit our data [5,14,68].
The comparison between sympatric vs. allopatric populations provide further evidence for the occurrence of introgression between the sympatric species. If the sympatric pair alone had been analyzed, it might have been more difficult to determine whether the low divergence observed in some loci was merely a consequence of retention of ancestral polymorphisms or was the result of gene flow. However, as in some other studies (e.g., [23,69,70], the sympatric/allopatric comparison provided a clearer picture of the divergence process and indicated the occurrence of gene flow, which was confirmed by the IM analysis [12].
Interestingly, the allopatric sibling data also fit the isolation with migration model, which might reflect gene flow coming from populations of the same species subject to introgression in sympatry. The estimates of migration rates, time of divergence and mean F ST between the allopatric sibling species, all satisfy the tentative threshold criteria for species diagnosis proposed by Hey and Pinho [71]. However, in the case of the sympatric species, the situation is a bit less clear with introgression affecting the mean F ST and migration rates.
The observed level of introgression is highly asymmetrical with Sobral 2S receiving about seven times more migrants than Sobral 1S. Asymmetry in gene flow might be caused by differences in population density, and the gene flow direction is often predominantly from the most abundant species to the least [72][73][74]. In fact, our estimates of the h population size parameters suggest that Sobral 1S is larger than Sobral 2S. The evidence for differential introgression among loci fits the mosaic model of speciation (e.g. [8,22,75,76] in sympatric and allopatric species and can be used to identify genomic regions containing genes involved in speciation [77]. Some genes may have freely crossed the species boundaries; these genes include housekeeping and other genes that may not be associated with reproductive isolation or species-specific adaptation (e.g., sod2 and tfIIAL) and were therefore not subjected to selective pressure against introgression [78]. Among the loci used in this study, three (per, cac, and para) are known to control courtship songs in Drosophila [79] and, interestingly, two (per and para) are among those with highest F ST values. Future analysis of other song genes may reveal whether these genes tend to show higher levels of divergence and lower levels of gene flow between the L. longipalpis sibling species.
Assuming that the speciation process between the Burst/ Cembrene and Pulse/Germacrene lineages occurred in allopatry (,0.5 MYA), followed by secondary contact in localities such as Sobral, the evidence for introgression provided by the data suggests that, at first, the period of separation was not long enough to ensure the appearance of full reproductive isolation mechanisms. Following the secondary contact and a period of stronger hybridization and introgression, reinforcement of reproductive isolation might have promoted the evolution of more efficient mate discrimination, and other mechanisms of isolation could have taken place [67,80]. Differences in male sex pheromones [35], copulation songs [39], locomotor activity rhythms [81] and life cycle traits [82] between Sobral 1S and 2S indicate that selection might be playing an active role in the divergence process of the two sibling species. In fact, there is evidence that the reproductive isolation between the sympatric Sobral siblings is stronger than between allopatric siblings of the L. longipalpis complex [29,40], and gene flow might therefore be currently diminished.
In addition to Sobral, two other pairs of sympatric species in Jaíba (Bahia State) and Estrela de Alagoas (Alagoas State) localities were previously described in Brazil. In both cases, a Burst-type song and Cembrene-1 population, coexist in sympatry with a different Pulse-type song sibling [33]. Unlike Sobral, in Jaíba the The network of the loci fcop and rpL17A from the sympatric species of the Lutzomyia longipalpis complex shows mixed a haplotype distribution, unlike the two wellseparated cluster for the allopatric species. This is in agreement with the low degree of divergence between sympatric and high divergence between allopatric species. (B) The loci per, sec22, up and para presented sympatric and allopatric networks with haplotypes separated by species group, which also corroborates the high values of pairwise F ST (see Table 2). Each L. longipalpis population is represented by one color: Sobral 1S (orange), Sobral 2S (green), Lapinha (red) and Pancas (light green). Colored circles represent unique haplotypes, and their sizes are proportional to their frequencies. Black circles are hypothetical haplotypes. Curved lines represent alternative branching between haplotypes. doi:10.1371/journal.pntd.0002495.g005 Table 4. Estimates of parameters from IM analysis of sympatric and allopatric comparisons. two siblings differ in the type of diterpene isomers that they carry and in Estrela de Alagoas the siblings share the same type of pheromone [83]. Future work might reveal semipermeable species boundaries which would not necessarily involve the same genes, intensity and/or direction of gene flow. Nevertheless, our results indicate that past gene flow has affected several areas of the L. longipalpis s.l. genome. In some cases, genomic regions with suppressed or reduced recombination, as a result of chromosomal rearrangements such as inversions, might have less introgression than colinear regions [84,85]. Genes located in such regions might be involved in adaptive differences or reproductive divergence between siblings, and therefore be filtered by selection, while genes that are not linked to such regions might introgress more readily [8,66]. In L. longipalpis s.l., putative chromosomal inversions between some siblings have been identified [86], but future mapping experiments are needed to reveal which genome regions of this species complex have undergone more introgression than others. Furthermore, the knowledge of introgression patterns throughout the genome is important to understand whether loci related to vectorial capacity can influence the transmission dynamics of Leishmania parasites by the different L. longipalpis sibling species. This will be particularly interesting under an epidemiological point of view considering the potential for increased introgression caused by human-made changes to the environment [87] and the current trend of visceral leishmaniasis urbanization in Brazil [88].

Supporting Information
Table S1 Chromosome positions of the 21 loci in D. melanogaster and A. gambiae and their biological functions and processes. (DOC)     Figure 7. Distribution of marginal posterior probabilities of the migration rates at each locus. Probability densities (P) of migration estimates for 21 nuclear loci for two comparisons, sympatric -Sobral 1S (orange) vs. Sobral 2S (green) -and allopatric -Lapinha (red) vs. Pancas (light green). doi:10.1371/journal.pntd.0002495.g007