Introduction

In southern Africa, the archaeological record of the Middle Stone Age (MSA), dated to 280,000–50/25,000 YBP (McBrearty & Brooks, 2000), displays variations in technological and cultural attributes and shifts in ecological niches (d’Errico et al., 2017; Mackay et al., 2014; Wurz, 2016). These variations were particularly associated with Marine Isotope Stage 5 (Greenland Stadials and Interstadials 25–20) at ca. 130,000–71,000 YBP (Lisiecki & Raymo, 2005; Railsback et al., 2015; Rasmussen et al., 2014; Shackleton et al., 2003; Wolff et al., 2010). Scholars have explored the role of climatic and vegetation changes as potential catalysts for human adaptability, changing mobility patterns, and increased cultural complexity during this time (Compton, 2011; Mackay et al., 2014; Marean et al., 2020). Reconstruction of environmental conditions is crucial for our understanding of resource procurement strategies, technological and cultural innovations, and site choices made by these early humans. Analyses of micromammals incorporated into archaeological sediments have the potential to elucidate past local environmental conditions (e.g., Andrews, 1990; Avery, 1979, 1982, 2002; Matthews, 2004; Matthews et al., 2009, 2011, 2020; Stoetzel et al., 2011, 2018). Micromammals are referred to here as rodents, shrews, and bats with a live weight of < 200 g. These small mammals are suitable palaeoenvironmental informants as they have small home ranges (usually less than 1-km radius). Some have precise ecological requirements and often exhibit high intrinsic rates of population increase and the need for rapid response to environmental change within fine spatiotemporal scales (Heisler et al., 2016; Korpimäki et al., 2004; Reed et al., 2019). Analyses of modern micromammal samples have demonstrated a close correlation between the relative abundance of species and composition of vegetation substrate near sample sites in Europe, East Africa, and South Africa (Andrews, 1990; Avery et al., 2005; Avery, 1992a; Reed et al., 2019). D. M. Avery has been prominent in establishing a South African methodological record, demonstrating the suitability of micromammals as palaeoenvironmental informants of vegetational and climatic conditions through the application of biodiversity indices, analyses of taxonomic composition in relation to habitat requirements, and relative abundance of species (Avery 1979, 1982, 1987, 1992a, b, 2001, 2007). Her research pioneered the analyses of micromammal assemblages from archaeological sites in South Africa during the past two decades (e.g., Faith et al., 2019; Matthews, 2004; Matthews et al., 2005, 2009, 2011, 2020; Nel, 2013; Nel & Henshilwood, 2016; Nel et al., 2018). Recent research has emphasized the need for improved information about modern habitat requirements and distributions of some micromammal species to better understand past and present ecological niche behaviors (Avery, 2007; Matthews & Nel, 2021; Matthews et al., 2020). This issue is discussed further in this article. Simultaneously, new analytical developments are in the early stages of exploration, such as analysis of stable carbon and oxygen isotopes of micromammal tooth enamel from archaeological contexts in South Africa (e.g., Leichliter et al., 2016, 2017; Williams et al., 2020).

In this article, we present environmental implications derived from analyses of the micromammal assemblage from the MIS 5c–a, ca. 104,000–74,000 YBP sequence at Blombos Cave (BBC), situated on the southern Cape coast of South Africa (Fig. 1; Lisiecki & Raymo, 2005; Railsback et al., 2015; Rasmussen et al., 2014; Shackleton et al., 2003; Wolff et al., 2010). The sequence is divided into the M2 Lower and M3 phases, based on human occupation, and is characterized by on-site production of lithics such as blanks, blades, and flakes using multiple core reduction methods with some variations throughout (Douze et al., 2015). There are indications of ochre processing as evidenced by the liquefied ochre-rich mixtures produced in two abalone shells (Halotis midae), ca. 100,000 YBP (M3 phase; Fig. 1; Henshilwood et al., 2011). Ochre, bone, charcoal, grindstones, and hammerstones form a composite part of this production toolkit. Chemical analyses of the toolkit components and the contents of the containers show a complex process involving heat treatment of bone to extract marrow (fat) which was then mixed with ground ochre powder and charcoal (Henshilwood et al., 2011). We aim to establish a local palaeoenvironmental record for BBC during MIS 5c–a by reconstructing local vegetation and climatic conditions to provide an environmental context for the manufacturing of the ochre toolkits and contribute to the understanding of the site utilization.

Fig. 1
figure 1

Locations mentioned in text (with rainfall zones indicated): BBC—Blombos Cave, BBS—Blombos Schoolhouse, BKP—Byneskranskop, BP—Boomplaas, DHN—De Hoop Nature Reserve, DK—Die Kelders, KEH—Knysna Eastern Head Cave 1, KRM—Klasies River main site, NBC—Nelson Bay Cave, and PP—Pinnacle Point

Full size image

Regional Setting

Blombos Cave (BBC) is a 55-m2 limestone cave situated in the Blomboschfontein Nature Reserve (Hessequa Municipality) on the southern coast of South Africa, approximately 300 km east of Cape Town (Fig. 1). The cave formed due to solution action and wave cutting of the calcarenite and calcrete cliff situated above a basal layer of Table Mountain Sandstone of the Cape Supergroup (Henshilwood et al., 2001). It is 35 m above sea level and about 100 m from the ocean (Henshilwood et al., 2001).

Climate and Vegetation at Present

The climate system in the southern coastal region of South Africa is influenced by both tropical/subtropical and temperate atmospheric systems. The latitudinal location of South Africa, between 22°S and 34°S, is at the intersection between these atmospheric pressure systems (Braun et al., 2017; Chase & Meadows, 2007; Tyson, 1986; Tyson & Preston-Whyte, 2000). The complex interactions of these systems result in the South African rainfall seasonality patterns. The tropical synoptic systems affect the summer rainfall region to the north and east of the country, and the temperate westerly systems bring rain to the winter rainfall region along the west coast. Their interaction can produce significant rainfall events at any time of the year in the aseasonal rainfall region on the southern coast (Braun et al., 2017; Engelbrecht et al., 2015; Roffe et al., 2019; Fig. 1). Based on this seasonal division of rainfall, the regions are known as the summer rainfall zone (SRZ), winter rainfall zone (WRZ), and year-round or aseasonal rainfall zone (ARZ) (Tyson & Preston-Whyte, 2000; Fig. 1). The WRZ receives predominantly frontal-induced rainfall, driven by westerly-derived mid-latitude cyclones from April through September (Philippon et al., 2012; Roffe et al., 2019). In the SRZ, tropical easterly winds bring moisture from the Indian Ocean during summers (March to October), while the winter months are arid, with less than 30% of the annual rainfall occurring during this period (Chase et al., 2015; Crétat et al., 2012). Across the ARZ, rainfall is received from both the winter and summer rainfall systems. Thus, there is no pronounced seasonality; however, there are small peaks in rainfall during April and September (Engelbrecht & Landman, 2016; Engelbrecht et al., 2015).

BBC is situated in the ARZ and, at present, the Blombos area receives approximately 56% of the annual rainfall during the winter half-year (Fick & Hijmans, 2017). The current mean annual precipitation (MAP) at BBC is 537 mm, and temperatures range from a mean daily temperature of 20.7°C to a maximum temperature of 27.0°C in the summer. While during winter, the mean daily temperature is 13.4°C, with a minimum temperature of 6.5°C. The average temperature for the whole year is 17.08°C (Fick & Hijmans, 2017). The landscape near BBC comprises mosaic vegetation, consolidated dunes, numerous rocky outcrops, and patches of loose sandy soil with sparse or no vegetation cover (Table 1). Freshwater is available from two springs 300 and 600 m east of the cave. The current vegetation characteristics within a 10-km radius of BBC are summarized in Table 1.

