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, 1981; Coaton & 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. Coaton and 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 (Bourguignon et al., 2012), 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 (Bourguignon et al., 2012). 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 (Lo et al., 2004) under the PCR conditions described by Legendre et al. (2008) for the gene COI. A second reverse primer was self-designed C1R1PsR (5′-TAGTATGTGTCGTGTAGTAC-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 within-group 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).

Results

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 parsimony-informative. 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).

Fig. 1
figure 1

Phylogeny of 65 P. allocerus collections inferred by the Bayesian analysis of COI and COII markers. Support values are given in posterior probability. Clades are coloured according to the genetic group. Dark blue: Succulent Karoo; Light Blue: Southern Namib; Dark green: East Gariep; Light Green: Southwestern Kalahari; Yellow: Nama; Ochre: Western Kalahari Basin; Red: Northern Namib. Abbreviations of study sites are shown in Fig. 2

Full size image
Fig. 2
figure 2

Distribution map of all Psammotermes allocerus collections and their assignment to the genetic groups found in the phylogenetic analyses (colours). Circles: Collections used for the phylogeny. Black lines: Country borders. Northern Namib: (AH) Aba Huab rivier, (GV) Giribesvlakte, (HD) Hartmann Dunes, (HO) Hoada, (IO) Iona, (MF) Marienfluss, (PR) Purros, (ROE) Rössing mountain, (SO) Sorris Sorris, (TM) Tomakas; Western Kalahari Basin: (AHV) Alt-Hartebeestvlei, (AM) Alex Muranda, (GU) Gariganus, (KM) Katima Mulilo, (RI) Rimini, (RS) Rundu South, (RV) Ravenna, (SA) Samehaling, (SK) Swartkop, (WF2) Warmfontein; Nama: (BY) Barby, (FRC) Fish river Canyon, (GO) Goageb, (SH) Seeheim, (WV) Witvley; Southwestern Kalahari: (GM) Goedmoed, (GR) Gurus, (KFE) Kalkfontein East, (NK) Neikop, (SKN) Swartkop North, (TT) Tranental, (UK) Ukamas, (WF1) Warmfontein; East Gariep: (AP) Akadispass, (BD) Belda, (BH) Bruinheuwel, (DH) De Hoop, (KFE) Kalfontein East, (KFW), (KV) Koeroegabvlakte, (NO) Norachas, (NU) Numees, (RB) Rooiberg, (TB) Tatasberg, (TBQ) Tatasberg Quarzfield; Southern Namib: (DV) Dieprivier, (GA) Garub, (KH) Keetmanshoop, (KW) Keerweder, (RO) Rostock, (SR) Sesriem; Succulent Karoo: (HN) Holhat North, (KTV) Kortdoornvlakte, (LUE) Lüderitz, (RDE) Red Dune East, (SP) Springklipplain, (YD) Yellow Dune, (YDRD) Yellow Dune Red Dune Transect

Full size image

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’).

Table 1 Estimates of evolutionary divergence for the concatenated alignment (COI and COII) for mean p-distances (± SE) between the groups
Full size table
Table 2 Estimates of evolutionary divergence for the concatenated alignment (COI and COII) for mean p-distances (± SE) within the groups
Full size table

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).

Fig. 3
figure 3

Median-joining network of the combined COI and COII P. allocerus sequences. Coloured circles represent the observed haplotypes, and the size is proportional to the number of collections. Black circles represent missing haplotypes. Marks show the number of mutation steps. Haplotype numbers are gained from DnaSP v6. Dotted lines and colours mark haplotypes according to the genetic group of the phylogeny

Full size image
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’s D value (nonsignificant = ns, for p > 0.10)
Full size table

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.

Fig. 4
figure 4

Differences in the tapetum colour of three P. allocerus colonies and the royal pair. A Whitish tapetum of a nest from the Springklipplain from the ‘Succulent Karoo’ group (South Africa, 26 September 2016). B First image of the king and queen of P. allocerus from Yellow Dune (‘Succulent Karoo’, South Africa, 08 March 2015). C Blackish tapetum of a nest from Dieprivier (‘Southern Namib’, Namibia, 05 April 2017). D Blackish tapetum and chambers filled with foraged grass from Iona (‘Northern Namib’, Angola, 26 September 2016). Images taken by Norbert Jürgens, Felicitas Gunter

Full size image

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 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 calculated by Hausberger et al. (2011), p-distance values within the putative species clusters of P. allocerus were higher (Table 2), indicating a higher genetic diversity within at least some putative species clusters of P. allocerus. This should be investigated in future with other methods, e.g., AFLP or microsatellites.

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., 2013, 2015; Schyra 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, Stipagrostis garubensis, St. seelyae). The ‘East Gariep’ group range is shared by Aloe gariepensis, Callistigma inachabense, Commiphora gracilifrondosa, Maerua gilgii, and Zygophyllum decumbens. The combined ranges of the ‘Nama’ group and the ‘Southwestern Kalahari’ too contain localised plant taxa (Acrotome pallescens, Monechma calcaratum, Monsonia parvifolia, Petalidium parvifolium, and Wellstedia dinteri). The ‘Northern Namib’ group range is home to Commiphora wildii, Helichrysum roseo-niveum, Hermbstaedtia spathulifolia, Sarcocaulon mossamedense, Stipagrostis damarensis, and Stipagrostis uniplumis var. intermedia. The ‘Western Kalahari 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.