Fine‐scale bee species distribution models: Hotspots of richness and endemism in South Africa with species‐area comparisons

While global patterns of bee diversity have been modelled, our understanding of fine‐scale regional patterns is more limited, particularly for under‐sampled regions such as Africa. South Africa is among the exceptions on the African continent; its bee fauna (ca. 1253 species) has been well collected and documented, including mass digitising of its natural history collections. It is a region with high floral diversity, high habitat heterogeneity and variable rainfall seasonality. Here, we combine a South African bee species distributional database (877 bee species) with a geospatial modelling approach to determine fine‐scale (~11 × 11 km grid cell resolution) hotspots of bee species richness, endemism and range‐restricted species. Our analyses, based on the probabilities of occurrence surfaces for each species across 108,803 two‐minute grid cells, reveal three bee hotspots of richness: Winter rainfall, Aseasonal rainfall and Early‐to‐late summer rainfall. These hotspots contain large numbers of endemic and geographically restricted taxa. Hotspots with particularly high bee diversity include the Fynbos, Succulent Karoo and Desert biomes; the latter showing 6–20 times more species per unit area than other biomes. Our results conform with global species‐area patterns: areas of higher‐than‐expected bee density are largely concentrated in Mediterranean and arid habitats. This study further enhances our knowledge in identifying regional and global hotspots of richness and endemism for a keystone group of insects and enabling these to be accounted for when setting conservation priorities.


