Cryptic subterranean diversity: regional phylogeography of the sand termite Psammotermes allocerus Silvestri, 1908 in the wider Namib region

Psammotermes allocerus Silvestri, 1908 is the only described species representing the genus Psammotermes Desneux, 1902 in Southern Africa. The large geographical range of this subterranean termite covers both summer and winter rainfall regimes. Deadwood is the preferred food when available, but in more arid habitats, both live and dead grasses form the major dietary component. Along the Namib Desert margins, the species’ localised herbivory creates circular bare patches known as fairy circles. For a regional phylogeographic study of this species, we sampled 65 sand termite populations within drier parts of Namibia, South Africa, and Angola. Based on combined molecular and ecological data, we found considerable genetic diversification within P. allocerus. Analyses of two mitochondrial markers (COI, COII), including a Bayesian inference tree, haplotype analysis and genetic distances suggest a delineation into seven highly differentiated genetic groups. The ‘Succulent Karoo’ group is additionally characterised by unique features of the royal chamber, nest and tunnel system. In conclusion, our data suggest that P. allocerus should be not regarded as one species but as a species complex. Termites of each analysed group ‘Northern Namib’, ‘Western Kalahari Basin’, ‘Nama’, ‘Southwestern Kalahari’, ‘East Gariep’, ‘Southern Namib’ and ‘Succulent Karoo’ should be considered as distinct species. The species name P. allocerus should be used for termites of the ‘Succulent Karoo’.


