Old origin for an European‐African amphitropical disjunction pattern: New insights from a case study on wingless darkling beetles

The origin of the amphitropic Mediterranean Basin and southern African disjunction (European–African amphitropical disjunction; EAAD) pattern is generally attributed to recent dispersal events. However, our knowledge is limited because the origin of the EAAD pattern has been almost exclusively studied in plants. Here, we investigate the origin of this wide‐ranging disjunction pattern in a group of wingless insects, consisting of two major clades, both of which have EAAD distributions.


| INTRODUC TI ON
Wide-ranging geographic disjunctions between seemingly closely related groups of organisms were noted by Darwin (1859). Since then, researchers have documented many additional lineages exhibiting disjunct distributions across the tree of life. Understanding the mechanisms responsible for the origin of these striking distributional patterns has long been regarded as one of the main goals of biogeography (e.g. Cox & Moore, 2010;Lomolino et al., 2010). One of the most studied trends in biogeographic research of wide-ranging disjunction patterns concerns arid Mediterranean-type regions (e.g. Buerki et al., 2012;Coleman et al., 2003;Cowling et al., 1996;Simpson et al., 2017).
According to the available biogeographic classifications of the world, five Mediterranean-type regions are recognized: The Cape Floristic Region, Central Chile, Southwest Australia, the California Floristic Province and the Mediterranean Basin (Buerki et al., 2012;Coleman et al., 2003;Cowling et al., 1996). These areas are believed to have derived from temperate or tropical ancestral biomes (e.g. Raven, 1971) thanks to environmental changes associated with the consecutive drops of atmospheric carbon that started after the Middle Eocene climatic optimum, and after the Eocene-Oligocene transition (EOT) (Zachos et al., 2001). The Mediterranean-type regions are characterized by high physiographic diversity and hot, dry summers, which contrasts with cold, wet winters (Médail & Myers, 2004). Together these abiotic factors create a wide variety of environments suitable for many different groups of organisms. Although the Mediterranean-type regions occupy <5% of the Earth's land surface, they are inhabited by nearly 20% of all known species of vascular plants (Cowling et al., 1996). Moreover, across these five regions at least 40% of plant species are believed to be endemic (Raven, 1971), whilst in the Cape Floristic Region this number reaches up to 80% (Goldblatt, 1978). Due to their high levels of endemicity and the threat of human-driven habit change, all five Mediterranean climatic regions are considered biodiversity hotspots (Myers et al., 2000).
Despite their geographic distance, Mediterranean-type regions share many phylogenetically related taxa. This pattern has most frequently been reported for various groups of plants (e.g. Drew et al., 2017;Jürgens, 1997;McGuire & Kron, 2005;Wen & Ickert-Bond, 2009). Of particular interest is the origin of the amphitropic Mediterranean Basin and southern African disjunction (European-African amphitropical disjunction; EAAD). This pattern was first attributed by Balinsky (1962) to dispersal events facilitated by the expansion of a corridor of drier habitats during the Pleistocene glacial periods. However, this corridor may be much older and may have opened multiple times. For instance, dispersal between the Mediterranean Basin and southern Africa could be associated with the well-documented late Miocene/Pliocene spread of C 4 grasslands (Edwards et al., 2010). An increase in aridification is also associated with the East African uplift, which started some 8 million years ago (Ma) (Couvreur et al., 2008); the resulting East African Plateau also constituted a major migration route between southern Africa and the Mediterranean Basin (Galley et al., 2007). Dispersal between the Arabian and Eurasian plates was also facilitated by the Miocene closure of the Tethys (ca. 14-20 Ma; Hamon et al., 2013), leading to major faunal interchanges via the Arabian plate (Couvreur et al., 2021). Although the role of dispersal in generating EAAD patterns is often advocated, one should not discard the potential role of vicariance, especially in relation to the fragmentation of the pan-African rainforest (Morley, 2000) during drier periods of the Late Oligocene, mid-Miocene and Pliocene (Couvreur et al., 2021;Pokorny et al., 2015). Vicariance events associated with the repeated fragmentation of major forest blocks have often been suggested as the origin of major trans-African disjunct distributions, including the EAAD (see Couvreur et al., 2021 for a review).
Most biogeographic studies to date have been done in plants, and support a predominant role of dispersal from southern Africa to the Mediterranean Basin through an eastern African corridor (e.g. Calviño et al., 2006;Coleman et al., 2003;Désamoré et al., 2011;Goldblatt et al., 2002;del Hoyo et al., 2009;Kandziora et al., 2017;Sanmartín et al., 2010; but see Carlson et al., 2012;McGuire & Kron, 2005 for patterns of southward dispersals). Although recent dispersal events are sometimes inferred during the Pleistocene (e.g. Coleman et al., 2003;Désamoré et al., 2011) or the Pliocene (Kandziora et al., 2017) yet, these studies did not cover groups with extant EAAD patterns.
In insects, little is known as only a handful of studies have tackled the issue of the origin of EAAD taxa. A review on the origin of patterns of EAAD taxa in meloid beetles was conducted by Bologna et al. (2008), who suggested that the disjunct distribution exhibited by several meloid genera resulted from both northern and southern dispersal events between the Miocene and Pleistocene. Another study focusing on two weevils genera (Hernández-Vera et al., 2013) recovered a northern dispersal/Late Miocene scenario, attributed to distinct aridification events.
Adults of this species-rich tribe (ca. 400 species in 38 genera) range in size from 2 to 20 mm, are saprophagous, and inhabit various arid and semi-arid ecosystems in Africa and Western Europe ( Figure 1); however, the highest species diversity is observed in the Cape Floristic Region and the Mediterranean Basin (Kamiński, 2011;Koch, 1956;Medvedev, 1968).
Based on these recent studies this tribe is subdivided into two major clades (interpreted as subtribes-Dendarina and Melambiina), both of which have amphitropical African distributions ( Figure 1). What seems to be unique in comparison to other known EAAD groups (e.g. Coleman et al., 2003;Biondi & D'Alessandro, 2008;Bologna et al., 2008), is that the Afrotropical and Palaearctic complexes do not share similar genera, which could provide evidence of earlier dispersal and/or vicariance events. Six distinct generic complexes are known: they group morphologically similar genera (see Appendix S1 in Supporting Information) that also have similar distribution areas  Figure 1).
Since the historical biogeography of Dendarini has never been investigated, this group can be considered as a potentially rich and unexplored source of data to understand the origin of the EAAD.
In order to analyse this information, our study relies on dated phylogeny and historical biogeography to infer the ancestral range evolution of all major lineages of Dendarini (i.e. generic complexes; see Appendix S1). The following questions are addressed: (i) Did This dataset is primarily based on a recent study focused on pedinoid beetles, i.e., Dendarini, Opatrini, Pedinini and Platynotini (see ; it was modified to retain one specimen per species to better accommodate the historical biogeography analyses. In the matrix, Dendarini is represented by 27 species from 14 genera. This Dendarini sampling encompasses all major generic/biogeographic complexes traditionally designated among F I G U R E 1 Global distribution and species diversity of Dendarina (a) and Melambiina (b) in Sub-Saharan Africa and the Mediterranean region. After Iwan et al. (2011), Kamiński (2011), Kamiński, Kanda, et al. (2018), , Iwan et al. (2020), Kamiński, Lumen, et al. (2021). The range of sampled taxa is highlighted in blue. For classification details see Appendix S1. Selected representatives of the six generic complexes (c) are presented for illustrative purposes (pictures by M. Kamiński). Abbreviations: g.c., generic complex. Projection: Lambert cylindrical equal area projection both subtribal components ( Figure 1). Our sampling also includes a representative of each of the three Dendarini genera endemic to the Canary Islands (see above). Since no fossil specimens of Dendarini are known, this study relied on external calibrations. Members of the closely related tribe Opatrini are relatively well-represented in the fossil record (Nabozhenko, 2019), which enables molecular dating analyses based on primary fossil calibrations. Representatives of Pedinini and Platynotini were also included in the analyses to provide a comprehensive sampling from a taxonomic viewpoint. According to recent phylogenetic studies, Dendarini, Pedinini, Platynotini and Opatrini form a group of closely related tribes (coined as the 'opatrinoid' clade;  within the newly defined subfamily Blaptinae (Kamiński, Lumen, et al., 2021).

| Phylogenetic analyses
Phylogenetic analyses were conducted using both Bayesian inference (BI) and maximum likelihood (ML), using MrBayes v3.2.6 (Ronquist et al., 2012) and RAxML v8.2.8 (Stamatakis, 2014), respectively. For both analytical approaches, we carried out partitioned analyses to improve phylogenetic accuracy. The concatenated dataset was divided a priori into 14 partitions: we used three partitions (one per codon position) for each of the four coding genes and one partition for each of the ribosomal RNA genes. Best partitioning scheme and substitution models were determined using PartitionFinder v2.1.1 (Lanfear et al., 2017) using the default search algorithm ('greedy' option) and the most complete set of models ('model=all' option); we also used the 'linked branch lengths' option, where each subset has its own rate multiplier parameter but only one underlying set of branch lengths. For both partition and model selection, we relied on the Bayesian information criterion following Ripplinger, and Sullivan (2008).
For ML partitioned analyses, the best-scoring tree was obtained using a heuristic search implementing 100 random-addition replicates. Clade support was first assessed using standard nonparametric bootstrap support (BV) values, with 1000 replicates; nodes supported by BV ≥70% were considered as strongly supported following Hillis and Bull (1993). In addition, we implemented the transfer bootstrap expectation (TBE) method, which provides a better measure of branch repeatability, or robustness (Lemoine et al., 2018). For TBE values, we also used a 70% threshold following Lemoine et al. (2018). For BI partitioned analyses, two independent runs with eight MCMC (one cold and seven incrementally heated chains) were conducted: they ran for 50 million generations, with trees sampled every 5000 generations. A conservative 25% burn-in was applied after checking for stability on the log-likelihood curves and the split-frequencies of the runs. Support of nodes for BI analyses was provided by clade posterior probabilities (PP) as directly estimated from the post-burn-in majority-rule consensus topology; nodes supported by PP ≥ 0.95 were considered as well-supported following Erixon et al. (2003).

| Divergence time estimation
Divergence times were estimated under BI with BEAST v1.10.4 . For this study, we relied on two complementary approaches.
First, we implemented a node dating approach based on external fossil calibrations. Following the recommendations of Parham et al. (2012) for fossil calibrations, we conservatively retained two fossils (see Appendix S2). Priors for the two constrained nodes were conservatively enforced using uniform distributions with minimum ages provided by the fossil constraints and a maximum age conservatively set to be twice the age of the oldest known indisputable fossil of  et al. (2014) was also enforced on the root for the analyses to converge on a more restricted range of solutions (Sauquet, 2013). We used distinct molecular clocks for each of the mitochondrial and the nuclear genes, with clocks and substitution models selected with the 'beast' set of models of PartitionFinder (Lanfear et al., 2017).
Second, we carried out a secondary calibration approach based on substitution rates. Here, we relied on substitution rates inferred by the fossil-based study of Andújar et al. (2012), which provided detailed rates (ucld.mean estimates with 95% credibility intervals) for different mitochondrial and nuclear gene fragments. These rates have been used for secondary calibration in diverse beetle groups for which unequivocal fossil data is lacking (e.g. Caterino & Langton-Myers, 2018;Haran et al., 2021). To properly implement the corresponding rates of evolution we relied on a calibration strategy based on two clocks, one for the mitochondrial and one for the nuclear genes. Following Andújar et al. (2012), we tested combinations of uncorrelated lognormal (UCLN) and strict clocks (SC), resulting in four distinct calibration strategies. The fit of each corresponding calibration strategy was compared using Bayes factor (B F ); for this, we performed a marginal likelihood estimation using a path sampling procedure, with default settings (100 path steps, chains running for 1 million generation with a log-likelihood sampling every 1000 cycles).
For all analyses, the tree model was set to a birth-death speciation process (Gernhard, 2008) to better account for extinct and missing lineages. To limit the risk of over-parameterization, we also enforced a fixed topology corresponding to the topology with the highest support (this topology corresponds to the best-score tree obtained under ML). All corresponding BEAST analyses were per-
We used the BioGeoBears (Matzke, 2014) implementation, which can run and compare Dispersal-Vicariance Analysis (DIVA; Ronquist, 1997), Dispersal-Extinction-Cladogenesis (DEC; Ree & Smith, 2008) and Bayesian inference of historical biogeography for the discrete area (BayArea; Landis et al., 2013) approaches. For each approach (referred to as DEC, DIVAlike and BayArealike under BioGeoBears/ RASP), we also carried out additional analyses with models implementing founder-event speciation (+J parameter); the latter is particularly appropriate when dealing with oceanic islands or when wide areas are considered (Matzke, 2014 was set for the analyses. To account for major periods of paleoenvironmental changes, we used four distinct slices (see Appendix S3). Model selection followed Yu et al. (2020) and was conducted with likelihood ratio tests (LRT) for the nested models (null models vs. models with a +J parameter) whereas non-nested models (DEC vs. DIVAlike vs. BayArealike) were compared using the weighted F I G U R E 2 Phylogeny of Dendarini and the most closely related tribes inhabiting Afrotropical, Nearctic, and Palearctic regions. For details, see . Results of the ML analyses of the combined dataset. Support values are provided on nodes (TBE values in dark blue and BS values in black; PP are also provided in red for nodes shared with the post-burn-in majority-rule consensus topology resulting from the BI analyses), asterisks are used to indicate maximum support values. On the top left a picture of Litoborus moreleti (Lucas) is displayed for illustrative purposes (picture by M. Kamiński) Akaike information criterion corrected for sample size (AICc_wt).
Four additional analyses were also carried out with models (see also Appendix S3) modified to favour temporally distinct dispersal corridors corresponding to the alternative hypotheses presented in the Introduction: (i) A recent dispersal corridor during the Pleistocene,

| Phylogenetic analyses
The best-fit partitioning scheme recovered by PartitionFinder consists of eight partitions (see Appendix S4). ML and BI analyses yielded almost identical topologies, however, support values for the ML tree were higher than those of the BI tree (respectively, 69 and 60 of the 71 nodes were supported by TBE and BS ≥70% whereas 58 of the 71 nodes were supported by PP ≥ 0.95). The ML best-score tree (L = −46973.15) is represented in Figure 2 (see also Figure S1.5 for the post-burn-in majority-rule consensus topology resulting from BI analyses). The only difference between ML and BI analyses relates to the placement of Gridelliopus subsquamosus, which is recovered in a more derived position in the ML topology. All major tenebrionid lineages (at the tribe level) are recovered as monophyletic, generally with high support. Within the Dendarini, the two subtribes and all generic complexes are recovered monophyletic with high support. With the exception of Pedinus and Zadenos, the monophyly of all genera for which more than one species was included is also recovered. The three genera that are endemic to the Canary Islands are also recovered monophyletic with maximum support (TBE and BS of 100% and PP of 1.0).

F I G U R E 3 Dated phylogeny and historical biogeography of
Afrotropical and Palearctic components of Dendarini. Ancestral area reconstructions correspond to the result of a time-stratified RASP analysis with the best-fit model (DIVAlike+J model). Circles on nodes display the relative probabilities of ancestral ranges (ranges with probabilities below 5% are hidden and lumped together; they are reported in black). Most likely states are figured inside the circles using the coding of areas presented beforehand. Inferred dispersal and vicariance events are highlighted on nodes using the initials 'd' and 'v'. The four time slices used in this analysis are displayed using dotted lines. Black arrows underline the diversification of lineages representatives of the six Dendarini generic complexes; coloured rectangle squares associated with each arrow highlight the distribution range of each generic complex TA B L E 1 Pairwise comparison of the results of RASP historical biogeography analyses. The best-fit model is highlighted using a bold font. Nested models (i.e. alternative models vs. null models) were compared using a likelihood ratio test (LRT), whereas non-nested models (DEC vs. DIVAlike vs. BayArealike) were compared using the weighted Akaike information criterion corrected for sample size (AICc_wt)

| Ancestral areas reconstruction
Likelihood ratio tests indicate that +J models provide the best fit across all analyses (see Appendix S8). AICc_wt model selection also further supports the DIVAlike+J model over the DEC+J and BayArelikea+J models (Table 1). The ancestral areas reconstruction associated with the DIVAlike+J model (see Figure 3 and Figure

| Evolutionary relationships of Dendarini
The conducted analyses support the monophyly of Dendarina and Melambiina and their major geographic components/generic complexes, confirming that the observed amphitropical African distribution of both subtribes is not an artefact caused by flawed taxonomic concepts. This is not surprising as the distinctiveness of all major generic complexes employed in this study is clearly grounded in morphology (e.g. Iwan et al., 2011;Kamiński, 2017;Kamiński, Lumen, et al., 2021;Kamiński, Kanda, et al., 2018;Koch, 1956). However, whilst the monophyly of Dendarina and its main components could have been predicted based on the available morphological data alone, the phylogenetic relationships between specific lineages within the tribe remained unexplored until recently. Unsurprisingly, the results obtained here are almost identical to the ones presented by . However, the small differences in terms of data treatment and the reliance on a dataset including only one specimen per species (vs. a dataset including several specimens for some species in ) yielded a single-yet potentially important-difference.
Specifically, our ML analysis recovered Gridelliopus subsquamosus (representative of the Silvestriellum g.c.) sister to the Melambius g.c. (see Figure 2), whilst the BI analysis ( Figure S1.5), as well as  results, place it sister to all remaining Melambiina. The uncertainty in this placement (supported by a TBE value of 95% under ML and a PP of 0.95 under BI) is likely partially driven by the fact that we were only able to recover two gene fragments (12S and COII) from the 59-year-old processed specimen of G.
subsquamosus. Nevertheless, in light of available morphological data, the recovered ML topology (Figure 2) appears as the most reliable one; indeed, both G. subsquamosus and representatives of Melambius g.c. exhibit an emargination of the fifth abdominal ventrite, which is absent on all members of the Zadenos g.c. (Koch, 1956).

| An ancient origin for the EAAD pattern
The results of this study illustrate the potential of darkling bee- young lineage with respect to other darkling beetle tribes (Kergoat et al., 2014), they are comparatively much older than many of the insect and plant groups that have previously been used to tackle the question of the origin of EEAD patterns (e.g. Carlson et al., 2012;Coleman et al., 2003;Désamoré et al., 2011;Goldblatt et al., 2002;del Hoyo et al., 2009;Kandziora et al., 2017;Pokorny et al., 2015;Sanmartín et al., 2010). It is only thanks to the age of the lineage that we were able to test and unravel an interesting scenario related to the origin of EAAD.
The main novelty recovered here concerns the temporal context of the inferred biogeographical events (Figure 3).   in Singinda, Tanzania (e.g. Jacobs & Herendeen, 2004). The climate at that time in this area was estimated to be dry and comparable to those of modern Miombo woodlands (Jacobs & Heredeen, 2004), where a few species of Dendarini can be found nowadays (Kamiński, 2011). When the two main Dendarini lineages originated, it is possible that the EOT cooling event played a major role in their biogeographic history, accounting for by the disjunction resulting in their extant EADD patterns.
The EOT was one of the most extreme climatic events during the Cenozoic Era (Zachos et al., 2001), known for having generated or accelerated major shifts of floral communities worldwide (Pound & Salzmann, 2017). However, its precise influence on African climate and floral communities is still unclear. As underlined by Couvreur et al. (2021), there are uncertainties in atmospheric carbon dioxide concentration reconstructions during the Eocene and Oligocene (Steinthorsdottir et al., 2016). A recent study of palynological records also reported contradictory trends for East Africa (Pound & Salzmann, 2017), with either gradual changes predating the EOT or no changes during the EOT. However, as underlined by the review of Couvreur et al. (2021), both fossil data (e.g. Morley, 2000) and the results of diversification analyses on plant phylogenies (e.g. Bouchenak-Khelladi et al., 2010Faye et al., 2016) consistently indicate a major turn-over of floral communities during the EOT. We hypothesize that these environ- following the EOT (Couvreur et al., 2021). To our knowledge, the Dendarini model-system constitutes the oldest known example of EAAD origin as existing studies mostly point towards Pleistocene (e.g. Coleman et al., 2003;Désamoré et al., 2011), Pliocene (e.g. Kandziora et al., 2017) or Late Miocene events (e.g. Hernandez-Vera et al., 2013;Miguez et al., 2017). This probably explains whycontrary to other known EAAD taxa that usually belong to similar genera (e.g. Bologna et al., 2008)-the different geographical components of Dendarini are morphologically well-differentiated ( Figure 1); the latter is also clearly reflected in the classification of the group, as no common genera have been recognized for the Palaearctic and southern African components, while the Pythiopus g.c. was until recently treated as a separate-phylogenetically enigmatic-subtribe (Kamiński, Kanda, et al., 2018). Furthermore, the reliability of the hereby presented biogeographic scenario is also supported by the convergence of the results originating from different analytical methods. Finally, as the distribution of other Blaptinae representatives is relatively well-studied, and several species representing closely related tribes were reported from Equatorial Africa, this eliminates the possibility that the distributional pattern observed for Dendarini is artificial and caused by collection bias (Kamiński, 2015(Kamiński, , 2016Kamiński & Iwan, 2017).

| Distributional patterns
Comparison of the distributional patterns of Palaearctic components of Dendarina and Melambiina reveals that the first subtribe is more widespread in Europe (Figure 1). Even though the northernmost localities reported for Dendarina (Baltic and the North Sea coasts) concern a single extremely widespread species (i.e. Phylan gibbus Fabricius), it can be concluded that the greatest species diversity of that subtribe is located on the northern coast of the Mediterranean Sea (including the variable islands). Furthermore, unlike Melambiina, Dendarina is also widely represented in Transcaucasia (Figure 1). On the contrary, the greatest species diversity of Palaearctic Melambiina is located in Morocco. Contrary to Dendarina, the group is also present on the Canary Islands with various endemic species. Only a few Melambiina species are recorded from Europe (Kamiński, 2011).
When analysing the distribution of the sub-Saharan components of the two subtribes, the absence of Melambiina on Madagascar and the lack of Central African Dendarina species is notable (Figure 1).
Due to the absence of ecological contributions on the tribe, the aforementioned differences in distributions cannot be easily attributed to any specific variations in niche preferences (Kamiński, 2011). The representatives of both subtribes are known to co-occur in some European and southern African habitats, whilst the general appearance of the areas where Melambiina and Dendarina were collected seems largely convergent (Appendix S10). Furthermore, the larval morphology of only seven species of Dendarini is currently known, precluding any insights from this perspective . Instead, the currently gathered ecological and biological data point to similarities between Dendarina and Melambiina. Namely, taking into consideration the xerophilic nature of Dendarini, it can be safely assumed that all representatives of this tribe have soil-inhabiting larvae (Matthews et al., 2010). Adults and larvae of xerophilic species within darkling beetles are usually generalist detritivores; therefore, they are less impacted by competition for resources and can disperse and thrive in a wide variety of environmental conditions (e.g. Fattorini et al., 2020). Xerophilic tenebrionid species can thus be considered as more opportunistic colonizers, contrary to insect lineages with distinct feeding habits which are usually constrained by the availability of a specific feeding source, as is the case for most phytophagous insect lineages (e.g. Adler & Dudley, 1994), where patterns of delayed colonization are sometimes evidenced due to the contrasting biogeographic histories of insects and their plant hosts (e.g. Kergoat et al., 2015).
As underlined beforehand, the majority of known Dendarini are flightless, which constitutes a potential driver of allopatric speciation (e.g. Ikeda et al., 2012). Interestingly, and also perhaps counterintuitively, loss of flight did not preclude oceanic dispersal in multiple darkling beetle lineages (e.g. see Juan et al., 1995;Condamine et al., 2013 (Kergoat et al., 2014). Moreover, there are no disputes over the ownership of the data presented in the paper and all contributions have been attributed, via co-authorship or acknowledgement, as appropriate to the situation.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available in permanent open repositories (GenBank for the molecular data and