There are more than 3,000 described species of termites (Blattaria: Isoptera) (Krishna et al., 2013), and they are one the most successful terrestrial organisms on Earth. This success is due to a sum of ecologically relevant characteristics. Termites may be poor land dispersers and weak flyers (< 2 km; Hu et al., 2007; Tonini et al., 2013; Mullins et al., 2015), however, transoceanic dispersal events mostly via flotsam (Bourguignon et al., 2017) facilitated their near world-wide distribution. Coupled with their ability to decompose cellulosic compounds (Jouquet et al., 2016), these eusocial insects thrived in new ecological opportunities and became keystone species and ecosystem engineers with an exceptional biomass in tropical ecosystems (Tuma et al., 2020).

Despite the extension of insects across the world and being the largest and most diverse group of organisms on Earth (Chapman, 2009), their phylogenetic relationships have only recently become evident due to advances in next generation technology (Trautwein et al., 2012; Yeates et al., 2016). Not much is known about the evolutionary history of species and their populations in South America, which is poorly studied in comparison to North America and Europe (see Tembrock et al., 2019), and what little is known mostly focuses on agricultural pests due to their economic relevance (i.e., Tembrock et al., 2019; Raszick et al., 2021; Vilardi et al., 2021). Phylogeographic studies would help to clarify the evolutionary history of species and their populations through the use of several resources such as population genetics and phylogenetics to understand the geological and climatic processes that determined the current geographic distribution of genealogical lineages (Avise et al., 1987; Avise, 2000).

Little is known about termites’ species evolutionary history, with most articles just focusing on species of the Rhinotermitidae family (e.g., Kutnik et al., 2004; Szalanski et al., 2004; Park et al., 2006; Tripodi et al., 2006; Jenkins et al., 2007; Luchetti et al., 2007; Austin et al., 2008; Lefebvre et al., 2008; Li et al., 2009; Scicchitano et al., 2018; Hyseni and Garrick, 2019). Although being the evolutionary most recent and most diverse family, the phylogeography of the Termitidae remains scarcely studied (Ozeki et al., 2007; de Faria Santos et al., 2017, 2022; Singham et al., 2017). Among these, there are only two phylogeographic studies for Nasutitermes species: N. corniger and N. ephratae (de Faria Santos et al., 2017, 2022), up to date. In both species, significant degrees of genetic structure were found across the different domains of the wide Neotropical region. However, another species of the genus, N. kemneri (Figure 1), seems to be restricted to a smaller area known as the South American “diagonal of open formations” (DOF) (Mathews, 1977).

The DOF consists of three phytogeographic domains (or biomes) crossing the South America from southwest to northeast: Caatinga, seasonally dry tropical forests of northeastern Brazil; Cerrado, the central Brazilian savanna; and Chaco, seasonally dry subtropical forests in northern Argentina, Paraguay, and Bolivia (Vanzolini, 1963; Zanella, 2010). All these domains are characterized by marked seasonality with severe droughts, nevertheless, they are different physionomically and in vegetation, harboring unique biotas well adapted to such fluctuations (Werneck, 2011).

The current geographical arrangement of the DOF is the result of the Central Brazilian Plateau uplift since the Cenozoic (Martins-Ferreira et al., 2020), changing the extensions of the Caatinga, Cerrado, and Chaco domains, in addition to the periods of climatic fluctuations during the Quaternary, when the Earth experienced ice ages interspersed with warmer periods. During these fluctuations, dry vegetation would have expanded and contracted, respectively, during glacial maximums and minimums, while tropical forest regions would have followed the reverse pattern (Pinheiro and Monteiro, 2010; Zanella, 2010). Specifically, dry regions were connected during glacial periods, while when interglacial periods occurred, tropical forests expanded and disconnected dry regions (Pennington et al., 2004; Werneck et al., 2012). Consequently, these events not only shaped the current geomorphology of the region, but became important drivers of the diversification and distribution of biota across the DOF.