INTRODUCTION
Bees are keystone pollinators for both wild flowering plants (Neff & Simpson, 1993;Kearns & Inouye, 1997;Ollerton et al., 2011) and agricultural crops (Garibaldi et al., 2013;Klein et al., 2007).Despite their ecological and economic importance, a major drawback in developing effective conservation priorities is our limited understanding of the patterns of bee diversity at finer scales, including local hotspots of richness and range-restricted taxa.To remedy this, there has been an ongoing effort to mine data from natural history collections (Bartomeus et al., 2018;Bystriakova et al., 2018;Orr et al., 2021).Such data have also been used as an effective tool in assessing longterm changes in bee diversity and species declines due to global change, and in identifying priority areas for their conservation (Bartomeus et al., 2013(Bartomeus et al., , 2018;;Critchlow et al., 2022;Kougioumoutzis et al., 2022;Soroye et al., 2020;Zattara & Aizen, 2021).
Before the mass digitization of natural history collections, bee biologist Charles Michener, using his accumulated expertise and interpretation of the literature, first described global patterns of bee species richness as highest in warm temperate, xeric regions of the world (Michener, 1979).This pattern contrasts with the diversity of many other taxonomic groups, which are concentrated in the tropics following a latitudinal gradient of diversity (Gaston, 2000).A recent global analysis modelled bee species richness using digitised natural history collection data (Orr et al., 2021).It confirmed Michener's descriptive patterns, finding bee species richness to be highest at mid-latitudes, especially in warm temperate and xeric regions, showing a bimodal latitudinal richness gradient (Orr et al., 2021).
In addition to this study, several macroecological analyses of bee richness have been done across global or biogeographic regions (Bystriakova et al., 2018;Kaloveloni et al., 2018;Leclercq et al., 2023;Michez et al., 2009;Patiny & Michez, 2007;Williams, 1998).Surprisingly, few studies investigating fine-scale (approximately <30 Â 30 km grid cell) patterns of bee richness and endemism have been undertaken at the country level (i.e., regional scale) (but see Kuhlmann (2009); Kougioumoutzis et al. (2022); Patiny and Michez (2007)).Analysis at this geographic scale requires data of sufficient quantity and precision using the widest taxonomic coverage.Although regions of interest could potentially be extracted from currently available global species richness maps (e.g., Bystriakova et al., 2018;Orr et al., 2021), there are several limitations to using these datasets to assess local patterns within regions.These include using coarse-resolution data (e.g., ranging from 110 Â 110 to 550 Â 550 km grid cells), using a taxonomic subset of the data (e.g., family Colletidae) and projecting modelled Western Hemisphere bee data on to Eastern Hemisphere regions.Coarse-resolution data are likely to obscure fine-scale richness patterns, particularly in regions with high habitat heterogeneity (Cowling et al., 2017;Cowling & Lombard, 2002).The subfamily Colletinae has been put forward as a taxon to represent richness patterns for all bees (Bystriakova et al., 2018); however, this could introduce spatial biases.For example, the richness pattern retrieved for South Africa by Orr et al. (2021), which used a more comprehensive dataset of bees, strongly contrasts with Bystriakova et al.'s (2018) richness pattern.Projecting Western Hemisphere bee models onto Eastern Hemisphere regions (e.g., Orr et al., 2021) may be problematic due to inherent biotic differences, distinct species compositions and different evolutionary histories (De Palma et al., 2016;Kreft & Jetz, 2007;Michener, 1979;Ricklefs, 2004).In addition, the diverse environmental variability between hemispheres, encompassing climate disparities and variations in habitat structures (Kreft & Jetz, 2007;Wiens, 2011), poses a challenge when extrapolating models based on geographically and taxonomically restricted data to other regions and taxa (De Palma et al., 2016;Elith & Leathwick, 2009).These limitations become less significant in data-rich regions, such as South Africa (Melin & Colville, 2019), where it would be more parsimonious to use the source data to model potential species distributions and quantify patterns of diversity.
The main aim of this paper is to explore patterns of bee species richness and endemism within South Africa, at a fine-scale resolution ($11 Â 11 km grid cell), and using a uniquely comprehensive digitised bee dataset that included representatives from all six bee families that occur in the region.South Africa is a region of high species richness (ca.1253 species (Eardley & Coetzer, 2023)) and the Mediterranean and arid habitats in the southwest of the country are particularly rich in bees (Bystriakova et al., 2018;Kuhlmann, 2009;Michener, 1979Michener, , 2007;;Orr et al., 2021).The South African bee fauna is characterised by high levels of endemism (Bystriakova et al., 2018;Kuhlmann, 2009) and is regarded as evolutionary significant since it includes some of the most early diverging bee lineages (Danforth et al., 2006;Litman et al., 2011), as well as genera showing specialised pollinator adaptations (Steiner & Whitehead, 1990;Whitehead & Steiner, 1993).
Uniquely within the African continent, the bee fauna of the region has received a substantial amount of focused collection and taxonomic effort over the past 250 years, including mass digitization of its natural history collections (Melin & Colville, 2019).In this study, we use publicly available natural history collection data and employ a geospatial modelling approach (Colville et al., 2020;Mecenero et al., 2015).This allows us to produce the most detailed assessment of South African bee diversity and distribution to date, with the aim of delimiting local hotspots of bee species richness and endemism.We also used richness-weighted indices to identify patterns of richness for restricted and widespread species (Crisp et al., 2001).Furthermore, we determined the extent to which there is spatial congruence among patterns of species richness, national endemic richness and range-restricted species richness, as well as between national endemic richness and range-restricted endemic richness.
The high bee diversity of Mediterranean and arid climate regions has attracted much attention over the past few years (Kaloveloni et al., 2018;Leclercq et al., 2023;Marshall et al., 2023;Meiners et al., 2019;Orr et al., 2021).We place South Africa's bee diversity within this global pattern and compare how richness scales with area (Gotelli & Colwell, 2001) within South Africa-taking into account different biomes and rainfall seasons-and with studies from other arid and Mediterranean regions worldwide.Our approach involves estimating species density of bees using a species-area relationship, building on recent studies that have employed this methodology (Meiners et al., 2019;Minckley & Radke, 2021).Comparing bee density across these different geographic areas can provide insights into understanding global patterns of bee diversity (Minckley & Radke, 2021).
Detailed, modern biogeographic studies on South Africa's insect diversity are lacking (Colville et al., 2014;Melin & Colville, 2019).Consequently, insects are underrepresented in conservation planning, primarily due to the large numbers of species involved, and the outdated taxonomy for many groups, compounded by a decline in taxonomic expertise (Chown & McGeoch, 1995;Hamer, 2012;Scholtz & Chown, 1995).The results of this study should further enhance our knowledge in identifying important areas of insect richness and endemism in South Africa at a fine geographical scale and enabling these hotspots to be accounted for when setting conservation priorities.
Given the idiosyncratic patterns of bee richness (e.g., highest in Mediterranean, xeric and temperate habitats) (Bystriakova et al., 2018;Michener, 2007;Orr et al., 2021), finer-scale data may also help guide future studies, which aim to understand the macroecological factors driving bee diversity.

Determining hotspots of bee richness from species distribution models
We obtained species locality records from digitised natural history collections, mostly (86%) from South African institutions: Albany Museum, Iziko Museums of South Africa, the National Collections of Insects and Ditsong Museums of South Africa.The remainder of the data is derived from institutions outside of South Africa; publicly available on GBIF.org (GBIF.org, 2023a(GBIF.org, , 2023b(GBIF.org, , 2023c(GBIF.org, , 2023d, 2023e, 2023f), 2023e, 2023f).The South African datasets constitute over 250 years of collecting bees including ad hoc and targeted collecting by museums and entomologists.Melin and Colville (2019) and Coetzer and Eardley (2019) provide detailed historic information including a summary of institutional and collector contributions.Together this original dataset is comprised of 49,498-point locality records, which underwent a range of quality checks to improve its reliability and accuracy.For common geospatial issues found in natural history collections (Robertson et al., 2016) such as point-localities outside of their range (including points out at sea) for South Africa, we verified their accuracy and either corrected or removed them if georeferencing was not feasible.Similarly, disjunct distributions were identified as any data points >200 km from its nearest neighbour.These were flagged as potential outliers and manually checked and corrected.
We restricted the data to include only valid species by removing records without a species name (e.g., spp.and sp.).To maintain taxonomic identity, we followed the published regional checklist for South African bees (Eardley & Coetzer, 2023) and concerning synonymy and spelling we followed Eardley & Urban (2010) and Coetzer and Eardley (2019).Our cleaned dataset consisted of 26,474-point locality records across 877 bee species (ca.70% of South African species), of which 234 species had less than five-point locality records.
For South African bee species that also occur in Lesotho and Eswatini, we retained the occurrence records in our database for spatial continuity.
Several biases are inherent in the presence-only or 'atlas-type' data used here (Bartomeus et al., 2018;Graham et al., 2004).For example, collector effort might be biased towards areas that are accessed more easily (e.g., towns, cities and road verges) as opposed to areas that are difficult to access.Also, collectors may target species of special interest over common taxa (Bartomeus et al., 2018;Graham et al., 2004;Mecenero et al., 2015;Melin & Colville, 2019).Another limitation of natural history collection data is that it provides information on where species are present, not where they are absent.However, biases can be reduced by combining such data with appropriate statistical analyses (Araújo et al., 2019;Barbet-Massin et al., 2012;Colville et al., 2020;Pearce & Boyce, 2006).To account for these biases, we employed a geospatial modelling technique (Mecenero et al. (2015) and recently used by Colville et al. (2020)) to interpolate the distribution records for each bee species and to calculate a continuous probability of occurrence surface for each species at a 2-min grid cell scale ($11 Â 11 km).Grid cells of this scale are considered an appropriate spatial scale for capturing the highly heterogeneous habitats of South African landscapes (e.g., Colville et al. (2020)), such as steep upland/lowland gradients and floristic and edaphic transitions.
Using interpolated species distributions offered advantages over using raw presence-only data, as our measure of richness and endemism were not overly biased by gaps in the data (i.e., false absences).
We followed the same modelling approach ('Spatial Model 1') and code described in detail by Mecenero et al. (2015) and that had built on earlier models by Hengl et al. (2009).In summary, for each bee species, we built a model at 2-min resolution combining point pattern analysis methods with environmental niche information, to account for ecological similarity, inferred observer effort and geographical distance.This process involved two stages.The first stage involved selecting a sample of non-focal species records to act as pseudo-absences (reflecting the pattern of observation in the dataset), and the second stage involved interpolating distributions based on presence and pseudo-absence records.The first stage consisted of several separate steps: (1) mapping all records of the focal species and generating a kernel density estimate for records of this species, (2) identifying all records of all other bee species >100 m from records of the focal species and generating similar kernel density estimates, (3) computation of the difference in density estimates between focal and non-focal species (an approximate index of the probability of encountering the focal species), (4) computation of an environmental envelope within a principal component analysis of rainfall (mean annual rainfall and rainfall season) (Schulze, 2007) and temperature variables (mean winter and mean summer temperature) (Schulze, 2007), (5) computing the environmental distance between all 2-min raster cells and the centroid of the environmental envelope occupied by the focal species and (6) sampling records of the non-focal species using the environmental distance and geographic probability of encountering the focal species to bias selection towards locations where absence was most likely.Once pseudo-absence records were selected, the second stage of analysis used regression kriging of the presence/absence points onto the 2-min raster surface, using the rainfall and temperature covariates.For 234 bee species recorded from <5 locations, we were unable to accurately interpolate distribution and we therefore generated a raster map with presence (1) and assumed absence (0) directly from the raw data.One of us (A.M.) verified the modelled species distributions for groups they have detailed knowledge on by checking the maps against known distributions.Mecenero et al. (2015) undertook the same approach and found that modelled species distributions matched expert opinion, stressing the general accuracy of their modelling approach.
Our final modelled bee dataset consisted of probabilities of occurrence surfaces for each bee species across 108,083 2-min grid cells.Species richness and endemism were calculated for each 2-min grid cell as the summed probability surfaces for all our bee species.

Determining hotspots of endemism and rangerestricted species
To determine endemic species richness, we compiled a national list of endemic bee species (646 species) from Eardley and Coetzer (2023).All national endemic species are not known to occur in Lesotho or Eswatini.However, our grid cells fell across political boundaries, and we therefore allowed our models to predict species distributions across these artificial boundaries.Of the 646 listed endemic bee species, our modelled species distributional dataset contained 341 ($53%).As with overall species richness, we summed the probability of occurrence for each endemic species in a grid cell to obtain a measure of endemic bee species richness.
We determined hotspots of range-restricted species richness for all species and for national endemics using two indices (weighted endemism [WE] and corrected WE [CWE]) proposed by Crisp et al. (2001) and modified to accept our probability of occurrence for each species.CWE accounts for the correlative role of species richness in determining endemic and range-restricted richness patterns.To estimate WE, we weighted the probability of occurrence for each species by the inverse of their range (1/number of cells occupied by a species, summed across all species in a cell), and for CWE we correct WE for overall species richness in a cell (WE/cell richness).

Congruence of species richness with endemic richness and range-restricted species richness
We used Pearson's correlations to compare the spatial similarity of overall species richness with national endemic richness and rangerestricted species richness.Similarly, national endemic richness was correlated with range-restricted endemic richness.

Bee density: Biomes and rainfall seasonality
We calculated bee richness for South Africa's biomes (Mucina & Rutherford, 2006) and for rainfall seasonality (Schulze, 2007) by performing spatial overlays.Biomes and rainfall seasonality have been considered important drivers of biogeographic patterns of plant and insect diversity (Colville et al., 2014;Procheş & Cowling, 2007).
We calculated the species density for South African biomes and rainfall seasons by dividing the number of species by the respective area of the biome (n = 7; Figure 3a) or rainfall season (n = 6; Figure 3b).To compare density among biomes or rainfall seasons, we calculated the ratio by dividing density of one biome/rainfall season by another.

Bee density: Global comparisons
To compare South Africa's bee biome density with areas sharing similar habitat and climate conditions, we collated data for arid and Mediterranean sub-regions from nine published studies and one unpublished dataset (see Data S1).To also make comparisons at a global scale, we took advantage of recently collated (Orr et al., 2021) regional-level data for 170 countries with areas over 5000 km 2 (Data S1).To ascertain how richness scales with area, we used a species-area relationship as a double logarithmic equation (Rosenzweig, 1995, and as used recently in bee density studies by Meiners et al., 2019;Minckley & Radke, 2021) based on the species richness and area data presented in Data S1.
To assess whether studies conducted in arid and Mediterranean areas (n = 17, sub-region and biome data in Data S1) recorded more or fewer species than expected, we followed Meiners et al.'s (2019) and Minckley and Radke's (2021) approaches and measured the distance from each study's species-area point to the overall log-log regression line.With these measurements, we generated a bar plot to rank and compare how much the observed number of bee species deviates, either above or below, relative to the predicted number of bee species.
All analyses were carried out using the program R version 4.1.2(R Core Team, 2021) and packages: 'ggplot2' (Wickham, 2009) was used to plot the species-area relationship for different countries, biomes and sub-regions, and 'sp' (Bivand et al., 2013) 'spatstat' (Baddeley et al., 2015), 'rgdal' (Bivand et al., 2018) and 'gstat' (Graler et al., 2016) were used to perform species distribution models and spatial overlays.Species richness maps, as in Figures 1-3, were rendered in QGIS version 2.18.14 (QGIS Development Team, 2023).National endemic richness (unweighted) is concentrated in the western parts of South Africa, around the Winter rainfall bee hotspot and the Aseasonal rainfall bee hotspot (Figure 2b).Hotspots of unweighted range-restricted species (Figure 2c) and range-restricted endemic species (Figure 2d) are also concentrated around these two richness hotspots, although for the latter category, the Winter rainfall bee hotspot is extended as a hotspot of range-restricted endemic species.

Species
After weighting for richness of national endemics, both the Winter rainfall bee hotspot and the Aseasonal rainfall bee hotspot become more conspicuous, whereas the Early-to-late summer rainfall bee hotspot is less so, with only the eastern KwaZulu-Natal Drakensberg Mountains being rich in both species and national endemics.When correcting for the richness of range-restricted species, the Winter rainfall bee hotspot is emphasised as an area with a high proportion of bee species with narrow geographic ranges (Figure 2e).However, when correcting for range-restricted endemic richness a somewhat different pattern is retrieved (Figure 2f).The winter rainfall zone is retained but extends further into the central arid interior and includes elements of the Aseasonal rainfall bee hotspot as well as areas around the Early-to-late summer rainfall bee hotspot.
There was a significant positive correlation between species richness and national endemic richness (Pearson's R = 0.64, p < 0.001), as well as between species richness and range-restricted species richness (Pearson's R = 0.87, p < 0.001).Despite the high correlations among these three indices, there are some differences in the locations of the hotspots (Figure 2a-c).We also found a strong positive correlation between the spatial patterns of national endemic richness and rangerestricted endemic richness (Pearson's R = 0.93, p < 0.001; Figure 2b,d).
Comparisons across biomes and rainfall seasons showed that the number of bee species was similar (Table 1); however, when values were weighted on area size, species density was higher in the smallersized biomes and rainfall seasons (Figure 4a and Tables 2 and 3).The Fynbos, Succulent Karoo and Desert biomes, making-up $16% of the total area of South Africa (Mucina & Rutherford, 2006), had almost six times as many species per unit area than the three largest biomes (Savanna, Grassland and Nama Karoo), which make-up the bulk of South Africa (Table 2).The Desert biome showed the highest species density, with approximately six times more species per unit area than the Fynbos and Succulent Karoo biomes, and >20 times more species per unit area than the Grassland, Savanna and Nama-Karoo biomes.A similar pattern was seen for rainfall seasonality, with the winter and aseasonal rainfall areas ($15% of the total area of South Africa) having approximately two and five times more species per unit area than the summer rainfall areas, respectively (Table 3).
Area was a poor predictor of bee species richness explaining 6.6% of the variation in bee richness (F 1,184 = 14.13, p < 0.001) (Figure 4a).All South Africa's biomes fell above the regression line and outside of the 95% confidence interval indicating that they have higher bee species density than predicted by area (Figure 4a,b).Although all biomes fell above the regression line, the smallest biomes had the largest differences in the number of species observed relative to the number of species predicted (Figure 4b).These small-sized biomes are those with arid and Mediterranean climates.

DISCUSSION
Using the widest taxonomic coverage of South African bees (877 species across 85 genera) to date, we identified three main bee richness hotspots-two are closely aligned with the winter and aseasonal rainfall regions and the third with the summer rainfall region (Figures 1 and 3b).Hotspots of endemics were also predominantly located within the winter and aseasonal rainfall areas of the Greater Cape Floristic Region (GCFR), except for a small high-elevation hotspot of endemism found in the summer rainfall Drakensberg Mountains.Patterns of bee diversity have been assessed within South Africa to some extent, where Kuhlmann (2009) focused primarily on the bee families Apidae and Colletidae, Bystriakova et al. (2018) emphasised Colletidae, and Orr et al. (2021) used projected data for all families.This is the first fine-scale study using geospatial modelling to assess patterns of richness, endemism and range-restricted taxa.In addition, it is the first regional study to do this for any South African insect group.
We attribute South Africa's remarkably rich bee fauna in the Winter rainfall bee hotspot to the Succulent Karoo and Fynbos biomes (Figures 2a-f and 3a), and the Desert biome (Figures 2d-f and 3a).
Bee density in these winter rainfall biomes is considerably higher compared with the larger summer rainfall South African biomes (Figure 4a,b).This western region is characterised by having a high proportion of distinct bee genera: Scrapter, Fidelia and Patellapis (Danforth et al., 2013;Hedtke et al., 2013;Litman et al., 2011;Michener, 1979), which elevates this hotspot's profile as an important area for bee conservation.Added to this, the region is also of evolutionary significance for bees as a centre of diversity for Melittidae, a phylogenetically deep lineage of bees and sister to all other bees T A B L E 1 Number of bee species and national endemic species found in South African biomes and rainfall seasons.F I G U R E 4 Relationship between area and the number of bee species (log-transformed): (a) Each data point (n = 186; Data S1) represents a country, biome, and sub-region.The solid black line represents the linear regression fit, and the grey shaded region around the line shows the 95% confidence interval.Data points that fall above the line show higher than expected bee richness for their area size.Labels for South African biomes (circle with a star) and selected sub-regions have also been included for comparison and aesthetics; two countries have been highlighted: South Africa (this study; circle with a cross) and Israel (region with the highest bee density; circle with a dot).(b) Barplot of the difference in the number of bee species observed relative to the number of bee species predicted by the regression line plotted in panel A. Areas that recorded more bee species than expected are shown in blue with the exception of one area shown in red that recorded less bee species than expected.South African biomes are outlined in black.The following have been abbreviated: GSE National Monument, Grand Staircase Escalante National Monument; SA, South Africa; for states in North America: AZ, Arizona; CA, California; UT, Utah.Full study details are given in Data S1.(Danforth et al., 2013;Michez et al., 2008).The life history traits, such as the ancestral state of not lining their nests as seen in Fidelia, may explain the restriction of some lineages to arid environments (Litman et al., 2011).The winter rainfall Succulent Karoo and Desert (desert climates) and Fynbos (Mediterranean climate) biomes are recognised globally for their exceptional plant diversity, which is disproportionally high in comparison to its physical size and to other Mediterranean and desert regions (Cowling et al., 1998(Cowling et al., , 1999;;Manning & Goldblatt, 2012;Mittermeier et al., 1998;Myers et al., 2000).Ecological factors associated with the high plant diversity within these biomes include fine-scale climatic, edaphic and topographic heterogeneity, fire regime, and biotic interactions (e.g., plant-pollinator coevolution) (Cowling et al., 2017).Areas within this region have also been strongly influenced by historical climatic stability, which is thought to have lowered rates of extinction and allowed for higher niche specialisation in plants (Colville et al., 2020;Cowling et al., 1998Cowling et al., , 2015Cowling et al., , 2017)).This area is also known as a global centre for species richness for several other insect groups, for example, monkey beetles (Hopliini (Colville et al., 2018)), grasshoppers (Lentulidae (Otte, 2020)), weevils (Curculionidae (Hévin et al., 2022)), heelwalkers (Mantophasmatodea (Predel et al., 2012)) and several others (Colville et al., 2014).Findings to date suggest that plant species richness (Procheş et al., 2009;Wright & Samways, 1998), plant phylogenetic diversity (Procheş et al., 2009), soil-type preferences for nesting (Gess, 1992;Gess & Gess, 2014) and temperature (Botes et al., 2007;Stuckenberg, 1969) are important ecological factors in explaining insect diversity in this western region.The generally strong relationship between plant and insect species diversity is likely due to both direct interactions and where both groups are responding to similar environmental factors (Colville et al., 2014;Procheş & Cowling, 2006).
Similar ecological factors to the Winter rainfall bee hotspot are likely to influence the diversity of bees in the western part of the Aseasonal rainfall bee hotspot because it is composed of Fynbos and Succulent Karoo vegetation (Cowling et al., 2017).This southern hotspot is characterised by having a high proportion (>50%) of species of Compsomelissa and Allodapula; both have distinctive life history traits related to their sociality and progressive feeding (Michener, 1979;Schwarz et al., 2003).The eastern part of the Aseasonal rainfall bee hotspot is also likely to be influenced by a high degree of habitat heterogeneity due to the transitional nature of several floristic elements found here: Succulent Karoo, Albany Thicket, Nama-Karoo, Savanna and Grassland biomes (Figure 3a).This region, also known as the Albany Centre, has a very old phytogeographic origin (Clark et al., 2009;Cowling & Procheş, 2005;Van Wyk & Smith, 2001), indicating an area of historical climatic stability (Colville et al., 2020;Potts et al., 2013), which may also contribute to its higher than expected bee diversity.Vertebrates (Perera et al., 2011) and several insect groups, for example, butterflies (Carcasson, 1964), damselflies (Dijkstra, 2007) and grasshoppers (Dirsch, 1965;Gordon et al., 2023) also have high species diversity in this region.After adjusting for species richness, this hotspot maintains a higher concentration of geographically restricted taxa than would be expected based solely on species richness (Figure 2c,d).Although we have factored in collector effort through geospatial modelling, the Aseasonal rainfall bee hotspot might be the result of intensive collecting carried out by hymenopterists Sarah K. Gess and the late Friedrich W. Gess over a period of approximately 40 years (Gess & Gess, 2014;Melin & Colville, 2019).
However, these hymenopterists conducted widespread sampling throughout the western parts of South Africa during this timeframe (Gess & Gess, 2014).
Between the extremes of the Winter rainfall bee hotspot and the Aseasonal rainfall bee hotspot, we see a large proportion of the GCFR as seemingly depauperate of bee diversity.Given that the Winter rainfall bee hotspot and the Aseasonal rainfall bee hotspot largely fall within the GCFR, with similar habitat and processes, it is conceivable that bee diversity is most likely higher in this area, for example, high species richness has been retrieved in well-sampled insects such as butterflies (Mecenero et al., 2013) in this area.This gap may be due to poor sampling (Melin & Colville, 2019) and therefore targeted surveying in this area, particularly in protected areas as seen in recent bee surveys (e.g., Meiners et al., 2019;Messinger Carril et al., 2018), could be used to establish if the Winter rainfall bee hotspot and the Aseasonal rainfall bee hotspot are contiguous across the south-eastern GCFR.
The Early-to-late summer rainfall bee hotspot is dominated by more widespread generalist bee species and shares many of the taxa with other Afrotropical areas (e.g., Eucerini (Eardley, 1989) and Melectini (Eardley, 1991)).The low numbers of endemic and range-restricted species richness in this hotspot, particularly away from the Drakensberg area, confirm this pattern (Figure 2b-d).In contrast, Kuhlmann (2009) using a restricted species dataset and at a coarse resolution ($200 Â 200 km grid square) found the proportion of endemic bee species in the early-to-mid summer rainfall area to be almost equivalent to the winter rainfall area, although recognising that the true pattern is most likely considerably lower in the early-to-mid summer rainfall area and more in line with the results of this study.In the Early-to-late summer rainfall bee hotspot, other insects, such as butterflies (Mecenero et al., 2013(Mecenero et al., , 2015)), dragonflies (Deacon et al., 2020), and dung beetles (Davis, 2002) exhibit high species richness.Similar to bees, these insects demonstrate overall low values of endemism and share numerous widespread taxa with other Afrotropical areas.
As observed in other parts of the world, our bee richness and endemism patterns conform to the general pattern that bee faunas are highest in arid or Mediterranean climates (Michener, 1979;Orr et al., 2021).Across these regions, specific climatic factors and habitat heterogeneity including high floristic diversity rather than area size per se appear the more important variables in determining bee richness (but see Kaloveloni et al. (2018)).This scenario has also been hypothesized for other areas of exceptional bee richness, such as Pinnacles National Park in Central California (Meiners et al., 2019).Our results also show that these areas are not only disproportionally high in bee species richness but also high in endemics and species with narrow geographic ranges.This suggests that rates of bee diversification and specialisation, and species turnover (beta diversity (e.g., Colville et al. (2002)), are associated with the high number of fine-scale floristic and edaphic niches of these areas (see also Cowling et al. (2015)).
In summary, the Mediterranean climate and arid areas in South Africa contain a high density of bee species and a large number of endemic and geographically restricted bee taxa.Although our diversity maps are based on 70% of South African bee species, adding records for missing species would most likely strengthen our hotspots.
These hotspots have frequently been recognised as areas of high richness for a range of insects and plants.This study goes further by illustrating the entomofaunal distinctiveness of these hotspots in terms of disproportionally high numbers of species, endemics and rangerestricted taxa.The strong correlations observed between diversity metrics suggest that species richness is a suitable indicator of endemic and range-restricted bees in South Africa.

CONCLUSION
Understanding bee diversity at a regional scale can give insights into steps in understanding these regional patterns will require an assessment of the underlying processes that shape them (Kaloveloni et al., 2018), particularly in the context of specialisation and the extraordinarily rich flora and the role of historical climatic stability.
Future analysis, using our species distributional models, can also be extended to calculate other patterns of diversity such as betadiversity (species turnover) and functional diversity (e.g., Leclercq et al., 2023;Marshall et al., 2023).Because of the high numbers of range-restricted taxa confined to small numbers of grid cells, habitat specialisation is most likely driving high species turnover and functional diversity in the Winter rainfall bee hotspot and the Aseasonal rainfall bee hotspot.As undertaken for plants (Cowling et al., 2015), comparisons across other global Mediterranean and arid areas of exceptional bee richness would also provide deeper insights into the drivers of bee diversity and the role of habitat heterogeneity, including floristic resources, and historical climate stability.In addition, determining the extent to which the existing protected area networks in South Africa, based predominantly on plant-based data (e.g., Skowno et al. (2018)), adequately cover the bee hotspots warrants further scrutiny.This could be expanded by taking into consideration the vulnerability of the South African bee fauna under several scenarios of global change (Kuhlmann et al., 2012).Finally, refining the hotspot boundaries presented here would require long-term inventories to ground truth these patterns.
richness concentrates into three distinct hotspots (Figures 1-3): Winter rainfall bee hotspot the arid and mesic winter rainfall zone along the west coast and incorporating parts of the Fynbos and Succulent Karoo biomes; Aseasonal rainfall bee hotspot the mesic aseasonal rainfall zone extending along the south coast incorporating parts of the Fynbos, Succulent Karoo and Albany Thicket biomes and encroaching inland towards the drier interior of the very late summer rainfall zone; and Early-to-late summer rainfall bee hotspot a narrow band along the east coast (subtropical) extending inland into the high-altitude montane region of the eastern side of the Drakensberg Mountains, and comprising the north-eastern extent of the summer rainfall zone and including parts of the Grassland and Savanna biomes.

F
I G U R E 1 South African bee richness hotspots at a $11 Â 11 km grid cell resolution: (1) Winter rainfall bee hotspot: along the west coast and incorporating parts of the Fynbos and Succulent Karoo biomes, (2) Aseasonal rainfall bee hotspot: extending along the south coast incorporating parts of the Fynbos, Succulent Karoo and Albany Thicket biomes, (3) Early-to-late summer rainfall bee hotspot: a narrow band along the east coast extending inland into the high-altitude montane region of the eastern side of the Drakensberg Mountains, and including parts of the Grassland and Savanna biomes.
Comparison of areas with similar arid and Mediterranean climates revealed that the San Bernardino Valley (Arizona/ Sonora) holds the highest rank, followed in descending order by Pinnacles National Park (California), Lesvos (Greece), Yosemite National Park (California), and Grand Staircase Escalante National Monument (Utah) (Figure 4b).These rankings surpassed those of all South Africa's biomes.Notably, the Albany Thicket and Fynbos biomes outperformed analogous areas F I G U R E 2 Spatial patterns of South African bee diversity (a-f) at a $11 Â 11 km grid cell resolution.such as Mediterranean Chile and Southwestern Australia (Figure 4b).Additionally, the more arid Succulent Karoo and Desert biomes secured higher rankings than the Mojave National Preserve, San Rafael Desert and Desert Chile.

F
I G U R E 3 Biomes (a) and rainfall seasonality (b) of South Africa.
local scale patterns and can be highly informative for understanding broader global patterns of bee diversity.It can also provide important information on local areas of insect diversity, such as hotspots for conservation, areas of evolutionary interest where certain lineages are concentrated, and in understanding the historical biogeography of a region.The spatial patterns of species richness and endemism in South African bee fauna presented here corroborate but vastly extend the results of previous studies and conform with global patterns of unusually high bee diversity in Mediterranean and arid areas.The next T A B L E 2 Comparisons of number of species density across South Africa's biomes.Comparisons of number of species density across South Africa's rainfall seasons.Comparisons are read left to right, and values indicate the ratio after dividing species density of one season by another.Rainfall season area size (km 2 ) from smallest to largest: All-year (52272.04481);Winter (136589.81);Early-summer (149991.16);Late-summer (257857.34);Very-late-summer (305383.55);Mid-summer (318943.09).