Introduction
More than 3000 termite species in 281 genera and seven families are presently distinguished (Chouvenc et al., 2021;Engel & Krishna, 2004;Inward et al., 2007a;Kambhampati & Eggleton, 2000;Krishna et al., 2013), of which 165 species of 54 genera in five families have been documented for southern Africa (Uys, 2002). Among these, the sand termites (Psammotermes Desneux, 1902;Roonwal & Bose, 1962) have developed extraordinary plasticity with regard to habitat tolerance (Coaton, 1958(Coaton, , 1981Coaton & Sheasby, 1972, 1973, 1974, for both summer and winter rainfall regions as well as humid and arid climates. Zeidler (1997) identified Psammotermes as one of only four genera (also including Microhodotermes, Baucaliotermes, and Hodotermes) that tolerate the hyperarid habitats with only 30 mm MAP. Grube (1993) and Grube and Rudolph (1995) studied the water relations of Psammotermes in a laboratory experiment and reported on an extraordinary ability to take up capillary soil humidity using hypopharynx hairs. Coaton and Sheasby (1973) and Zeidler (1997) described Psammotermes as generalists that feed on wood and litter. Zeidler (1997) assumed that Psammotermes in the Central Namib Desert feed solely on subterranean organic material. Juergens (2013) confirmed this, and reported that Psammotermes fed on the roots of live grasses. Along the Namib Desert margins, the species' localised herbivory in grassland creates circular bare patches known as fairy circles (Juergens, 2013;Juergens et al., 2015;Jürgens et al., 2020).
In Africa, the core range of Psammotermes is in the western parts of southern Africa (Coaton & Sheasby, 1973). The wider range also includes localities in the eastern parts of South Africa, Zimbabwe, and Mozambique. Despite the large geographical range, only one species, Psammotermes allocerus Silvestri, 1908, is recognised from southern Africa.
The regional genetic variability of P. allocerus in southern Africa has not yet been studied. Sheasby (1972, 1973) published a detailed distribution map based on roughly 800 documented localities and considered all P. allocerus collections as belonging to one species. Similarly, only a few samples of P. allocerus from South Africa were integrated into continental-wide molecular phylogenetic studies (Inward et al., 2007a, b;Schyra et al., 2018). Samples from the type locality (Lüderitz, Namibia) have not yet been used for morphological (Grube & Rudolph, 1995) or genetic studies.
The remaining five species of Psammotermes Desneux, 1902(Roonwal & Bose, 1962 occur in northern and eastern Africa, Madagascar and India. Previous studies focused on morphological analyses of Psammotermes hybostoma, e.g., chemistry and anatomy of the frontal gland in soldiers (Krasulová et al., 2012), sex and trail pheromone (Sillam-Dussès et al., 2011), developmental pathways , or the foraging activity in Egypt (Salman et al., 1988). So far, little genetic work has been done on this species (Inward et al., 2007b;Jürgens et al., 2020), and for P. voeltzkowi Wasmann (1910), a few sequences were included by Inward et al. (2007a) in a broader phylogeny.
The recent study by Schyra et al. (2018) and others (Fuller, 1921;Hausberger et al., 2011) showed that a robust species delimitation on species level based only on morphological characters sometimes remains difficult, also due to variation in soldier size, as in P. allocerus (Uys, 2002) and P. hybostoma . Therefore, applying new molecular methods such as sequence-based species delimitation (Fujita et al., 2012;Pons et al., 2006) and DNA barcoding methods (Hebert & Gregory, 2005) has become more relevant for species identification of termites. However, this has led to an increase in the number of newly discovered and putative or cryptic species (Austin et al., 2007;Bourguignon et al., 2013;Hausberger et al., 2011;Roy et al., 2006). The two mitochondrial cytochrome oxidase subunits COI and COII are the most suitable markers to identify species due to easy, cheap processing, high substitution rate, occurrence in high copy numbers, and their role as DNA barcode (COI) (Brown et al., 1979;Cheng et al., 2014;Hebert & Gregory, 2005;Hebert et al., 2003). Both markers can distinguish species, in contrast to certain nuclear markers (Hausberger et al., 2011;Huang et al., 2013;Schyra et al., 2018), and they can also reveal cryptic or putative species (Bickford et al., 2007;Hausberger et al., 2011;Roy et al., 2006).
We aim to clarify the evolutionary status of P. allocerus populations from its core range within arid regions of Namibia, South Africa and Angola, including specimens from the type locality Lüderitz (Namibia). Additionally, we include data from nest structures (royal chamber, nest, and tunnel system) and compare them with the genetic results.

Material and methods
In total, Psammotermes allocerus termites were collected (mainly in the early morning hours) from 65 study sites comprising 50 Namibian, 14 South African, and one Angolan study site between 2011 and 2020 (Fig. 2, Table S1). Samples were collected from subterranean nests or underneath sand galleries around foraged grass and dead wood, or directly from the top layer of sand. The material was directly preserved in 75% ethanol. The sampling includes 29 sites with fairy circles and 36 without obvious fairy circles. For details regarding the collected samples and GenBank accession numbers, please refer to Table S1.

DNA extraction, PCR amplification, and sequencing
The DNA was extracted according to the protocol 'hair, nails, and feathers' from the 'E.Z.N.A.® Forensic DNA Kit' (Omega Bio-Tek, Inc. Norcross, USA) and as described in Jürgens et al. (2020). Each termite was dried and homogenised with a ball mill MM 400 (Retsch GmbH, Germany) at 30 Hz for 60 s. The final DNA elution step was repeated with 50 µl elution buffer for each step. Afterwards, two mitochondrial markers were sequenced and analysed: cytochrome oxidase I (COI) and cytochrome oxidase II (COII). The first marker was partially amplified with the primers C1F0 and C1R1  under the PCR conditions described by Legendre et al. (2008) for the gene COI. A second reverse primer was self-designed C1R1PsR (5′-TAG TAT GTG TCG TGT AGT AC-3′) and used under the same PCR conditions as C1R1. The PCR mixture described by Lo et al. (2000) was adjusted for a final volume of 10 µl, including 1 µl 10 × buffer, 2 mM MgCl2, 0.2 mM dNTP, 1.25 u Taq polymerase, 1 µM of each primer, and 10 to 20 ng DNA. The second marker, COII was amplified with the primer pair ModA-tLeu (Miura et al., 1998) and B-tLys (Beckenbach et al., 1993) under the PCR conditions listed by Legendre et al. (2008) but adjusted to a hot start Taq polymerase. The first denaturation step was elongated at 94 °C for 14 min. The PCR mixture for COII included 1 µl 10 × buffer, 0.1 mg/ ml BSA, 2 mM MgCl2, 0.2 mM dNTP, 1.25 u hot start Taq polymerase, 1 µM of each primer, and 20 ng DNA. Each PCR product was purified with FastAP and exonuclease I (Thermo Fisher Scientific, Waltham, USA) and sequenced to the conditions of the 3500 Genetic Analyser (Thermo Fisher Scientific, USA) as described in the manufacture instructions. Each PCR product was sequenced forward and reverse.

Phylogenetic analyses
The sequences were visualised and manually edited with FinchTV version 1.4.0 (2006 Geospiza Inc.). The P. allocerus consensus sequences were aligned with sequences from three Rhinotermitidae species, Termitogeton planus (KP026298, Bourguignon et al., 2015), Prorhinotermes canalifrons (KP026256, Bourguignon et al., 2015), and Schedorhinotermes breinli (JX144935, Cameron et al., 2012), downloaded from NCBI GenBank. All sequences were aligned under the default settings of the ClustalW algorithm in BioEdit version 7.0.5.3 (Hall, 1999). Both partly sequenced mitochondrial markers (COI 362 base pairs (bp) and COII 448 bp) were concatenated (810 bp). Missing sequences in the concatenated alignment were filled with N's as for COI sequences of P. allocerus collection collected in Akadispass (AP, South Afrika) and Keerweder (KW, Namibia). The Bayesian inference (BI) was conducted with MrBayes version 3.2.6 (Ronquist & Huelsenbeck, 2003) according to the instructions by Rohwer et al. (2014). Model selection was performed for the alignments using default settings in MEGAX version 10.0.4 (Kumar et al., 2018). The model Hasegawa-Kishino-Yano (HKY) with a discrete gamma distribution (+ G) and a proportion of invariable sites (I) was selected based on the Akaike information criterion (AIC) (Akaike, 1974) and the Bayesian information criterion (BIC) (Schwarz, 1978) for both alignments. The concatenated alignment was not partitioned because phylogenetic results, topology, and support values are similar to partitioned analyses. From five million generations and two simultaneous runs for four Metropolis-coupled Monte Carlo Markov chains (MCMCMC), one tree was saved from 100 generations. The number of chains was set to one 'cold' and three heated chains. The 'burnin' was selected for the first 25% of all trees. Paup version 4.0a (Swofford, 2001) was used to calculate the posterior probability values for a majority-rule consensus tree. The heuristic search and the tree bisection-reconnection (TBR) were implemented.
The sequence replicates were set to 100. The retaining of all minimum-length trees and the collapse of zero-length branches were chosen. Maximally 1000 trees were saved per replicate. The number of branch exchanges was set to 5 million per bootstrap replicates. Between-group and withingroup mean p-distances were calculated in MEGAX version 10.0.4 (Kumar et al., 2018) for the concatenated alignment.

Haplotype analyses
The median-joining network (Bandelt et al., 1999) was performed for the concatenated sequences of COI and COII P. allocerus sequences to visualise shared haplotypes using POPART (Leigh & Bryant, 2015). The genetic diversity indices were calculated for each group in DnaSP v6 (Rozas et al., 2017): number of haplotypes (NH), number of sequences (NS), number of segregating sites (S), nucleotide diversity (π), and haplotype diversity (Hd). The Tajima's D-value (Tajima, 1989) was calculated for each genetic group. This value indicated the natural selection as well as the type of selection pressure. A negative value indicated an expansion in population size, whereas a positive value indicated a recent population bottleneck. The null hypothesis assumed that the gene region of interest is neutral. If the test is significant (p < 0.05), the null hypothesis of neutrality will be rejected. All parameters were calculated for the concatenated alignment as well for the COI and COII sequences shown in the appendix. Assignment of the sequences to the haplotypes is shown in Supplementary Information (Table S2).

Phylogeny analyses
The mitochondrial sequences were visually examined according to the presence of NUMTs (nuclear mitochondrial sequences). Characteristics such as stop codons, frameshift mutations and insertions-deletions were not detected. Finally, the COI marker (362 bp) contained 137 constant, 70 uninformative, and 155 parsimony-informative positions. For the second marker, COII (448 bp), 259 positions were constant, 89 were uninformative, and 124 were parsimonyinformative. The concatenated alignment (810 bp) resulted in 507 constant positions. A low number of 92 positions were uninformative, and 211 were parsimony-informative.
All mitochondrial COI and COII sequences were concatenated in one alignment, and the resulting Bayesian Inference (BI) analysis shows a strongly supported differentiation of the 65 P. allocerus samples from the outgroup (1 posterior probabilities (pp)). We identified seven genetic groups which may represent putative species clusters according to the high posterior probability values of the nodes (pp > 0.95), higher values of p-distances between these clusters as within the clusters as well haplotype analyses according to Roy et al. (2006) and Hausberger et al. (2011). The seven strongly supported clades (0.99 pp, Fig. 1), roughly fit a differentiation from south to north within the study area (Fig. 2).
The most basal clade within Psammotermes (dark blue, Fig. 1) was the 'Succulent Karoo' clade, as all seven samples are located within the Succulent Karoo biome in the north-western part of the winter rainfall region of South Africa (Fig. 2). It included the collection of the topotype near Lüderitz (LUE, Namibia). The next clearly distinguishable genetic group was the 'Southern Namib' clade, with the majority of the samples located at the eastern margin of Southern Namib (Fig. 2). Three localities are very well supported (> 0.99 pp), and only two are sisters to one other (RO and DV). The next higher positioned clade was split into two larger subclades at an equal level of differentiation, one occurring in south-eastern Namibia, adjacent to South Africa. The other clade includes samples from northern Namibia and Angola. The well-supported subclade ( Fig. 1) includes collections from south-eastern Namibia and the eastern Richtersveld (South Africa). Phylogenetically it is subdivided into three genetic groups: 'East Gariep', 'Southwestern Kalahari' and 'Nama'. The 'East Gariep' group is represented by 13 sites, of which six are located in the northeastern Richtersveld and the adjacent larger Orange River valley in Namibia (Fig. 2). The weakly supported (0.52 pp) combined groups 'Nama' and 'Southwestern Kalahari' are sisters to the 'East Gariep' group. Both individual groups (Southwestern Kalahari and Nama) are very well supported (1 pp). The 'Southwestern Kalahari' clade comprises eight collection sites situated in dunes and sand valleys east of the Karas mountains. The 'Nama' group (yellow) includes six study sites in the southern Namaland of the Namibian Karas Region. The second larger subclade is well supported (0.92 pp) and included the groups 'Western Kalahari Basin' and 'Northern Namib', both very well supported (1 pp). The 'Western Kalahari Basin' comprises 11 collection sites ordered in a polytomy. This group includes sites that are located from south-eastern Namibia (SK3, SA) to the Waterberg (RI and RV) and the north-easternmost site Katima Mulilo (KM). The 'Northern Namib' group encompassed all remaining study sites, located from the Rössing mountain (ROE) to the Iona National Park in Angola (IO).
Overall, the mean p-distances are higher between the seven groups than within each group (Tables 1 and 2). The mean p-distance between the groups ranges from 0.0508 ('Nama' vs. 'East Gariep') to 0.1007 ('Kalahari Basin' vs. 'Succulent Karoo'). The values of the mean p-distance are lower and ranged from 0.0100 ('Succulent Karoo') to 0.0363 ('Northern Namib').

Haplotype analyses
The median-joining network based on the concatenated alignment comprised 31 observed haplotypes, including 18 unique haplotypes (which includes a single sample) and 12 calculated hypothetical haplotypes (Fig. 3). Eight hypothetical haplotypes are situated in the centre of the median-joining network. The observed haplotypes of the seven genetic groups are linked over this centre. The highest value of unique haplotypes was calculated for the 'Southern Namib' group with six haplotypes. In contrast, the 'Succulent Karoo' comprised two haplotypes calculated for six sequences. The two groups, 'Western Kalahari Basin' and 'Succulent Karoo', are connected to the centre with a high number of mutation steps (8 and 10, Fig. 3). The number of haplotypes differed slightly between the median-joining network analysis (Fig. 3) and the calculated number of haplotypes gained by DnaSP v6 (Table 3). The nucleotide diversity (π) ranged from 0.00416 ('Western Kalahari Basin') to 0.01597 ('Northern Namib', Table 3). The lowest haplotype diversity was calculated for the 'Succulent Karoo' with 0.78571 and the highest for the 'Southern Namib' with 1.0000 (Table 3). The number of segregating sites (S) was highest within the 'Southern Namib' group with 22 S and lowest within the 'Succulent Karoo' and 'Western Kalahari Basin' with each six S ( Table 3). The genetic groups undergo a different selection pressure. A positive Tajima's D value was calculated for the groups: 'Northern Namib', 'Western Kalahari Basin' and 'Southwestern Kalahari' and a negative Tajima's D value was estimated for the remaining groups ('Nama', 'Southern Namib', 'East Gariep', 'Succulent Karoo', Table 3).
Two additional features distinguish the 'Succulent Karoo' group from the other groups. The tapetum of the tunnels within the nest system of P. allocerus colonies in the 'Succulent Karoo' is whitish. In all other groups, it is blackish (Fig. 4). Only in the 'Succulent Karoo' nests was it possible to locate the royal chambers with paired kings and queens at a depth of 15 to 30 cm below the soil surface (Fig. 4). In all other groups, one of the authors (NJ) never locate royal chambers or kings or queens in the first 100 cm below the soil surface, despite some 100 nest excavations.

Discussion
Our molecular results show that the studied populations of the sand termite Psammotermes allocerus can be subdivided into seven putative species clusters (Fig. 1), following the definition of Bickford et al. (2007). This differentiation is confirmed by the Bayesian Inference analysis of the COI and COII markers (Fig. 1), high genetic p-distances between the seven clusters and lower p-distance within the clusters (Tables 1 and 2) and observations regarding the termite  Fig. 2 nests, colour of the tapetum, and position of the royal couple. A similar genetic differentiation into putative species clusters was investigated by Hausberger et al. (2011) for Trinervitermis and Odontotermes species from Benin and by Roy et al. (2006) for Cubitermes species from Gabon. In comparison to the within-cluster p-distance values  Our results based on mitochondrial DNA markers alone clearly allow the detection of robust putative species clusters, as has been previously shown in several other studies (Hebert et al., 2003;Hebert & Gregory, 2005;Hausberger et al., 2011;Bourguignon et al., 2013Bourguignon et al., , 2015Schyra et al., 2018;Korb et al., 2019). We are aware that results based on mtDNA markers alone can, in some cases, lead to misinterpretation due to an overlap of sequence similarities (Meier et al., 2006;Rubinoff et al., 2006). We could not identify this problem in our alignments, and we did not find nuclear mitochondrial sequences (NUMTs) that might result in misleading phylogenies (Sorenson & Quinn, 1998).
The putative species clusters are also supported by the high numbers of mtDNA haplotypes (Fig. 3) and haplotype diversity values (Table 3) between the clusters. These results also suggest the subdivision of P. allocerus into seven putative species clusters. Similar values and haplotype differentiation were investigated by Austin et al. (2004) for populations of Reticulitermes flavipes in the USA and by Nobre et al. (2006) for Reticulitermes species in Portugal. In contrast to these studies, we calculated an overall higher nucleotide diversity (π) for each cluster. Likewise, the high numbers of hypothetical haplotypes, unique haplotypes, and mutation steps (Fig. 3) support the hypothesis that the putative species clusters are not recently evolved and have low gene flow between them, as the presence of unique haplotypes points to a limited current level of genetic exchange (Nobre et al., 2006). The differentiation into seven putative species is underlined by Tajima's D-value results. We analysed a different type of selection pressure for each cluster, shown by a positive or negative Tajima's D-value (Table 3). This was also shown for R. flavipes populations by Austin et al. (2004). These results support the hypothesis that P. allocerus should no longer be considered one species. The phylogenetic results and ecological differences show that P. allocerus should be regarded as a species complex composed of seven species. The species name P. allocerus is tied to collections of the 'Succulent Karoo' clade due to the location of the type specimen near Lüderitz (Namibia). The molecular results strongly support a further differentiation of the remaining analysed samples into the six species: 'Northern Namib', 'Southern Namib', 'Southwestern Kalahari', 'Nama', 'East Gariep' and 'Western Kalahari Basin'. Formal taxonomic decisions can only be made with a comparison of morphological data from soldiers of all putative species.
Most clearly, the well-known distribution ranges of the defined clades support the hypothesis that Psammotermes shared the evolutionary history of many other austro-temperate taxa (Bakkes et al., 2018;Colville et al., 2015;Davis, 1994;Jürgens, 1991;Jürgens et al., 1997;Naskrecki & Bazelet, 2012;Pitzalis & Bologna, 2010;Pitzalis et al., 2016) which adapted to arid conditions and summer rainfall seasonality during the Neogene. The distribution ranges of the defined phylogenetic clades of Psammotermes are paralleled by the range types of plant species (Jürgens, 1991;Jürgens et al., 1997). The most basal group ('Succulent Karoo') in the Psammotermes tree falls within the northern parts of the Succulent Karoo biome in the winter rainfall zone of Southern Africa. This region is characterised by localised plant taxa (Brownanthus marlothii, Dregeochloa pumila, Stipagrostis dregeana, Cephalophyllum ebracteatum, Stipagrostis geminifolia, Dracophilus, Juttadinteria, and Psammophora, Jürgens, 1991;Jürgens et al., 1997). The 'Southern Namib' group covers a neighbouring geographical space. However, this shift in geographical space requires a very important evolutionary step, the adaptation to the environmental conditions of the tropical summer rainfall zone. This range is shared by the various plant taxa (Hermannia minifolia, Monsonia ignorata, Sesamum abbreviatum, Table 3 Results of the genetic diversity indices calculated for the concatenated alignment of COI and COII markers and seven clades within P. allocerus: number of haplotypes (NH), number of sequences (NS), number of segregating sites (S), nucleotide diversity (π), and haplotype diversity (Hd), Tajima Basin' group region has not been explored completely and is therefore not compared here with known area types. In summary, the ranges of almost all Psammotermes clades match the distribution patterns of various plants within the arid regions, suggesting a shared response to underlying physico-geographic factors. The observed spatial pattern suggests an important role in allopatric speciation. In fact, the reproductive system of Psammotermes, in combination with the clear preference for sandy habitats, should drive allopatric speciation. Regarding the reproductive (and dispersive) phase, Coaton and Sheasby (1973) documented a long alate season from March to December for P. allocerus. As with other termite species (Nutting, 1969), alates of P. allocerus fly mainly during the rainy season and after the first rainfall events (Juergens, 2013;Juergens et al., 2015), but these rainfall events are scarce, irregular, and localised within western southern Africa (Coaton & Sheasby, 1973) and especially in the Namib Desert (Eckardt et al., 2013;Henschel et al., 2005). During the last two decades, all swarming events observed by us support the assumption that dispersal distances are limited to several hundred metres or a few kilometres at most. This would suggest that alates do not disperse further than the freshly moistened sand allows after the local rainfall events. This spatial limitation of gene flow may be important because, under hyper-arid conditions, the survival of Psammotermes is often limited to spatially isolated sandy habitats or woody vegetation along dry river beds. Therefore, we propose the hypothesis that a lack of genetic exchange between the ranges of species resulted from geographic distances too large for gene transport by alates.
To resolve the relationship between the putative taxa of the two larger groups (group one including 'East Gariep', 'Nama', and 'Southwestern Kalahari' and group two with 'Northern Namib' and 'Southern Kalahari', Fig. 1), the analysis of additional nuclear markers, e.g. ultraconserved elements (UCEs) (Buček et al., 2022;Hellemans et al., 2022) or the method double digest restriction-site associated DNA Sequencing (ddRadseq) (Blumenfeld et al., 2021) may help to investigate the time since speciation. The current study includes collections from more arid regions of southern Africa. In the future, an extension of the study area to disjunct populations in more humid areas within South Africa (Eastern Cape, Western Cape, St. Lucia) should take place. Morphological analyses of soldiers from all investigated seven species and an extension to more eastern taxa and additional Psammotermes species, e.g., P. hybostoma, will hopefully allow formal taxonomic determinations.
Acknowledgements In Namibia, we thank the Ministry of Environment, Forestry, and Tourism (MEFT), the National Commission on Research Science and Technology (NCRST), the National Museum, and the National Botanical Research Institute (NBRI) for research permission. We thank the National Museum for the curation and loan of specimens. Our thanks also go to the Gondwana Collection, Mannfred Goldbeck, and the management of Namib Desert Lodge for the generous support of research activities at Dieprivier over many years. We are also grateful that Namibrand Nature Reserve (NRNR) and numerous conservancies gave permission to work in their area. In Angola, we thank Fernanda Lages at Instituto Superior de Ciências de Educação (ISCED, Lubango) for help in logistic and administrative matters. Bruce Bennett participated in and logistically arranged fieldwork in remote parts of the Iona National Park. In South Africa, we thank the Northern Cape Province Department of Environment and Nature Conservation (DENC), South African National Parks (SANParks), and especially Brent Whittington and Pieter Van Wyk for permission and logistic support. We thank Martin Husemann (Centrum für Naturkunde, CeNak) for helpful comments on an earlier draft of this manuscript. Jürgen Deckert (Museum für Naturkunde Berlin) for the loan of type specimen of P. allocerus. We sincerely thank Barbara Rudolph for years of support and advice on this project, for initial molecular analyses and general interpretations of the data, and for managing the molecular laboratory. We gratefully thank Imke Oncken for her intensive and great support in the preparation of and during several field surveys and laboratory work. Vilho Mtuleni, Ingo Homburg, and several others helped during numerous field surveys. Clärin Bohn carried out preliminary molecular studies on Psammotermes.
Author contributions FG performed the molecular and morphometric work and drafted the text. All authors contributed to fieldwork and collecting, the bulk fieldwork was done by FG and NJ. FG and JO did the statistical analyses. All authors contributed to the interpretation of the results and improved the text. NJ conceptualized the project.
Funding Open Access funding enabled and organized by Projekt DEAL. The research was carried out in the framework of Southern African Science Service Centre for Climate Change and Adaptive Land Management (SASSCAL) and was sponsored by the German Federal Ministry of Education and Research (BMBF) under promotion number 01 LG 1201 N.

Data availability
The DNA sequences will be available on GenBank. All other datasets generated and/or analysed during the current study are available from the corresponding author on reasonable request.

Declarations
Ethics approval NCRST approved our research and granted the permit number RPIV00032018.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.