Nasutitermes kemneri Snyder and Emerson, 1949 is a termite species endemic to the DOF that feeds on wood in early stages of decomposition in a diversity of vegetation types (Mathews, 1977). Due to its well delimited area of distribution, this termite is a relevant study model to understand how species respond to the vegetational, elevational, and climatic shifts between the Caatinga–Cerrado–Chaco phytogeographic domains. Accordingly, this study aims to determine if there are divergences between the populations distributed across the DOF, and if these form a spatial structure in genetic differentiation. Secondarily, through demographic and time estimation analyses, we tested if populations of N. kemneri contracted in size and were disconnected during interglacial periods but survived in climatic refugia. For this purpose, we used both nuclear (ITS) and mitochondrial molecular markers (COI, COII, and16S rRNA).

Sequencing efforts resulted in varying lengths: 321, 640, and 735 bp for 16 rRNA, COII, and ITS, respectively. ITS sequences were highly conserved in the 43 sampled individuals, with only 3 polymorphic sites. On the other hand, 14 polymorphic sites were present in the 55 16S rRNA sequences, 43 polymorphic sites were present in the 55 COII sequences, and the final alignment of 961 bp containing both 16S rRNA and COII genes revealed 53 polymorphic sites.

We used different approaches to estimate the degree of population differentiation. The genealogical relationships reconstructed by the ML phylogenetic tree evidenced high bootstrap support for the two mtDNA genetic groups of N. kemneri within the South American DOF, one at the Northeast (NE), and another one at the Southwest (SW; Figure 2).

The Median-Joining haplotype network of the two mitochondrial gene regions (16S rRNA + COII) also revealed the presence of these two groups, and that they are separated by 23 mutations, with no haplotypes in common (Figure 2). A total of 11 haplotypes were from the NE group, while 9 haplotypes were from the SW group. The ITS haplotype network illustrated one genetic group with only 4 closely related haplotypes (Supplementary Figure S1).

Furthermore, BAPS v6 and the correlated frequency model in GENELAND analyses used to infer population structure separated N. kemneri in the same two spatial clusters (Figure 3), further affirming the strong genetic differentiation found between these two genetic groups.

Notably, some genetic substructure in N. kemneri was hinted by the uncorrelated frequency model from GENELAND and BAPS non-spatial clustering of individuals. These results illustrated 4 clusters, distinguishing the Caatinga from northern Cerrado within the NE group, and separating the Cerrado from samples near the Paraná biogeographic dominion (Morrone et al., 2022; Supplementary Figure S2). Moreover, the phylogenetic tree shows high bootstrap values(= 97) which support the subtle divergence of the Caatinga and Cerrado within the Northeastern group (NE), while the easternmost individuals are the most divergent samples within the SW group (Figure 2). A 2-level AMOVA distinguishing both NE–SW and these four clusters further evidences the presence of sub-structure (Supplementary Table S5). However, since Cluster B (N = 5) and Cluster C (N = 3) are comprised of few individuals, further sampling is needed to firmly ascertain the presence of genetic sub-structure within Nasutitermes kemneri.

Considering the highest level of differentiation broadly supported by previous analyses, the AMOVA partitioned in two groups (NE–SW) revealed a high variation among groups (82.74%) contrasting with a low variation within groups (17.26%; Table 1). The high Fst (0.827, p < 0.0001) confirms a strong genetic population structure as supported by the phylogenetic tree, haplotype network and Bayesian spatial clustering. Regarding the genetic diversity across the South American DOF, haplotype and nucleotide diversity was overall high (Hd = 0.867, π = 0.00107), and high but slightly lower within both genetic groups (NE Hd = 0.767, π = 0.00734, SW Hd = 0.723, π = 0.00922).

Demographic inferences through Tajima’s D neutrality test showed no signal of population expansion or retraction for the entire species range (D = 0.395, p = 0.69) nor for each separate group (NE, D = −1.496, p = 0.13/SW, D = −1.883, p = 0.059). Finally, the dated phylogenetic analyses inferred that the divergence between NE and SW lineages occurred around 2.41 Mya (3.59–1.47, 95% height posterior density), during the Pliocene–Pleistocene boundary.