Table 1 Characteristics of vegetation components currently present within an approximately 10-km radius of BBC (data from Biodiversity Geographic Information Systems (BGIS) by South African National Biodiversity Institute (2012) at http://bgisviewer.sanbi.org/, Mucina et al. (2006), Mucina and Rutherford (2006), Rebelo et al. (2006), and Bergh et al. (2014))
Full size table

It should be noted that the sea level fluctuated during MIS 5. Fischer et al. (2010) proposed that the shoreline would have been 1.45–2.33 km from the site during the M3 phase, with a potential further regression during M2 Lower to between 2.33 km (minimum) and 9.42 km (maximum). Due to the offshore topography, BBC would thus have been situated inland from the coast during some periods of MIS 5, as lower sea levels would have exposed the Palaeo Agulhas plain (PAP) (Cawthra et al., 2020; Fischer et al., 2010; Haaland et al., 2020; Marean et al., 2020). Climatic fluctuations on glacial–interglacial timescales are likely to have impacted the southern Cape region as changes in polar ice volume, and interhemispheric temperature gradients shift the location of atmospheric circulation systems (Braun et al., 2020). Based on climate simulations, using models with different levels of complexity, Broccoli et al. (2006) found cooling in the northern hemisphere to be stronger than in the southern hemisphere during glacial phases, which led to a southward shift of the Intertropical Convergence Zone (ITCZ) and tropical easterlies. Greater frequency in the tropical temperate troughs could have been the cause of more humid conditions in the southern Cape during MIS 4, as Chase (2010) suggested. Speleothem analyses from southern coastal locations indicate that winter and summer rainfall intensity fluctuated considerably during MIS 5–3 (Bar-Matthews et al., 2010; Braun et al., 2019, 2020). Meanwhile, the exposed Palaeo Agulhas Plain would have facilitated a mosaic of known vegetation from the southern Cape area, as Cowling et al. (2020) have modeled the scenario for the Last Glacial Maximum. This likely enabled the establishment of woodland and grassland on fertile, alluvial soils, creating savanna-type vegetation that is not commonly found in the region today. Stratiform winter rainfall has been found to correlate to high abundances of C3 plants, while convective summer rainfall was generally associated with the spread of C4 grasses (Braun et al., 2019). The foregoing illustrates the range of potential climatic and environmental variabilities that could have affected the micromammal population during MIS 5 in the southern Cape region.

Stratigraphy and Chronology

The MSA levels at BBC have been divided into four phases (M1, M2 Upper, M2 Lower, and M3). In this article, we discuss the oldest phases—M2 Lower and M3. We use the Bayesian chronological framework constructed by Jacobs et al. (2020), based on 37 single-grain Optically Stimulated Luminescence (OSL) ages. Although additional ages are available for BBC (Fig. 2), the Bayesian model is the most comprehensive published chronology for the site. The model was created in OxCal (OxCal 4.2.4; Ramsey, 2009; Ramsey & Lee, 2013), modeling each MSA phase as a Phase, separated by Boundaries and ordered into a Sequence informed by stratigraphy (Jacobs et al., 2020). A Transitional Boundary was placed between each Phase, assuming that the BBC sediments have accumulated continuously without any significant breaks (Jacobs et al., 2020). Consequently, the ages of each of the four MSA phases (M1, M2 Upper, M2 Lower, and M3) are constrained by the modeled probability distributions of the Boundaries defining their start/end, and modeled ages are not available for individual layers. The Bayesian modeled OSL chronology is presented in Fig. 2, along with the stratigraphy of the MSA sequence at BBC and published optically stimulated luminescence (OSL) (Jacobs et al., 2020), thermoluminescence (TL) (Tribolo et al., 2006), and uranium-series (U/Th) dates (Henshilwood et al., 2011), and their relationship to the MIS (sub) stages.

Fig. 2
figure 2

The Middle Stone Age stratigraphy and dating sequence for Blombos Cave (courtesy of Magnus M. Haaland)

Full size image

Jacobs et al. (2020) estimate the transition from M2 Upper to M2 Lower at 76.0 ± 5.6 kya (95.4% probability, total uncertainty), and the M2 Lower to M3 at 82.5 ± 6.2 kya. The M2 Lower thus corresponds to MIS 5b/5a (93–74 kya) (Lisiecki & Raymo, 2005; Railsback et al., 2015; Rasmussen et al., 2014; Wolff et al., 2010). The M2 Lower layers are characterized by low frequencies of anthropomorphic material, although some basin hearths and flakes, blades, and cores of quartz, silcrete, and quartzite are present (Thompson & Henshilwood, 2011). The reduced intensity of anthropogenic material defines the transition from M3 to M2 Lower in the M2 Lower phase (Table 2; Haaland et al., 2020). The start of the M3 phase occurs at 97.7 ± 7.6 kya (95.4% probability) based on the Bayesian age estimates (Jacobs et al., 2020). The M3 phase corresponds to MIS 5c (104–93 kya) and MIS 5b (93–84 kya) (Lisiecki & Raymo, 2005; Railsback et al., 2015; Rasmussen et al., 2014; Wolff et al., 2010). The M3 phase (CH–CR) has several culturally sterile levels rich in microfauna, interspersed with levels containing a high frequency of ochre pieces and food waste, especially shellfish. There are eight deliberately engraved ochre slabs from the M3 phase (Henshilwood et al., 2009), in addition to the ochre toolkits previously mentioned. However, only CH/CI, CJ, and CJh1 in the upper M3 phase have high human occupation intensity deposits, while CH and the levels below CJ are low-intensity deposits (Table 2).

Table 2 Division of phases, layers, and squares sampled from BBC for this study (anthropogenic material frequency is based on Thompson & Henshilwood, 2011, Nel, 2013, and Haaland et al., 2020)
Full size table

Material and Methods

Material

The materials for this study are from layers CG–CR, excavated during 2000–2010. Both the M2 Lower and M3 phase have high densities of fossil micromammals. Therefore, sub-quadrants containing a continuous sequence of M2 Lower and M3 layers were selected to obtain representative samples (sample columns) (Table 2). Materials from 16 layers were analyzed, although layers CG–CGAA, CGAB/h1, and CQ also contained sub-layers that were grouped for this study (Table 2).

Methods

Establishing the Accumulator and Post-Depositional Processes

The micromammal remains were extracted from sieved fragments larger than 1.5 mm. Limb bones and cranial elements were examined under a × 40 binocular microscope. Taphonomic analyses are necessary to establish the predator(s) or accumulator agent(s) of the micromammal assemblage and the chemical and mechanical post-depositional processes. The analyses follow the standard methodology described in detail elsewhere (Andrews, 1990; Nel & Henshilwood, 2016) by assessing micromammal fragmentation patterns, digestive marks (grade and frequency), and post-depositional alterations, such as polishing by water, weathering, trampling, rootmarks, and burning, among others.

The various proportions in which the skeletal elements of prey occur in a predator assemblage can distinguish the different predator types (Andrews, 1990). The skeletal representation is Ri = Ni/(MNI × Ei) where Ri is the relative abundance of the element i, Ni is the number of elements in the assemblage, MNI is the minimum number of individuals, and Ei is the number of elements in the prey skeleton (Andrews, 1990). We quantified the fragments of limb bones, represented by humeri (expressed in percent), to provide information of potential destruction by predators and post-depositional trampling. The average relative abundance (ARA) of skeletal elements (means of the relative abundances for all skeletal elements except isolated teeth) was also calculated (Andrews, 1990). ARA would be 100% if there were no loss of elements from a skeleton, and low percentages would indicate loss of elements due to digestion, fragmentation, or post-depositional transportation by water or wind.

Protocols of previous research on micromammals in South Africa (Avery, 2002; Matthews, 2004; Matthews et al., 2009, 2011) and new research (Fernández-Jalvo et al., 2016) led to minor changes in digestion classification compared to the standard developed by Andrews (1990) especially with regards to differentiating light digestion caused by barn owls (Tyto alba) and light to moderate digestion caused by spotted eagle owl (Bubo africanus). Rodentia incisors, both in situ and isolated, humeri and femora, were microscopically analyzed for predator-related digestion. Their suitability for digestion analysis is further described in Nel and Henshilwood (2016).

The Taxonomic Record

Taxonomic identification was based on molars, mandibles, maxillae, and dental morphology, following standard methods (Nel, 2013) with the aid of comparative specimens from Iziko South African Museum in Cape Town and identification keys developed by Avery (1979) and de Graaff (1981). Specimens not identified to species level were assigned to either family or sub-family, while post-cranial elements were all registered as micromammals. No further identification to species level was made except for obvious elements such as limb bones of bats, shrews, and golden moles. Systematic classification followed Wilson and Reeder (2005).

Palaeoecological Study

Biodiversity indices based on micromammals can be used to make inferences about local vegetation structures. Complex structures, which contain a range of plant species and comprise various habitats, are associated with the general abundance and diversity of micromammal species. In contrast, “simple” structures have an inverse relationship (Cuenca-Bescós et al., 2009). Species richness, diversity, and dominance of micromammal species thus serve as measures of vegetation structure near BBC. Species richness is expected to be greater where seasonal variability in thermal, energy production (i.e., vegetation), and precipitation regimes are limited (Andrews & O'Brien, 2000; Belmaker & Hovers, 2011). In southern Africa, small mammal species richness correlates strongly with the diversity of plant species, while the correlation to seasonal distribution of rainfall and maximum monthly precipitation (MMP) is less apparent but still significant. Thermal variation has a weak correlation, particularly concerning small mammals (Andrews & O'Brien, 2000).

Species richness (number of taxa in each layer) is affected by sample size. Rarefaction curves were calculated to investigate whether variation in taxonomic abundance results from different sample sizes (Belmaker et al., 2016; Hammer & Harper, 2006; Magurran, 2013). The curves were based on the minimum number of individuals (MNI) per layer.

General diversity considers both the number of taxa and the relative frequency (evenness of representation) of each taxon. General diversity was assessed by the Shannon–Wiener index, H, where H =  − Σ Pi (ln Pi), where Pi is the proportion (P) of taxon, i, in the assemblage. The Simpson index assesses the probability that two randomly picked individuals are of the same species, and evaluates if there is a dominance of single/few species (Hammer & Harper, 2006). The Simpson’s index of dominance, D, is given as D = Σ (p2i), where pi = ni/n (the proportion of species i). The value will be close to 1 if there is a single dominant taxon in the sample. The results of the general diversity and dominance indices were subject to a t test. All statistical calculations were completed with Paleontological Statistics (PAST) (Hammer et al., 2001).

The taxonomic habitat index (THI) developed for BBC is described in detail in Nel and Henshilwood (2016) and was also calculated for the M2 Lower and M3 small mammal assemblage. The index shows variations in local vegetation and substrate characteristics through the analyzed sequence. A cumulative index is developed by combining habitat indications of all species occurring in a micromammal assemblage (Andrews, 1990; Andrews & Evans, 1983). The index comprises categories based on specific habitats of extant micromammals in the study assemblage and characteristics of current vegetation types in the Blombos area (Table 3). The categories are described in detail in Nel and Henshilwood (2016) and accommodate local vegetation and substrate. It is based on present vegetation characteristics as summarized in Table 1. The weighting of species followed the method presented in Nel (2013), Nel and Henshilwood (2016), and Nel et al. (2018).

Table 3 The weighting of taxa is based on information of distribution, preferred habitats, and habits of the micromammal taxa*
Full size table

BBC is situated in the all-year rainfall zone (aseasonal rainfall), which means it is in an area that straddles the transition between winter rainfall to the west and summer rainfall to the east (Fig. 1). Potential changes in rainfall seasonality and aridity were examined using a backdrop of modern data derived from 123 micromammal assemblages collected from barn owl (Tyto alba) roosts in South Africa (Faith et al., 2019). The roost locations encompass a broad range of precipitation regimes (19–91% summer rain), occurring in the summer rainfall zone (n = 60), the winter rainfall zone (n = 41), and the aseasonal rainfall zone (n = 22) (Faith et al., 2019). In addition, based on the classification scheme for aridity developed by the United Nations Environmental Programme (UNEP, 1997), sites are characterized by climates that are arid (n = 41), semi-arid (n = 52), dry sub-humid (n = 28), and humid (n = 2) (Faith et al., 2019). For each of the locations, Faith et al. (2019) extracted data on rainfall seasonality from WorldClim global climate database (Hijmans et al., 2005), and aridity was taken from the Global Aridity Index database (Trabucco & Zomer, 2009; Zomer et al., 2007, 2008). Faith et al. (2019) have used canonical correspondence analysis (CCA), a constrained ordination approach in which faunal assemblages are ordinated relative to one or more known environmental variables, and the species abundances are considered a response to the environmental gradient (Legendre & Legendre, 2012), to express the ordination of the 123 locations.

Our application of the methodology follows Faith et al. (2019), and relative abundances from the BBC phases (Nel & Henshilwood, 2016) are added as supplementary points in the CCA. The supplementary points are projected in the same multivariate space as the modern assemblages, but they do not influence the original ordination. We amended Faith et al.’s (2019) method by excluding Cape dune mole rat (Bathyergus suillus) as it is likely that humans were the main accumulators of these species at BBC (Henshilwood, 1997; Nel & Henshilwood, 2016). The CCA displays taxonomic composition of the modern assemblages against three known environmental variables; precipitation seasonality as indicated by two variables; % warm-quarter and % cold-quarter rainfall; and an aridity index (AI) (Faith et al., 2019). The aridity index (AI) is expressed as AI = MAP/MAE, where MAP is mean annual precipitation and MAE is mean annual evapotranspiration. Although MAP is strongly correlated with AI in Faith et al. (2019) (r = 0.958, p < 0.001), the inclusion of AI is based on the notion that many climate proxies, including flora and fauna, are more directly influenced by moisture availability (humidity/aridity) rather than the total amount of precipitation (e.g., Chevalier & Chase, 2016). Following Thompson and Henshilwood’s (2011) and Haaland’s (2018) notations of occupation phases at BBC, the M3 phase was divided into M3 Upper (CH–CJ—high-intensity occupation) and M3 Lower (CK–CR—low-intensity occupation) for the application of this methodology.

Results

Taphonomy

Fragmentation and Skeletal Element Abundance (SEA)

Based on the analyses of breakage of humeri (Table 4), both M2 Lower and M3 have extensive fragmentation of limb bones where only 4.9–10.1% of humeri are complete. A high frequency of breakage was also observed for the cranial material (Nel, 2013). This is further illustrated by low to moderate percentages of average relative abundance of all skeletal elements, ranging from 25.9% in CGAB/h1 to 38.5% in CP (Table 4). Skeletal element abundances (SEA) are shown in Fig. 3. In CG–CGAA and CGAC, incisors are best represented, while humeri are most frequent in CGAB. Vertebrae, pelves, and scapulae generally are the least represented elements in all M2 Lower layers. In the M3 phase, incisors are the best-represented element in CH, CH/CI2, and CR. Humeri are most abundant in 9 of 13 layers, while in CP tibiae are most frequent. Similar to M2 Lower, vertebrae, scapula, and pelves all have low representation. The distribution pattern of SEA could indicate a mammalian carnivore as being the predator (Fig. 3). However, the SEA conforms to the site-typical pattern as seen in the younger phases (Nel & Henshilwood, 2016). Thus, it is likely that post-depositional processes affected SEA and ARA, which will be further discussed below.

Table 4 Recorded breakage of micromammal humeri throughout the studied sequence at BBC*
Full size table
Fig. 3
figure 3

Skeletal element abundance for a modern samples; b M2 Lower and M3 phases at BBC; c M3 Upper; and d M3 Lower [barn owl (BO), spotted eagle owl (SEO), white-tailed mongoose (WTM), small-spotted genet (SSG)]

Full size image

Digestion

In M2 Lower, 37.9% of incisors were digested in CG–CGAA, 48.0% in CGAB, and 39.8% in CGAC (Table 5). The majority of the observed digestion was within the light category. Heavy digestion was visible on 1.7–2.9% of the incisors, while extreme digestion was only present in CG–CGAA with 0.7%. Digestion on humeri varied between 46.7% and 54.1%. Femurs sustained more gastric etching, with 69.7–83.7% of the bones being digested.

Table 5 The percentage representation of digestion on incisors (combined in situ and isolated), femurs and humeri (the division of incisor digestion classes is based on Andrews (1990))
Full size table

Digestion of incisors in the M3 phase displays a relatively similar pattern throughout the phase, except for the two oldest layers, CQ and CR. For CH–CPA, the total amount of digested incisors range between 50% and 60.8%, all with the majority of the digested elements falling in the light digestion class. In comparison, CQ and CR have less frequent digestion with 37.5% and 35.9%, respectively. CPA has a relatively high percentage (18.4%) of incisors with moderate digestive marks present; the other layers range from 5.9% to 14.8%. Humeri were generally less digested than femurs, as the frequency of humeri digestion spans from 58.7% to 85.9%, while femurs range from 83% to 95.7%.

Post-Depositional Modifications

Etching, pitting, and rounding of breaks were the most frequent post-depositional modifications observed (Table 6). Rounding was particularly frequent on humeri in CG–CGAA (21.5%) and CGAB/h1 (24.5%) in M2 Lower, and in CH/CI2 (26.1%), CJ (29.9%), CK (36.9%), CL (21.6%), and CPA (30.6%). Voorhies categories have excluded water transportation as a cause of rounding (Nel, 2013). Pitting and etching occur in low ratios on the limb bones and incisors (Table 6), except for layer CPA, where 46.8% of humeri and 15.5% of femurs were pitted. In the same layer, 23.3% of incisors, 8.5% of femurs, and 47.7% of humeri had sustained etching. Burnt elements were not frequent in the analyzed sequence, except the hearth feature in sub-layer CGABh1 (here combined with CGAB but analyzed separately in Nel, 2013) where 27.9% of post-cranial (n = 365) and 26.9% of cranial (n = 364) elements were burnt.

Table 6 Percentage of etched femurs, humeri and incisors, rounding of breaks and pitting on humeri and femurs in M2 Lower and M3
Full size table

Taxonomy

The study sample contained 1,477 minimum number of individuals from 23 different species, including golden moles, shrews, mice, vlei rats, gerbils, and bats (Table 7). There were three taxa of Chrysochloridae, four species of Soricidae, thirteen species from the Muroidea superfamily, and three species of Chiroptera.

Table 7 Taxonomic percentage representation (relative abundance in percent), species richness, Shannon–Wiener index, H, and Simpson dominance, D, for the M2 Lower and M3 phases at BBC (the variation in the three latter values is also illustrated in Fig. 4)
Full size table

Chrysochloridae, Soricidae, and Chiroptera

Based on the Bayesian model developed by Jacobs et al. (2020), the transition from M2 Upper to M2 Lower is dated to 76.0 ± 5.6 kya (95.4% probability, total uncertainty), and the M2 Lower to M3 at 82.5 ± 6.2 kya. Golden moles (Chrysochloridae) are present throughout the M3 phase, corresponding to MIS 5c/5b, although absent in the uppermost M3 layer, CH, while they occur more sporadically during MIS 5b/a in M2 Lower before they disappear in the transition to M2 Upper (MIS 5a) (Table 7; Nel, 2013; Nel & Henshilwood, 2016). Although the golden moles are generally represented in low numbers (< 2.8%), Chlorotalpa duthieae is represented most frequently, especially in M3. It is present in all layers except in CH, CK, and CPA. Amblysomus hottentotus is present in MIS 5a/b and MIS 5b/c in layers CG–CGAA, CH/CI, CJh1, and CL. Chrysochloris asiatica is sporadically represented, present only in the M3 phase (97.7 ± 7.6 to 82.5 ± 6.2 kya). Amblysomus hottentotus is currently present in the grassland and woodland habitats of South Africa, especially in the Eastern Cape and KwaZulu-Natal Provinces (Bronner & Mynhardt, 2015). The species also occurs in afromontane forests, marshes, and coastal dune sands, and it extends to the fynbos and Nama-Karoo biomes in the southern limit of their range (Bronner & Mynhardt, 2015). Chlorotalpa duthieae is associated with alluvial sands and sandy loams in southern Cape Afrotemperate forests, particularly the coastal platform and scarp forest patches (Bronner, 2015) along the south coast in a limited area ranging from George to Uitenhage (Fig. 1). They are previously reported in archaeological assemblages on the south coast, within their current distribution area at Nelson Bay Cave (from 23 kya) and about 76 km north of their distribution range at Boomplaas (from 65 kya) (Avery, 1982; Faith et al., 2019). Chrysochloris asiatica is presently associated with a wide variety of vegetation types in a region broadly corresponding to the winter rainfall zone (Fig. 1), although the species’ primary necessity seems to be dry, sandy soils (Bronner, 2015; Skinner & Chimimba, 2005). According to Skinner and Chimimba (2005), the species restrict their breeding to the wet winter months in the Western Cape.

The shrews (soricids) have a continuous presence in the sequence and range from 17.2 to 32% representation. Crocidura flavescens is either absent or represented with one individual per occurrence throughout the sequence. In contrast, Myosorex varius is the most frequent species (together with Otomys irroratus) in two of the mid-M3 phase layers, with 22% in CN/CO and 18.8% in CJ, while ranging in frequency from 5.9% to 17.3% throughout the rest of the sequence. Suncus varilla is also present, in low to moderate frequencies, while Crocidura cyanae occurs throughout the sequence in all layers except toward the M3 Upper (CJ and CH/CI) and at the end of the M2 phase (CGAC), corresponding to the M3 to M2 phase transition at 82.5 ± 6.2 kya.

C. flavescens has been associated with mean annual precipitation (MAP) between 500 and 750 mm (Meester, 1963; Skinner & Chimimba, 2005). Although it may be more tolerant of arid areas than previously assumed, as it was found in modern owl roosts on the arid west coast of South Africa where MAP is between 150 and 275 mm (Avery et al., 2005; Matthews et al., 2011). Avery (1992b) has suggested that C. flavescens may utilize moisture from coastal fogs and underground water to meet its required moisture needs, based on her recording of the species at Spoeg River, a site dated to ca. 1.9 kya, with a current MAP of 50–100 mm (Avery, 1992b). M. varius is considered a generalist species, able to adjust its reproductive cycle to suit different environmental conditions (Baxter, 2005). Hence, the species is not associated with any specific habitat; it occurs in all the biomes of South Africa (Cassola, 2016). S. varilla has characteristics similar to M. varius, although it occurs in less abundance at archaeological sites (Matthews et al., 2020). C. cyanea occurs in a wide variety of montane grasslands and temperate and subtropical forests across southern Africa (Baxter et al., 2016). Bats (Chiroptera) occur in low frequencies (1.0–3.4%) in all layers except for CK. Some may be from natural occurrences, as there are no signs of digestion and bats would likely inhabit the cave, while others have digestive traces and would thus be occurring in the cave due to predators.

Muroidea

Gerbilliscus afra is endemic to the Western and Northern Cape provinces and has low, but rather stabile, relative abundances (1.7–5.7%) during M3 and M2 Lower. Gerbilliscus afra is associated with sandy soils necessary for burrowing. The mesic associated Dasymys incomtus/capensis is present in low numbers, although absent during some periods in MIS 5c/b (layers CR, CM, CL, CJ, and CH/CI2). It is associated with reed-beds, semi-aquatic grasses, and marshes and has a limited, yet persistent, occurrence at BBC (Monadjem et al., 2015; Pillay, 2003; Skinner & Chimimba, 2005). Acomys subspinosus and Myomyscus verreauxii have strong fynbos affiliations (Biccard & Midgley, 2009; Palmer et al., 2017; Relton et al., 2017) and range between 0.9% to 5.3% and 0.4% to 6.3% in representation, respectively. A. subspinosus is fynbos endemic, often found in rocky areas on mountain slopes in fynbos vegetation (Palmer et al., 2017), while M. verreauxii is associated with both lowland and montane fynbos vegetation (Relton et al., 2017). Both species decline in representation after 82.5 ± 6.2 kya (M2 Lower).

Climbing mice (Dendromus spp.) are present throughout (1.8–8.0%), although they decline in general abundance after M2 Lower. The two species identified in the assemblage are Dendromus melanotis and Dendromus mesomelas. D. melanotis occurs in grasslands and savanna, although common in a variety of habitats (Child & Monadjem, 2016a). D. mesomelas is a complex species; it has been recorded in forest/grassland mosaic habitats in South Africa (Child & Monadjem, 2016b). Both climbing mice are primarily associated with tall grasses (Skinner & Chimimba, 2005). In contrast, Steatomys krebsii appears in low frequencies (0.9–1.6%) during the upper M3 phase (CH/CI) and in M2 Lower at BBC. It is associated with dry, shrubland Mediterranean-type vegetation (Monadjem & Schlitter, 2016). Mystromys albicaudatus, a Highveldt grassland specialist, marginally occurs in fynbos and succulent Karoo-type vegetation (Avenant et al., 2019) and is only present in one layer, CH/CI. Mus minutoides, a versatile species with an extensive range throughout Sub-Saharan Africa (Monadjem et al., 2015), is present in low frequencies in all layers except CH/CI2.

Vlei rats (Otomyinae) are most abundant in the assemblage and range in frequencies from 30% to 51.5%. At least three vlei rat species are present—Otomys irroratus, O. karoensis, and O. laminatus. Otomys irroratus is the most abundant species in all layers (14.3–28.2%) and is associated with fynbos and Albany thicket biomes (Monadjem et al., 2015) particularly in wetlands (Taylor & Baxter, 2017). It generally occurs in dense vegetation cover and higher moisture contents, such as grasslands and marshes (Taylor & Baxter, 2017). Otomys karoensis co-occurs with O. irroratus, however, in significantly lower frequencies and is absent in layer CH/CI. O. karoensis is mountain fynbos, thicket, and grassland endemic species, common in Restio-dominated fynbos, predominantly on rocky mountain slopes and dense fynbos patches (Baxter et al., 2017). The laminate vlei rat (Otomys laminatus) associated with moist grassland and shrubland (Taylor & Baxter, 2019) occurs sporadically in M3 and M2 Lower. Its present distribution is recorded in two separate areas, KwaZulu-Natal and Mpumalanga Provinces and Paarl area in the Western Cape Province. Skinner and Chimimba (2005) have noted that there may be deficits with regard to the reported distribution areas for the species. Matthews et al. (2020) found the species occurring at six modern owl pellet locations in the forest, grassland, and thicket biomes along the Garden Route (Wilderness to Plettenberg Bay) in the southern Cape and also in 23 k-year-old layers at Knysna Eastern Heads Cave 1, while Faith et al. (2019) report it occurring at Boomplaas from ca. 65 kya, Nelson Bay Cave at 23 kya, and Byneskranskop at 17 kya. Even though it was not recovered from modern pellet collections in De Hoop Nature Reserve or BBC (Avery et al., 2005.), it seems to be a prevalent taxon in the southern Cape 100 kya. Rhabdomys pumilio is present in all layers with a consistent relative abundance, although with peaks in representation in M3 Lower (CR, 22.0%) before 97.7 ± 7.6 kya (MIS 5c), and late M2 Lower (CG–CGAA, 18.5%, MIS 5b/a). The species is a niche generalist, highly socially flexible, switching from group to solitary living in accordance with population density. Thus, it is able to adjust to environmental fluctuations and therefore occurs frequently in the archaeological micromammal assemblages on the south coast (Du Toit et al., 2012; Ganem et al., 2012; Matthews et al., 2020; Nel et al., 2018; Schradin et al., 2012).

Palaeoecology Measures

Species Richness and Rarefaction

Species richness for each layer is presented in Table 7 and Fig. 4a. The highest number of taxa, 20 species, occurred toward the transition between M2 Lower and M2 Upper in CG/CGAA during MIS 5b/a. In contrast, the uppermost layer CH toward the M3/M2 Lower transition at 82.5 ± 6.2 kya had 12 species. Rarefaction curves for the M2 Lower and M3 layers generally overlap, although species richness varies (Table 7; Figs. 4a and 5). The rarefaction results were standardized to the smallest sample based on MNI (CK—species richness = 14, MNI = 53) and compared statistically by a t test (Hammer & Harper, 2006).

Fig. 4
figure 4

Biodiversity indices for the analyzed layers within the M2 Lower and M3 sequence. a Species richness, S; b Shannon–Wiener, H; and c Simpson dominance, D, with SDs expressed as confidence bands

Full size image
Fig. 5
figure 5

Individual rarefaction curves of the statistically significant layers (CG–CGAA, CGAC, CH, CH/CI, CH/CI2, CJh1, CN/CO, CPA, and CR) that are comparable to the smallest sample based on MNI (CK)

Full size image

Layers CG–CGAA, CGAC in M2 Lower and CH, CH/CI, CH/CI2, CJh1, CN/CO, CPA, and CR in the M3 phase are statistically different when standardized to the smallest sample, CK, from the mid-M3 phase (Table 8 and Fig. 5). The observed variations are thus not due to differences in sample size. Rather, the factors influencing the species richness in these layers could result from environmental productivity (Belmaker & Hovers, 2011). All other layers from the assemblage are within expected variations of sample size (Table 8). However, it should be noted that CK has a high species richness comparable to sample size and thus may represent an outlier in this instance. Other biodiversity indices will shed light on this question.

Table 8 Standardization of individual rarefaction curves: results of t test and permutation test (n = 9,999) for the micromammal assemblage from the M2 Lower and M3 phases (significance at p < 0.05 is in bold)
Full size table

Diversity and Evenness

The least diverse micromammal assemblage based on the Shannon–Wiener index (H = 1.82) occurs just after the M3 to M2 Lower transition (CGAC), while CL has the highest diversity where H = 2.42 (Table 7 and Fig. 4b). There is relatively low dominance of single species throughout the assemblage, as shown by the Simpson index. However, CGAC also has the highest dominance index (D = 0.24). In contrast, CL and CK, with D = 0.12, have the lowest dominance of any given taxa throughout the analyzed sequence (Table 7 and Fig. 4c). Diversity and dominance indices were statistically compared for each layer by a diversity permutation test for equality where p < 0.05 (Table 9; Hammer & Harper, 2006). The statistically different layers (p < 0.05) are marked in bold italics in Table 9. Layer CGAC is statistically different from most other layers in the M2 Lower and M3 phases (Table 9).

Table 9 P values of the t test based on the results of the Shannon–Wiener H (diversity) and Simpson D (dominance) indices (significance at p < 0.05 is in bold)
Full size table

Taxonomic Habitat Index (THI)

The THI shows that grass is the most abundant vegetation component in the local area throughout the MIS 5c–a, although there are variations in the abundance (Table 10 and Fig. 6). Moist grass has the highest representation but varies throughout the sequence (18.5–28.4%). There are lower abundances during MIS 5c/b in layers CR, CM-CL, and CH/CI2. However, a rapid increase of 9.9% occurred from CL to CK. Dry grass is abundant throughout with only 4% variation. The bush/forest component is relatively consistent as well (variation 3.1%). Shrubland fluctuates, with the highest level of abundance occurring before the 97.7 ± 7.6 kya transition to layer CR (16.2%), and there is an overall decline in M2 Lower (MIS 5b/a). Both coastal scrub and the rocky component are relatively consistent, with a spike in coastal scrub in CK. Sandy areas expand during the upper M3 phase, before a decline in layer CH, and then increase in M2 Lower.

Table 10 Reconstructed microhabitat at BBC based on the taxonomic habitat index (the results are given as percentage distribution of the various vegetation and substrate types and cover the entire Middle Stone Age sequence at BBC
Full size table
Fig. 6
figure 6

The percentage representation of vegetation components from the M2 Lower and M3 sequence (for percentage representation by layers, see Table 9)

Full size image

Rainfall Seasonality and Aridity

The ordination of modern and archaeological sites and species is plotted along two axes in the CCA (Fig. 7). Summer–winter rainfall ordinates on the x-axis are divided into ARZ, aseasonal rainfall zone (ARZ), summer rainfall zone (SRZ), and winter rainfall zone (WRZ). The arid-humid conditions are reflected on the y-axis (humid = AI > 0.65; dry sub-humid = 0.5 < AI < 0.65; semi-arid = 0.2 < AI < 0.5; arid = 0.05 < AI < 0.2).The modern sites are positioned on the graph (Fig. 7) according to precipitation zones (WRZ, ARZ, and SRZ) and UNEP aridity classification (arid, semiarid, dry sub-humid, and humid). For clarification, only the archaeological phases and modern sites discussed in the text are indicated by names. The first axis, x, has 54.32% of the inertia and reflects the precipitation seasonality, while the second axis, y, has a 45.63% value and reflects the aridity index (Fig. 7; Supplementary Tables 1, 2, and 3). For the first axis, negative values are associated with a higher percentage of warm quarter rainfall (summer), while positive values are higher percentages of cold quarter rainfall (winter). As inferred from site locations and species, arid conditions are displayed as positive values along the second axis, while negative values indicate sites and species associated with more humid conditions. Based on the present-day site data, Faith et al. (2019) have demonstrated a significant inverse correlation between these two variables (r =  − 0.880, p < 0.001).

Fig. 7
figure 7

Canonical correspondence analysis (CCA) showing ordination of the 123 modern locations (Faith et al., 2019) with BBC phases as supplementary points

Full size image

For comparative purposes, previously published data from the two younger phases at BBC, M2 Upper (MIS 5a) and M1 (MIS 5a/4), are included (Nel, 2013; Nel & Henshilwood, 2016). The BBC sequence then spans the MIS 5c to MIS 5a/4 transition. The CCA ordination scores for species, modern locations, and archaeological phases at BBC are included in Supplementary Table 1. The information on the 123 locations and their aridity score and rainfall seasonality is given in Supplementary Table 2; the weighting of the micromammal species (precipitation seasonality and AI) is available in Supplementary Table 3. The M3 Lower, M3 Upper, and M2 Lower phases, corresponding to MIS 5c–a, are close to each other on the CCA and nestled with modern locations in the lower range of dry sub-humid on the aridity scale. Furthermore, the phases are close to sites that presently receive 30.7% (De Hoop Nature Reserve), 34.2% (Groot Hagelkraal), and 39.3% (Die Kelders) rainfall in the winter quarter (Faith et al., 2019: 4; Supplementary Tables 1 and 2).

In comparison, the Blombos area receives 27.6% at present (Fick & Hijmans, 2017, Supplementary Table 2). These suggest that there may have been slightly more precipitation during winter quarters in MIS 5c–a than in the present. Interestingly, the M2 Upper and M1 phases (MIS 5a/4) display greater climatic divergence from the lower phases. Data for the M1 and M2 Upper have been presented in Nel and Henshilwood (2016). Thus, MIS 5a/4 (M2 Upper and M1) were likely more humid and with rainfall seasonality shifting to become more aseasonal. Consequently, there was a reduction in winter-quarter precipitation compared to the older BBC phases (Fig. 7). This does not mean that there was overall more annual precipitation, but rather there may have been a seasonal shift in moisture availability (see Faith et al., 2020; Thackeray, 2020). Thus, the distribution of rainfall throughout the year seems more evenly dispersed. The increase in humidity suggests more moisture availability in the younger MSA phases of MIS 5a/4 than MIS 5c-a (M2 Lower and M3 phases) (Fig. 7).

Discussion

Establishing the Predator—Fragmentation and Digestion

It is essential to establish the predator responsible for accumulating the micromammals in order to assess any potential bias that may affect the interpretation of the palaeoenvironmental data derived from the assemblage. For instance, if the predator is known to only hunt particular prey species, this would not reflect an appropriate sample of the living micromammal community in the area and cannot be used for palaeoenvironmental reconstruction. Owl species are often the main accumulator of an assemblage in the South African context, especially in caves and rockshelters along the southern coast (Matthews et al., 2011, 2020; Nel & Henshilwood, 2016; Nel et al., 2018). The owls would have used these sites for roosting. Owls, such as spotted eagle owl (Bubo africanus) and barn owl (Tyto alba), often return to the same location over generations, allowing for a buildup of regurgitated pellets; a result of their inability to digest fur and bone. The pellets disintegrate over time, leaving bones to become part of the sediments at the sites. This line of evidence is based on modern analogs and indicates that micromammals from owl pellets are accurate proxies for understanding present vegetations (e.g., Heisler et al., 2016). Reed (2005), for example, found spotted eagle owls and barn owls produce similar taphonomic overprints. Both barn owls and spotted eagle owls are opportunistic and hunt a similar range of micromammal taxa on the south coast, with a preference for Otomyinae, Gerbillinae, and Soricids (Avery et al., 2005; Matthews et al., 2011, 2020). The co-occurrence of these two predator species is known from both Pinnacle Point and Klasies River (Matthews et al., 2011, 2020; Nel et al., 2018).

The fragmentation of the micromammal assemblage is a combination of destruction by predators and post-depositional processes. Thus, the skeletal element abundance pattern is not viable for accurate comparison of modern predator patterns since these are based on pristine pellet and scat assemblages. This is a known problem when analyzing archaeological assemblages (Andrews, 1990; Fernández-Jalvo et al., 2016; Matthews, 2004; Nel & Henshilwood, 2016; Nel et al., 2018). It highlights the need to revise skeletal element abundance based on modern experimental analogs where the material has been subject to replicated trampling and other site formation processes.

The digested incisors in the assemblage are mainly within the light and moderate categories. The digestion pattern is similar to modern comparative samples of spotted eagle owls (Andrews, 1990; Fernandez-Jalvo & Avery, 2015; Fernández-Jalvo et al., 2016; Matthews et al., 2011; Reed, 2005). However, in modern samples from spotted eagle owls, 25–26% of the incisors have moderate digestion (Matthews et al., 2011), while the assemblage from BBC has lower ranges, 5.4–18.4%. It should be noted that there is also a discrepancy in digestion frequency reported by different researchers for barn owl (Tyto alba) assemblages. Andrews (1990) found digestion on 8% of incisors from Makapansgat and Boomplaas in South Africa, and Stratton and Barton Turf in the UK. Also, Matthews et al., (2011, p. 217; Fig. 3) observed digestion on over 60% of incisors from modern barn owl pellets collected in South Africa. However, the abundance of prey available, how hungry the predator is, and the predator’s age can affect the digestion intensity and frequency (Fernández-Jalvo et al., 2016). The possible occurrence of both predators at different time intervals is thus likely. The spotted eagle owl was the main predator but not solely responsible for the accumulation. In particular, CQ and CR in M3 and CGAC and CG–CGAA in M2 Lower could have had another primary predator leaving lighter digestion on prey elements, and the contribution of barn owls may be greater than currently estimated for this analyzed assemblage. The large number of micromammals present in the M2 Lower and M3 indicate multi-year roosting in most of the layers sampled. Thus, both spotted eagle owl and barn owl could have used the cave for roosting for prolonged periods without occupying it concurrently.

Post-Depositional Processes and Intra-Site Spatial Signals

The observed rounding is consistent with rounding caused by abrasion, as demonstrated on modern samples by Fernandez-Jalvo and Andrews (2016, p. 179). The abrasion actions are likely caused by several actions: bioturbation, sediment and water movement, and trampling. There is no evidence of an increased abrasion rate during periods of intensified human occupation (Haaland, 2018; Thompson & Henshilwood, 2011). Humeri were more frequently abraded than femurs; the reasons for this are unknown and remain a topic for further taphonomic investigation.

The burnt elements in M2 Lower come from CGAB/h1. This situation is expected because CGABh1 is a localized hearth feature within layer CGAB. The lack of burnt elements in M3 is probably due to the low intensity of human occupation (Haaland, 2018; Thompson & Henshilwood, 2011). Although the upper layers of M3 have a greater frequency of hearths and human occupation intensity compared to the rest of the analyzed sequence (Haaland et al., 2020), this is not visible with regards to accidental burning of the micromammal assemblage. The general frequency of burnt elements is lower than the high-intensity anthropogenic deposits in the M1 and M2 Upper phases (Nel, 2013; Nel & Henshilwood, 2016). The lack of burnt samples in Upper M3 could be a function of sample location. Charred bones may be found in other cave areas where several hearth features have been identified during excavations (van Niekerk, personal communication, December 8, 2020, Bergen) or through micromorphology analyses (Haaland et al., 2020).

The pitting observed in M2 Lower, especially on ulnae and other limb bones in CGAB/h1, seems to be confined to one side of the element, illustrating that post-depositional actions may be caused by the surface where these elements were resting (Nel, 2013). The row I-sub-squares in M2 Lower (Table 3), located toward the cave wall, had more frequent pitted and etched material than the other sub-squares. It is possible that the occurrences of bacterial attacks in this area caused an etching effect on enamel and bones (Fernández-Jalvo et al., 2002; Fernandez-Jalvo & Andrews, 2016). Thus, the I-sub-squares may represent areas where organic materials were decomposing and creating increased bacterial activities. These are likely the defecation areas used by the animals while the cave was uninhabited during the low-intensity phase of human occupation. There is no apparent relationship regarding post-depositional etching and rounding in M2 Lower and M3, as these damages were not necessarily occurring on the same elements.

CPA, the layer on which the ochre toolkits were found, comprises a thin ochre layer (Haaland, 2018; Henshilwood et al., 2011). There was also a degraded roof spall that had fallen from the cave ceiling into the sub-squares (van Niekerk, personal communication, December 8, 2020, Bergen). The high ratios of etching, pitting, and rounding of micromammal limb bones in CPA is puzzling, particularly as they are all recorded from H6b; and humeri seem to be more affected than femurs and incisors (Tables 3 and 6). The pitting and etching resemble microbial activities, which might have derived from increased use of the cave by animals and humans. However, the rounding of the limb bones could be ascribed to abrasion caused by trampling. Pitting can also result from trampling, but the pitting observed on the bones, combined with etching, is most likely associated with microbial attacks (Fernandez-Jalvo & Andrews, 2016)

CPA comprises an ochre smear interspersed with bone and disintegrating roof spall from the limestone cave roof rich in calcium carbonate. It is this layer on which the ochre toolkits were resting (in G5b), and Haaland (2018) has described CPA as a thin, red lens with spongy bone and shellfish inclusions based on micromorphological analyses (in G7b). Collectively, the taphonomic traits on the micromammal bones and the archaeological evidence point toward a short period of increased use of the cave. The origin of the post-depositional damage is thus likely linked to increased activity, and it is plausible that decomposing organic matter associated with human utilization of the cave was causing the observed damage on the bones.

Taxonomic Distribution Patterns and Palaeoenvironmental Implications

The micromammal taxonomic representation shows variation in both occurrences and relative abundances throughout the M3 and M2 Lower phases (MIS 5c–a). Several species are recorded in the archaeological assemblage but do not seemingly occur in the contemporary Blombos area. However, the ongoing mapping of micromammal species distribution, based on both modern and archaeological materials, highlights the many grey areas of our knowledge regarding the distribution of species habitats in South Africa. In addition, the taxonomic status of several murid species is unclear as current monotypic species may contain other species or sub-species that require further research to determine their habitat requirements (i.e., Otomys spp., Steatomys krebsii) (Matthews & Nel, 2021; Monadjem & Schlitter, 2016,). Yet, while keeping these issues in mind, some inferences of taxonomic distribution patterns, occurrences, and palaeoenvironmental implications can be made based on the BBC assemblage. The presence of Amblysomus hottentotus at BBC, for example, suggests that their range may have reached further west during MIS 5c–a than their currently known distribution, in grassland and woodland habitats in Eastern Cape and KwaZulu-Natal Provinces (Bronner & Mynhardt, 2015; Skinner & Chimimba, 2005). The occurrence of Chlorotalpa duthieae at BBC confirms their presence on the south coast before 97.7 ± 7.6 kya and also a potential westward expansion of their preferred forest habitat with alluvial sands compared to the present.

C. asiatica is present at BBC from 97.7 ± 7.6 kya to 82.5 ± 6.2 kya in MIS 5c/b, and in modern barn owl pellets collected at a farm, about 4 km inland from BBC (Avery et al., 2005; Faith et al., 2019). It is not recorded from MIS 5 coastal sites further to the east, such as Pinnacle Point (Matthews et al., 2009, 2011, 2020) and Klasies River main site (Avery, 1987; Nel et al., 2018). It appears that during MIS 5c/b, there were suitable conditions for the species as far east as BBC, but seemingly not further east along the southern coast. The species disappears in upper M3 (CH/CI) in the BBC record (Jacobs et al., 2020; Nel & Henshilwood, 2016). The continuous presence of G. afra at BBC indicates prevailing alluvial conditions. Thus, a potential scenario for C. asiatica disappearance could have been due to reduced precipitation patterns, making the area less attractive for breeding during winters. However, a better understanding of C. asiatica reproduction patterns and indications of precipitation shifts are necessary to support this hypothesis.

The mesic D. incomptus/capensis indicates that there were wet areas in the vicinity of the cave. The continuous presence and dominance of O. irroratus also point toward sufficient high moisture areas close to the cave throughout MIS 5c–b/a. Thus, it is likely that there was a continuous supply of freshwater, which ensured moist habitats for these mesic species.

Similar to the forest shrew, Myosorex varius, and the southern African vlei rat, Otomys irroratus, Rhabdomys pumilio shares a divergent geographic line that cut across the grassland and fynbos biomes of South Africa (Engelbrecht et al., 2011; Rambau et al., 2003; Willows-Munro & Matthee, 2009, 2011). All these species have multiple lineages/clades, and there was a phylogeographical break at a known vegetation crossover zone around Alice in the Eastern Cape Province. The latter marks the division of grassland biome to the east and fynbos biome to the west (Engelbrecht et al., 2011). However, as Matthews et al. (2020) have also noted, this phylogenetics-based division does not manifest in dental morphology and, therefore, hampers placing the specimens from archaeological assemblages in their respective lineages/clades. The estimated time of divergence into the respective lineages/clades for O. irroratus is ca. 660–620 kya (Engelbrecht et al., 2011). It shows that niche specialization based on biome and rainfall patterns occurred before their appearance in the BBC micromammal assemblage in MIS 5c, before 97.7 ± 7.6 kya (Engelbrecht et al., 2011; Rambau et al., 2003; Willows-Munro & Matthee, 2009, 2011).

The presence of both Acomys subspinosus and Myomyscus verreauxi confirms persistent fynbos vegetation from MIS 5c (prior to 97.7 ± 7.6 kya), with peaks during MIS 5c/b, and a likely decline in fynbos vegetation after 82.5 ± 6.2 kya (Nel & Henshilwood, 2016). Meanwhile, Kreb’s fat mouse occurs at BBC in MIS 5c/b (CH/CI2). This is later than at Pinnacle Point (PP9), where it is present at ca. 130 kya (Matthews et al., 2011, 2020). Steatomys krebsii is absent from other southern Cape sites such as Boomplaas and Die Kelders in the MSA (Avery, 1982; Faith et al., 2019). However, it occurs in the Later Stone Age deposits at these sites. It is abundant in modern pellet collections from Blombos Schoolhouse (4 km from BBC) and Potberg in De Hoop Nature Reserve (Faith et al., 2019), but not in modern pellet assemblages further east on the south coast (Matthews et al., 2020). Monadjem and Schlitter (2016) have noted that the species may represent a complex of several similar species, and this needs further investigation. The occurrences of S. krebsii at both Pinnacle Point and BBC indicate that one or several of the potential sub-species were established on the southern coast by at least MIS 5, but with a seemingly patchy distribution compared to the present.

Matthews et al. (2020) have pointed out that there is a uniform suite of micromammal species (Myosorex varius, Otomys irroratus, Otomys karoensis, Mystromys albicaudatus, Gerbilluscus afra, and Rhabdomys pumilio) prevalent in larger micromammal assemblages from archaeological sites spanning both glacial and interglacial periods on the south coast, and most other locations, in South Africa. The prevalence of these species is also the case for BBC, except M. albicaudatus, which only occurs in one layer. Regardless, this suite of pervasive species occurs despite an extensive and documented palaeoenvironmental change in the region. Thus, Matthews et al. (2020) conclude that these taxa would appear to be of limited use for palaeoenvironmental reconstruction given their demonstrated adaptability and plasticity to diverse environments. However, it should be emphasized that the presence of these species in the assemblages may also reflect a prevalent predator bias. The established predators from most of these archaeological sites are primarily barn owls and spotted eagle owls. They both prefer all the prey mentioned above. A timely question is whether the presence of specialist species, such as Acomys subspinosus, may be underplayed as the generalist species dominate the assemblages and thus blur the nuances of variations in the relative abundance of specialist species. Nevertheless, the stable isotope analysis of these generalist/common species in the archaeological assemblages may yield palaeoenvironmental insights (Leichliter et al., 2016, 2017; Williams et al., 2020). Carbon and oxygen isotopes can reflect vegetation composition and humid/arid events, but this analysis is yet to be understood for most of the micromammal species in South Africa’s archaeological assemblages.

Biodiversity and Habitat Indices

The rarefaction analyses indicate a difference in small mammal community structure which cannot be ascribed to variation in sample size or different predators biasing the assemblage. This notion of change in community structure is supported by the diversity and dominance indices. However, the Simpson dominance index also shows that CK has low dominance, meaning that despite having the least number of individuals in the assemblage (MNI = 53), the layer contains a great diversity of species. Thus, the statistical differences of the individual rarefaction curves compared to CK show that CK is an anomaly with greater species richness than most layers in the sequence. This is confirmed by the diversity index, where only CL and CQ have greater diversity than CK.

On the other end of the scale is CGAC, which has both the least diversity and greatest dominance of one species—Otomys irroratus (Table 7 and Fig. 4). There is a steady decline in diversity from the peak in CL and CK during MIS 5c/b (< 97.7 ± 7.6 kya), culminating in MIS 5b/a (CGAC) (< 82.5 ± 6.2 kya) (Jacobs et al., 2020). In contrast, the dominance of single species correspondingly increases from MIS 5c/b to MIS 5b/a (Fig. 4). Given that micromammal diversity can be used as a measure of vegetation complexity (Cuenca-Bescós et al., 2009), there seems to be a general reduction in the diversity of vegetation, evidence of environmental productivity decline, in the local BBC area between during these two periods. However, the inability to differentiate species, sub-species, and clades is suggestive that biodiversity indices based on taxa are grossly underestimated (Engelbrecht et al., 2011; Matthews et al., 2020).

The THI calculated for BBC does not indicate that any of the major vegetation and substrate habitats vanished during MIS 5c/b–MIS 5b/a (Fig. 6). The reduction in vegetation diversity may not have been so severe that it eradicated any of these major vegetation groups. There are, however, fluctuations in their relative representation, with the increase of moist grass in CK potentially marking a temporary increase in precipitation. At the same time, there is an overall decline in shrubland and possibly reduced heather-like vegetation in M2 Lower (< 82.5 ± 6.2 kya, MIS 5b/a; Table 1; Nel & Henshilwood, 2016). During the M3 phase, the coastline is estimated to be 1.45–2.33 km from BBC, with a potential further regression during M2 Lower to 2.33–9.42 km (Fischer et al., 2010). Bar-Matthews et al. (2010) have suggested that fynbos may have followed the retreating coastline during the latter part of MIS 5. The decline in heather-type vegetation at the BBC during MIS 5b/a could be due to various shrubby fynbos, such as Proteaceae (Table 1), shifting toward the exposed Agulhas Plain.

The precipitation seasonality and aridity in CCA (Fig. 7), M2 Lower, M3 Upper, and M3 Lower are similar to the patterns in the modern pellet locations, 50–150 km west of BBC (Fig. 7). These locations currently receive a greater percentage of cold quarter rainfall, within the range of 30.7 to 39.3%, than in the Blombos area (27.6%) (Supplementary Tables 1 and 2; Faith et al., 2019). Meanwhile, the aridity index is higher (0.52 to 0.54) at Die Kelders and Groot Hagelkraal than at De Hoop Nature Reserve and Blombos (both 0.43) at present, thus suggesting that MIS 5c/b (M3) and MIS 5b/a (M2 Lower) were periods with slightly greater humidity than the present time. This is in line with the conclusion reached by Roberts et al. (2016), who analyzed stable carbon and oxygen isotopes from ostrich eggshell fragments from the M3 and M2 Lower phases. They found that there were humid winter rainfall conditions in early M3. However, there is a discrepancy in our two datasets for the latter part of MIS 5. Their data demonstrate higher δ13C and δ18O, indicating increasing aridity and increased year-round rainfall influence. Our analyses support the increased year-round rainfall influence. However, based on our canonical correspondence analysis (CCA), there are no indications of arid conditions but rather improved humid conditions. MIS 5a (M2 Upper) and MIS 5a/4 (M1) mark the greatest divergence from present rainfall patterns and humid/arid conditions during the BBC’s MSA sequence (Fig. 7). The period after 76.0 ± 5.6 kya show increasing moisture and more aseasonal rainfall than at present and during MIS 5c/b (M3 and M2 Lower), illustrating a shift that coincided with the transition from MIS 5 to MIS 4 (Nel & Henshilwood, 2016).

Conclusion

Our analyses of the micromammal assemblage from 100–77 kya contexts at BBC shed light on site formation processes, micromammal taxonomic distributions, and local palaeoenvironments on the southern Cape coast during MIS 5c–5a. The taphonomic analyses indicate that spotted eagle owls were the main predators responsible for accumulating the micromammals, but with significant contribution from barn owls. These predators may have added a bias to the relative abundance of species in the assemblage. There is a dominance of vlei rats, soricids, and R. pumilio, which are their favored prey. The consequence is that these generalist species are over-represented in abundance at the expense of the less common or specialized species. The assemblage has been subject to several post-depositional processes, although abrasion, pitting, and etching actions are most common. These processes suggest increased activity by occupants (both humans and animals) in certain areas of the cave at particular periods. However, focused taphonomic analyses that are combined with other disciplines such as micromorphology have the potential to reveal further insights into the spatial organization of the cave. The taphonomic signals from our material support the notion that the ochre toolkits were manufactured during an ephemeral visit to the cave.

A few species occurred in BBC archaeological assemblages that are not present in the area today, such as A. hottentotus, C. duthieae, and O. laminatus. As Matthews et al. (2020) and Matthews and Nel (2021) have noted, our knowledge of micromammal distribution patterns for the present and the past has voids and inaccuracies, especially for the biogeography of multiple taxa. These issues need to be addressed with phylogenetic and phylogeography research. Also, as long as we are basing our palaeoenvironmental interpretation on extant species, we need to refine our knowledge of present-day habitat preferences by different species. The biodiversity indices suggest that local vegetation was complex, consisting of different habitats suitable for micromammal species. The vegetation types varied in abundance throughout the sequence, but there are no indications of the disappearance of any particular vegetation component. Grasses dominated, although keeping in mind that the identified predators may overemphasize species associated with grasses such as R. pumilio and Otomyinae. There are indications of drier events during M3 Upper, and on a larger temporal scale, there was seemingly more rainfall during the winter season than at present. The amount of rainfall would have been similar to locations currently 50–150 km further to the west of BBC. Conditions were slightly more humid than the present, but overall there are unexpectedly small variations in the palaeoenvironmental record based on micromammal proxies for the entire 24,000-year period involved in this study. Similar interpretations of only slight changes in vegetation and precipitation during major glacial/interglacial shifts have been reported from micromammal assemblages elsewhere in the region (within 200 km of BBC), on the south coast at Pinnacle Point (Matthews et al., 2020), and inland at Boomplaas Cave (Faith et al., 2019). There are different explanations for these results:

  1. 1

    The reliability of micromammal assemblages as palaeoenvironmental proxies may depend on generalist vs. specialist species composition.

  2. 2

    There is a predator bias clouding the palaeoenvironmental information.

  3. 3

    The local vegetation and precipitation patterns on the south coast do not necessarily reflect global or regional climatic fluctuations on our temporal scale based on micromammal proxies.

    This paper has demonstrated that explanations 1 and 2 are likely for the BBC M2 Lower and M3. The third scenario could also be plausible and reflects the advantage of using micromammals as proxies for palaeoenvironmental reconstruction. The remains of these micromammals provide a geographically constrained snapshot of the existing habitats at the time the owls used the cave for roosting and how these habitats changed over time. These snapshots will likely improve with new methodologies and clarify regional records. A potential key to further refine is the stable isotope analyses of modern and fossil micromammals, and detailed accounts of the biogeography of micromammals and their adaptive mechanisms to climatic fluctuations also hold promise for refining the palaeoenvironmental data derived from micromammal assemblages in the archaeology of South Africa.