Data on the genetic diversity obtained for the N. ephratae populations showed high genetic variability for the species as evidenced by the large number of haplotypes in relation to the total number of samples as well as by the high haplotype diversity. As expected, the nucleotide diversity observed for the COII sequences was greater than the nucleotide diversity found for the 16S. Furthermore, the π value obtained for the COII synonyms sites was higher than the π value for the non-synonym sites (Table 3a). This is because synonym substitutions can produce deleterious modifications in the protein structures; hence, they tend to not settle in the populations. Non-synonym substitutions rarely produce important changes in the gene expression, and thus they suffer low selective pressure and, therefore, tend to be retained and increase the nucleotide diversity indexes for these sites in the populations.
The high genetic variability observed for N. ephratae seems not to be widely shared among the populations because the AMOVA results (Table 5) showed strong genetic structure for the species. That is, there are a lot of genetic variants that are exclusive of one or few regions suggesting low or moderated gene flow among distant populations. This limitation in gene flow may have occurred due to isolation by distance as confirmed by Mantel’s test. The clustering analysis results also corroborated the population structure inference because five clusters with a distribution limited to adjacent dominions were recovered. Cluster 1 shows a larger geographic distribution and occurs in all South American dominions. Great genetic variability with haplotypes little shared among distant geographic regions was also observed very similar to N. corniger [15].
Regarding to the demographic history of the N. ephratae populations, the neutrality tests (Table 3b,c), specifically Fu’s Fs, indicated demographic expansion for the species (considering the whole sample set) and particularly for each neotropical dominion sampled. This also can be inferred from the star shape of the haplogroups in the haplotype network, thus suggesting that a lot of descendent haplotypes had recently risen from the ancestor haplotypes located in the center of the haplogroups.
The haplotype network (Figure 2) also showed four haplogroups presenting clear differences about the geographical distribution of the haplotypes. Haplogroup 1 is entirely composed of haplotypes from the South America dominions while haplogroups 2 and 4 are mostly formed by Central American haplotypes. Some haplotypes included in haplogroup 4 can also be found in northern South America (Venezuela and Trinidad and Tobago), but were not often observed at latitudes below this. Haplogroup 3 is composed of haplotypes from northern South America (Trinidad and Tobago and French Guiana) and from the Antillean islands (Dominica and Guadeloupe).
These groups were also recovered by the clustering analysis (Figure 3a) and by the BI (Figure 5) except for a few differences among the three analyses. In general, the results of these three analyses are linked as following: haplogroup 1 corresponds to cluster 1 and to the blue clade; haplogroup 2 corresponds to cluster 3 and to the yellow clade; haplogroup 3 corresponds to cluster 5 and to the red clade; and haplogroup 4 corresponds to cluster 4 and to the green clade. These colors were used to reference the groups in Figures 1, 3, and 5.
The BI analysis (Figure 4) shows that a few nodes of the tree presented posterior probabilities above 0.50 possibly due to the difficulties of the algorithm, implemented in BEAST, in solving datasets containing very similar sequences [60]. Despite this, the larger clades of the tree recovered most of the network haplogroups and the clusters of the clustering analysis as we discussed previously. These conclusions help to support the phylogenetic inferences raised by BI.
The cladogenetic relationship among the haplogroups may also be observed in the BI performed with the haplotypes for the ancestral area reconstruction. In general, the larger clades of this tree include haplotypes from the same haplogroup (Figure 5). This suggests a relationship among the haplotypes and the ancestral area reconstructed. We inferred dispersal events occurred based on the ancestral ranges and then proposed a dispersal route for the N. ephratae populations (Figure 7).
The ancestral populations distributed in MES and PAC (“1”; Figure 7a) suffered a temporary separation that split the populations of South America (southern PAC) from the Central American (northern PAC + MES) populations (“2”; Figure 7b). The ancestor of the haplogroups 1, 2, and 3 occurred in MES during this separation (“3”; Figure 7c). This leads to the origin of haplogroup 2 that was restricted to Central America (northern PAC + MES). After the reconnection of the N. ephratae populations (indicated with an asterisk in Figure 7), the Mesoamerican ancestor dispersed to the South American portion of PAC (“5”; Figure 7e) from where there was a new dispersion to SOU (“6”, Figure 7f) and to ANT-BOR (“7”, Figure 7f). This last dispersion gave rise to haplogroup 3 composed mainly of French Guiana and Antillean populations. The ancestral populations of the haplogroup 1 that had arisen in SOU were widely dispersed to CHA and PAR (Atlantic Forest lato sensu) reaching to BOR and southern PAC (“8”, Figure 7.f), but remaining limited to South America.
Still during the temporary separation between the South and Central America N. ephratae populations, the South American ancestor of the haplogroup 4 arose in southern PAC (“4”; Figure 7d). After the populations’ reconnection, there was a dispersal from southern PAC to northern PAC and to MES originating as the late ancestor of haplogroup 4 and clade H25-A59, which were restricted to these dominions (“5”; Figure 7e).
Based on the dispersal route proposed, we inferred that the N. ephratae populations currently distributed in Central America arose from distinct dispersal events. This becomes clearer when we also observe the geographic distribution of haplogroups 2 and 4 (Figure 1) whose occurrence areas overlap, but who have different genetic groups. That is, even though these haplogroups are distributed in the same area, the results are from different evolutionary events. The ancestors of these haplogroups possibly diverged during a temporary split among the ancestral N. ephratae populations from South and Central American between the late Pliocene and early Pleistocene (around 5.06 and 2.78 My; Figure 6 and 7). Haplogroup 2 originated from populations that were restricted to Central America during this split while the haplogroup 4 originated after the reconnection perhaps arising from populations that were dispersed from the northern South America (southern PAC) to Central America (northern PAC + MES). The same direction of dispersion (South America to Central America) was also identified for the ant species Neoponera villosa [16] dated from 0.46 to 0.28 My; earlier dispersions have been observed for N. ephratae.
Haplogroup 1 is exclusively South American and arose from a Mesoamerican ancestor that dispersed to the northern South America and then to the SOU dominion. Therefore, the dispersal events of the N. ephratae populations among Central and South America occurred in both senses, thus shaping the genetic and phylogeographic patterns observed here.
Although the causes responsible for the temporary split between the South and Central American N. ephratae population are not clear, some hypotheses can be raised. A population isolation caused by the geographic distance between the populations may have led to this split: The occurrence of isolation by distance was suggested by the Mantel’s test. The loss of distribution area could also have caused this effect in the populations. In this way, the reconnection of the populations could have occurred due to the demographic expansion that was detected by the neutrality tests and can be suggested from the star shape of the haplogroups in the haplotype network.
Geological and/or geographic factors could also help to explain this split. Following this approach, it is possible that the separation is related to the tertiary's tectonic and paleographical reorganization movements (in the late Pliocene), which led to the emergence of barriers and changes in dispersal routes in South America [61]. Specifically, this split between the N. ephratae populations may have been caused by the momentaneous effects of the elevation of the Panama isthmus. These effects would have lasted between 4.6 and 2.6 My [62] – a date close to the one estimated for the separations detected (5.06 to 2.78 My). Moreover, the early Quaternary climate changes (early Pleistocene) was characterized by temperature and dryness oscillations in the continents and was also impacted the adaptability and the migration of species and populations [61].
Although dispersion by water (inside flotsam carried by ocean or river currents or ferried by vessels) helps to explain dispersion patterns for termite species [54, 63, 64], it is more likely that the paths taken by the populations of N. ephratae were overland, which makes the PAC dominion an obligatory passage for the dispersal of populations between South and Central America. This dominion harbor peculiarly haplotypes from all haplogroups and present a higher value of intra-dominion genetic distance (Table 3.b). These features help to infer the PAC intermediary position for the N. ephratae dispersions.
The dispersal route traced for the populations analyzed here are very similar to the dispersal route observed for N. corniger species [15] including the dispersion from Central to South America, the eastward dispersion on South America lands, and the late occupation of Atlantic Forest. Nevertheless, the dispersal from South to Central America and the temporary split among the population have not yet been detected for N. corniger unlike N. ephratae. However, it is important to note that only the 16S mitochondrial marker was used for the N. corniger analyses [15]: The addition of other markers can lead to a more robust comparison among the population patterns of the two species.
Crews and Esposito [65] also studied dispersal routes and identified that South America is most probably the origin of most of the Caribbean arthropod fauna including N. ephratae, N. corniger, and other species of the genus Nasutitermes. These data contradict the inferences raised by Santos et al. [15] although there are some important methodological differences among the studies. Crews and Esposito [65] analyzed 18 species of Nasutitermes using the same mtDNA genes used here as molecular markers (16S + COII). The dispersal route proposed herein is according to the inferences made by Crews and Esposito [65] because the N. ephratae Antillean populations arose from a South America ancestor according to the ancestral area reconstruction. Specifically, for these island populations, it is possible that the dispersion from northern South America to the Antilles occurred via flotsam or floating wood carried to the islands by oceanic currents. This kind of dispersion has already been detected for Caribbean termite species [64]. Most of the islands’ termite species analyzed here are compose of haplogroup 3 together with samples from French Guiana (BOR dominion) and Trinidad and Tobago (located in the transition between BOR and PAC). This data suggests high genetic similarity among populations of these areas – although a larger sampling can better clarify the relations among the two dominions.
Regarding the BOR dominion, we found that the samples from this area have haplotypes typically found in other areas, which stayed distant from each other in the haplotype network (except for the haplotypes from French Guiana). This distance was also observed in BI. Some authors argue that the Amazonian biota (composed of parts of BOR and SOU dominions and the southeastern Amazon dominion – the latter was not sampled herein) have unnatural biogeographic origins [66, 67, 68, 69, 70]. As a speculation, this hypothesis could help clarify the existence of haplotypes that are genetically distant from each other in Amazonian localities that are geographically close to each other. However, we point that it is necessary to include a greater sampling of this region to clearly detect the genetic population patterns and the evolutionary events involving N. ephratae in this dominion.
In general, there are many questions to be explored about the phylogeographic processes of the South American and Neotropical species especially for termites whose studies are still new. Better sampling of N. ephratae and/or including new molecular markers in future studies, as well as addressing phylogeographic issues of other species, can help solidify the inferences made in this work and can expand the understanding of the evolutionary history of the group and of the neotropics.