Results show a strong genetic and spatial structure within N. kemneri (Figure 2), evidenced by two highly differentiated groups, henceforth called populations. These populations are not restricted to the phytogeographic domain where they were sampled (Caatinga, Cerrado, and Chaco), instead, the NE population encompasses Caatinga and northern portions of Cerrado, while the Southwestern (SW) population ranges from central-south Cerrado to the northeast peripherical portions of the Chaco (Figure 2). This Northeast–Southwest pattern of spatial differentiation is shared with other taxonomic groups distributed along the DOF such as squamates (Werneck et al., 2012b; Recorder et al., 2014; Fonseca et al., 2018), amphibians (Oliveira et al., 2018), flies (Moraes et al., 2009), bees (Miranda et al., 2016), marsupials (Carvalho et al., 2011), trees (Resende-Moreira et al., 2017), and, to some extent, birds (Rocha et al., 2020).

Estimates suggest that these populations diverged ~ 2.5 Mya, during the Pliocene-Pleistocene boundary. During this period, the Cerrado plateaus were already uplifted to their present-day elevations while other areas within the Cerrado subsided, forming peripheral depressions (Del'Arco and Bezerra, 1989). The depressions inhibited gene flow between east and west of the plateau, and promoted diversification of subsequently isolated populations of flora and fauna. As discussed by Bonatelli et al. (2021), while several species across the DOF experienced range shifts, there is no single demographic response to Pleistocene climatic fluctuations. Some species experimented demographic expansion, others exhibited the opposite response, and some, like N. kemneri, underwent different processes since the separation of both populations (i.e., 59), where each experienced different evolutionary histories, illustrated by the different structures displayed by the haplotype network.

The NE population network has an intricate structure, with many “missing” haplotypes, which could be due to repetitive extinction-expansion events. Demographic inferences do not offer statistical support for this hypothesis, however, the geoclimatic history of northeastern Brazil does. During the Quaternary, major paleoenvironmental changes occurred during wet periods throughout the Pleistocene (Auler et al., 2004). Increased rainfalls facilitated the expansion of semi-deciduous forest to the Caatinga (Auler et al., 2004; Wang et al., 2004; Werneck, 2011), which possibly provided new ecological niches for N. kemneri to exploit besides xeric shrublands and an increase in biomass. On the contrary, arid periods led to regional extinction of taxa adapted to wet forests such as herpetofauna, as discussed by Gehara et al. (2017), and the cyclic nature of these humid-drier events left its mark on the NE population, accounting for the high genetic diversity and complex haplotype network. Moreover, recent studies have shown that termites are resistant to temperature shifts (Woon et al., 2019; Janowiecki et al., 2020; Woon et al., 2022), but precipitation would be the main factor in determining termite distribution, abundance, and survival rate (Woon et al., 2019; Pozo-Santiago et al., 2020). Although we did not measure this, it seems that N. kemneri is more abundant in the most humid phytogeographic domain, the Cerrado, than in the semiarid ones, Caatinga and Chaco. In the same way, elevation plays important roles in both restrictive and non-restrictive effects of temperature and humidity, as well as in the interaction between these factors, and while the Caatinga is composed mostly by flattened surfaces between 300 and 500 m above sea level (Silva et al., 2017), the Cerrado has higher plateaus in central and northeastern regions (Werneck et al., 2012).

The SW population displays a star-shaped haplotype group, which is characteristic of refugia with later population expansion, yet demographic inferences lack statistical support. Nevertheless, during the Quaternary climatic and vegetation fluctuations, the general area currently occupied by Cerrado appears to have been drier (Collevatti et al., 2020; Oliveira et al., 2020), paleoclimatic model reconstructions predicted areas of historical climatic stability (Werneck et al., 2012). These stable areas (the higher plateaus in central and north-eastern Cerrado) were important refugia for taxa such as trees (Ramos et al., 2007; Buzatti et al., 2017; Resende-Moreira et al., 2017; Camps et al., 2018), and may, consequently, also serve as refugia for wood feeding termites.

Finally, while the spatially explicit Bayesian clustering programs BAPS and GENELAND reached a consensus in dividing N. kemneri in two populations, the uncorrelated frequency model in GENELAND and the non-spatial BAPS clustering are also worth noting and briefly discussed. These less robust approaches hint at subtle population substructure within both the NE and SW populations, wholly distinguishing Caatinga (Cluster A) and northeastern Cerrado (Cluster B) in NE, and separating central western Cerrado (Cluster C) from the rest of Cerrado (Cluster D) in SW. Tree species with a similar distribution across the DOF (and some that extent to the Parana dominion), show fairly recent (~ 800.000 years BP) fine-scale genetic structure within the Chaco-Cerrado-Atlantic Forest, from 3 to up to 7 haplogroups/lineages (Collevatti et al., 2009; Novaes et al., 2010, 2013). These lineages present a similar pattern of distribution to Clusters A, B, C, and D of N. kemneri, and the recent divergence time between the main populations (NE/SW) at the Pliocene-Pleistocene transition (with confidence intervals spanning from Late Miocene to Middle Pleistocene) may explain the low resolution of these analyses, where subtle climatic or geographic differences within these phytogeographic domain may be driving this incipient or shallower population sub-structuring. In fact, biogeographic districts based on tree species (Françoso et al., 2020) partially match the distribution of these clusters, citing different temperatures, seasonal variation, radiation and the degree of cover transformation as responsible for the floristic differences. However, further sampling is needed to determine the extent and significance of these patterns in N. kemneri.

This study provides important inferences about the general panorama of the evolutionary history of N. kemneri in the DOF, its endemic region, offering important contributions to the understanding of biogeographic and phylogeographic issues in the Neotropics. Our results support an ancient divergence between two well-defined lineages (NE and SW) for N. kemneri, which occurred during the Plio-Pleistocene transition, reinforcing the role of climatic and vegetation variation events in the diversification of the DOF biota (Del'Arco and Bezerra, 1989; Wang et al., 2004). Further climatic fluctuations during the Quaternary had pronounced effects on the Caatinga biota (Bonatelli et al., 2021) and could have facilitated the expansion and diversification of N. kemneri across the NE region of the DOF, while refugia zones in the Southeastern region may account for the constricted diversity patterns within the Cerrado, and altogether explain the unique evolutionary history of these two separate populations. The findings presented in this work, however, came from a few mitochondrial markers, and the investigation using multi-locus nuclear data can reveal other relevant patterns in the population genetic structure in this species.

Most of the integrative phylogeographic and predictive studies in the DOF rely heavily on data from herpetofauna and plants, where insects are less than 15% of the species considered (Del'Arco and Bezerra, 1989). In this sense, data presented in this study not only aids the understanding of the evolutionary history of N. kemneri and add evidence to the response patterns described for this area (Del'Arco and Bezerra, 1989; Wang et al., 2004), but also partly accounts for our knowledge gap. Furthering our knowledge on biotic variables that influenced the demography and distribution of DOF endemic species is key, and especially insects (Wilson and Fox, 2021) for future decisions in conservation strategies of the DOF threatened phytogeographic domains and potential effects linked to climate warming.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary material.

KK and RGS made the laboratory work. KK, VP-O, RGS, and TC did the analyses. AV, EC, and RHS provided the specimens. RHS provided photos of the N. kemneri soldier. KK, VP-O, and TC did the writing. RHS, AV, and EC edited the manuscript. All authors contributed to the article and approved the submitted version.

This research was supported by Conselho Nacional de desenvolvimento Científico e Tecnológico (CNPq): productivity in research fellowship proc. #309820/2020-0 to AV and #308408/2019-5 to EC, and

São Paulo Research Foundation, Brazil (FAPESP): grants #13/20068-9 to EC and #13/03767-0 and #20/06041-4 to TC.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.