ARTICLE https://doi.org/10.1038/s42003-023-04908-w OPEN The earliest Ethiopian wolf: implications for the species evolution and its future survival Bienvenido Martínez-Navarro 1,2,3 ✉, Tegenu Gossa 4,5,6 ✉, Francesco Carotenuto 7, Saverio Bartolini-Lucenti8,9, Paul Palmqvist 10, Asfawossen Asrat11,12, Borja Figueirido 10, Lorenzo Rook 8, Elizabeth M. Niespolo13,14,15, Paul R. Renne 14,15, Gadi Herzlinger5 & Erella Hovers 5,16 ✉ In 2017, a hemimandible (MW5-B208), corresponding to the Ethiopian wolf (Canis simensis), was found in a stratigraphically-controlled and radio-isotopically-dated sequence of the 1234567890():,; Melka Wakena paleoanthropological site-complex, on the Southeastern Ethiopian High- lands, ~ 2300 m above sea level. The specimen is the first and unique Pleistocene fossil of this species. Our data provide an unambiguous minimum age of 1.6–1.4 Ma for the species’ presence in Africa and constitutes the first empirical evidence that supports molecular interpretations. Currently, C. simensis is one of the most endangered carnivore species of Africa. Bioclimate niche modeling applied to the time frame indicated by the fossil suggests that the lineage of the Ethiopian wolf faced severe survival challenges in the past, with consecutive drastic geographic range contractions during warmer periods. These models help to describe future scenarios for the survival of the species. Projections ranging from most pessimistic to most optimistic future climatic scenarios indicate significant reduction of the already-deteriorating territories suitable for the Ethiopian Wolf, increasing the threat to the specie’s future survival. Additionally, the recovery of the Melka Wakena fossil underscores the importance of work outside the East African Rift System in research of early human origins and associated biodiversity on the African continent. 1 ICREA, Pg. Lluís Companys 23, 08010 Barcelona, Spain. 2 Institut Català de Paleoecologia Humana i Evolució Social (IPHES-CERCA), Zona Educacional 4, Campus Sescelades URV (Edifici W3), 43007 Tarragona, Spain. 3 Universitat Rovira i Virgili, Departament d’Història i Història de l’Art, Avinguda de Catalunya 35, 43002 Tarragona, Spain. 4 Human Evolution Research Center (HERC), The University of California at Berkeley, Berkeley, CA, USA. 5 Institute of Archaeology, The Hebrew University of Jerusalem, Jerusalem, Israel. 6 Department of History and Heritage Management, Arba Minch University, Arba Minch, Ethiopia. 7 Department of Earth, Environment and Resource Sciences, University of Naples “Federico II”, Naples, Italy. 8 Earth Science Department, Paleo[Fab]Lab, University of Florence, Via G. La Pira 4, Firenze 50121, Italy. 9 Institut Català de Paleontogia M. Crusafont, Universitat Autònoma de Barcelona, E-08193 Cerdanyola del Vallès, Spain. 10 Departamento de Ecología y Geología, Universidad de Málaga, Universidad de Málaga, Campus de Teatinos, 29071 Málaga, Spain. 11 Department of Mining and Geological Engineering, Botswana International University of Science and Technology, Private Bag 16, Palapye, Botswana. 12 School of Earth Sciences, Addis Ababa University, P. O. Box 1176 Addis Ababa, Ethiopia. 13 Department of Geosciences, Princeton University, Princeton, NJ, USA. 14 Department of Earth and Planetary Science, University of California, Berkeley, CA, USA. 15 Berkeley Geochronology Center, Berkeley, CA, USA. 16 Institute of Human Origins, Arizona State University, Tempe, USA. ✉email: bienvenido.martinez@icrea.cat; tegenu.gossaaredo@mail.huji.ac.il; hovers@mail.huji.ac.il COMMUNICATIONS BIOLOGY | (2023)6:530 | https://doi.org/10.1038/s42003-023-04908-w | www.nature.com/commsbio 1 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-04908-w T he Ethiopian wolf (Canis simensis) is an endemic species of Anatomical and morphometric comparisons with other extant the Ethiopian fauna. Its census population of only ~500 and fossil canids. The dimensions and morphology of MW5- individuals, ~200 of them adults, is spread over six dif- B208 are more similar to the extant C. simensis than to other ferent isolated populations restricted to elevations over 3000 m extant African canids (both Canis and Lupulella species) or to above sea level in the Ethiopian highlands1–8. Given an estimated various fossil canids of Africa (Fig. 3, Supp. Table 1, Supp. minimum viable population size of 3876 individuals, with a 95% Note 4). The height of the corpus, the straight ventral margin, confidence interval of 2261–5095 individuals9, C. simensis is the position of the mental foramina (especially the distal one, currently one of the most endangered carnivores on the African under the diastema between p2-p3), the large diastemata continent. This animal and the large cercopithecid monkey between canine and premolars, and between premolars them- Theropithecus gelada are the most iconic faunal taxa of Ethiopia. selves are similar between MW5-B208 and the extant C. simensis. Although Theropithecus is now only present in the Ethiopian The shallowness and elongation of the mandible shows that highlands, the Pleistocene fossil record of this lineage is well- MW5-B208 falls within the range of variability of the extant known from eastern, southern and northern Africa, as well as in species C. simensis. The similarity extends to dental features: e.g., the Indian Subcontinent, the Levantine Corridor and the Iberian the oval shape of the p4, shortened in mesiodistal direction; the Peninsula10. In contrast, there is no Pleistocene fossil record of reduced distal accessory cuspulid on the p4; the slender m1 with C. simensis from Ethiopia or elsewhere. Thus, the origin of this reduced metaconid and entoconid; the similar-sized m2 proto- species and its time of arrival to Ethiopia remains unknown. conid and metaconid; and the squared distal shape of the m2 Extant C. simensis is more closely linked genetically with the (Figs. 2 and 3). Eurasian and North American gray wolf, Canis lupus, the North Morphometrically, MW5-B208 takes a position on the American coyote, C. latrans, the Asian jackal C. aureus, or the morphospace between the golden jackal and the Ethiopian wolf African small golden wolf C. lupaster, than with any other African (Fig. 4a). But the principal components analysis (PCA) canid such as the African painted dog, Lycaon pictus, or the distribution pattern (cluster) does not help in clarifying the African jackals (Lupulella mesomelas and Lup. adusta) (Supp. morphological affinities of the fossil. As a result, a linear Note 3, Supp. Figure 3). The absence of a fossil record of C. discriminant analysis between the three living species of jackals simensis led to the conclusion that the species arrived in the and the Ethiopian wolf based on six variables (i.e., Lp4, Bp4, Ethiopian highlands during the Late Pleistocene5,11,12. The dis- Lm1trig, Bm1tal, Lcm2, and JBp3p4) was employed (Fig. 4c). covery of a fossil of C. simensis (MW5-B208) in a The analysis reclassifies correctly all the specimens in their stratigraphically-controlled and radio-isotopically-dated context respective groups after cross-validations using the leave-one-out in the paleoanthropological site-complex of Melka Wakena13 method (see Methods subsection “Comparative sample, pro- undermines the latter hypothesis (Figs. 1 and 2). MW5-B208 grams, and software used for the morphological and morpho- represents the first—and (to date) the only—fossil of the Ethio- metric analyses”). This function classifies unequivocally the pian wolf, testifying to its appearance and existence in its Afro- Melka Wakena fossil as a jaw of C. simensis with a probability of montane habitat as early as the second half of the Early p = 0.997. The results of the bootstrapping cluster analysis Pleistocene. (BCA) confirm the general pattern of the PCA, with the The Melka Wakena site-complex is located at ~2300 m above hypercarnivorous canids (C. lupus, C. alpinus, and L. pictus) sea level (asl) (Supp. Figure 1) on the eastern shoulder of the on one side of the plot, whereas the species of the genus Lupulella Main Ethiopian Rift (MER) (Fig. 1a, b). Eight archaeological and (i.e., L. adusta, L. mesomelas), C. aureus, C. latrans, C. simensis two paleontological localities were identified to date, in a and the Canis from Melka Wakena are grouped in the negative chronostratigraphic sequence corresponding to the second half of side of the first principal component. Indeed, the BCA results the Early Pleistocene, ~1.6 to >0.7 Ma13 (Fig. 1c, d; Supp. Notes 1 identify two well-separated clusters composed of the same and 2, Supp. Figure 2). Currently, the vertebrate faunal list species (Fig. 4b). Within the larger cluster, the hemimandible derived from the stratigraphic sequence consists of Crocodylus cf. MW5-B208 from Melka Wakena is grouped with the extant C. simensis. The supporting value of the grouping is the highest niloticus, Hystrix sp., Canis sp., Elephas recki cf. brumpti, Equus of the analysis (percentage p-value = 97%), thus testifying to the sp., Metridichoerus modestus, Hippopotamus gorgops, Giraffa cf. strong similarity in the values of the two species. The results of jumae, Tragelaphus sp., Bovini indet. (Cf. Pelorovis sp.), Aepyceros the BCA fully confirm the classification of the discriminant cf. melampus, Gazella sp., Alcelaphini indet., Reduncini indet., function in clustering the fossil specimen close to the extant and Antelopini indet. This fauna is consistent with a late Early Ethiopian wolf. Pleistocene age. A lower left m3 of the proboscidean Elephas recki cf. brumpti, a Pliocene form, was found in a secondary deposi- tional context, embedded in sediments younger than 1.3 Ma13. Biogeography and phylogenetic relationships with other canids. Recent phylogenetic studies have demonstrated that C. simensis is closely related to species of the crown-group Canis, Results related to the extant gray wolf (C. lupus)14–17 (Supp. Note 3, During the 2017 field season, a right hemimandible (MW5-B208) Supp. Figure 3). This might hint at the ancestry of extant of a canid was found eroding out of a conglomerate layer (Unit II: Ethiopian wolf populations, since true Canis species are well- Fig. 1c, d, Supp. Note 2, Supp. Figure 2) at the base of a cliff-section represented in the Eurasian fossil record from 2.5 Ma and are of MW5-West, ca. 90 meters west of the archaeological locality commonly related to living canids, for example the lineage Canis MW5. This conglomerate layer is found sandwiched between two borjgali–C. orcensis–C. mosbachensis to C. lupus18–21 (Sup. Inf.), volcanic ash layers. The underlying ash flow (Unit-I) is dated by supporting new interpretations17 and contrasting with previous single grain 39Ar/40Ar geochronology to 1626.6 ± 6.4 thousand estimates based on mtDNA of a late Middle–Late Pleistocene years ago (ka) and the overlying ash flow layer (Unit III) dates to differentiation of C. simensis from a wolf-like ancestor at around 1372.5 ± 4.8 ka (Fig. 1d; all errors are reported at 1 sigma, including 100 ka12. An alternative hypothesis based on molecular data— systematic uncertainties)13. The fossil was initially classified as that C. simensis’ time of divergence from other crown-group Canis sp.13. It clearly belongs to a medium-sized canid similar to Canis species was around −2.5 Ma—has thus far lacked the extant endemic Ethiopian Canis simensis (Figs. 2 and 3). empirical evidence17 (Sup. Inf.). The discovery of MW5-B208 at 2 COMMUNICATIONS BIOLOGY | (2023)6:530 | https://doi.org/10.1038/s42003-023-04908-w | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-04908-w ARTICLE Fig. 1 Location and stratigraphy of the site of Melka Wakena (Ethiopia). a, b Relief map from Natural Earth (public domain): http://www. naturalearthdata.com/. c Stratigraphy of the MW5-West sequence. Photo taken by author (T.G.) in 2021. d Stratigraphic position of MW5-B208 in relation to the dated archaeological sequence at locality MW5. Stratigraphic log was created by author (A.A.). the time span of 1.6–1.4 Ma, with all the dentognathic peculia- Modeling the Canis simensis’ niche for past and future pro- rities of its extant relatives and very few differences, constitutes jections. The current regional topography of the Southeastern the first empirical support for the genetic interpretations of Ethiopian Highlands was largely shaped by late Pliocene tectonic C. simensis phylogenetic relationship with the crown-group uplift (ca. 4.5 Ma), after which the landscape stabilized to its Canis. current topographic configuration26. By modeling the species’ The Early Pleistocene dispersal of C. simensis’ ancestors from bioclimatic requirements in the present, we were able to project Eurasia to Africa appears to have been coeval with that of the its ecological niche onto simulated past and future environmental ancestors of the hypercarnivorous Lycaon pictus (African hunting conditions, creating Habitat Suitability Index (HSI) maps that dog), which arrived to Africa later than 1.8 Ma22,23. At the same identify the favorable bioclimatic conditions in Ethiopia through time, 2.0–1.8 Ma, yet in the opposite direction, hominins time for C. simensis (Figs. 5 and 6). Given the species behavior dispersed into Eurasia, as indicated by the record at the site of and ecological needs, this allows us to estimate potential fluc- Dmanisi24,25. In the absence of other earlier evidences, this tuations in its geographic range. The reliability of these projec- hypothesis fits with the new date of the fossil record of tions was guaranteed by checking for models’ transferability C. simensis, 1.6–1.4 Ma. (Supp. Table 2 and Supp. Note 5). COMMUNICATIONS BIOLOGY | (2023)6:530 | https://doi.org/10.1038/s42003-023-04908-w | www.nature.com/commsbio 3 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-04908-w Fig. 2 The fossil from Melka Wakena and the extant Ethiopian wolf. a Canis simensis drawing by S. Bartolini-Lucenti. b the silhouette of C. simensis and the photographs and 3D model of the right mandible MW5-B208 c in buccal, d lingual, e occlusal, and f detailed occlusal view of the p4-m2. Silhouette of Canis simensis made by authors (S.B.L.). Nowadays, C. simensis is adapted to live in the alpine grassland According to our models, changes in the species’ potential environments of the highest Ethiopian mountains, between 3000 geographic distribution were less dynamic up to 1 Ma, whereas it and 4500 m asl4 and with cold climates during most of the year27. contracted and expanded substantially and repeatedly with the The model predictions used in this study (see Bioclimatic onset of 100-ka glacial cycles (Brown-Forsythe Test29: 422.1953, variables in “Methods”) are consistent with the empirical p ≪ 0.001, Fig. 4b). Our tests confirmed a significant and negative information and indicate that the species occupies territories relationship between the species’ potential geographic range and with an annual mean temperature ranging from ~6 °C to 13 °C, both the mean temperatures (ARIMA regression: slope = −0.097, and an annual temperature excursion (i.e., highest temperature of p ≪ 0.001) and altitude sampled by the species (ARIMA the warmest month minus the lowest temperature of the coldest regression: slope = −4.160, p ≪ 0.001) during the last 2 Ma month) ranging from 15 °C to ~22 °C. Also, the annual mean (Supp. Note 5, and Supp. Figures 4 and 5). During periods of precipitation range is wide (from 890 to 1639 mm) with more warmer climate, its potential geographic range was drastically average rainfall occurring during the warmest months (Supp. reduced, forcing the species to move to higher elevations, and Note 5). Our models are consistent with the Early Pleistocene during periods of cooler climate its potential geographic range palynological record of the Melka Kunture archaeological series28, expanded into lower topographic elevations, favoring the also in the Ethiopian highlands (Supp. Note 1, Supp. Figure 1), connection of previously spatially fragmented isolated popula- which indicates that similar favorable environmental conditions tions (Fig. 5c). Interglacial periods, such as the Eemian (from for the survival of the Ethiopian wolf already existed when its about 130 to 115 ka) in our models, represent times of maximum ancestor arrived in Ethiopia (Figs. 5a and 6a, Supp. Note 5). deviation from the estimated favorable climatic conditions for the Within the radio-isotopic age interval for the Melka Wakena species (Fig. 5a) and, in turn, of maximum reduction and site-complex, our results suggest that the time interval with the population fragmentation within its potential geographic range highest percentage of highly suitable territories for C. simensis (Fig. 6d). Thus, these periods of warm climate may have driven was around 1657 ka (Fig. 5a), when favorable conditions occurred the populations close to a bottleneck, and possibly led some of throughout the Ethiopian highlands (Fig. 6b). Conversely, the them to extinction30. However, the general cooling of the climate model suggests a drop in habitat suitability at around 1437 ka during the Middle and Late Pleistocene increased the potential (Fig. 5a), leading to a highly reduced potential geographic range habitat area of C. simensis populations and may have allowed of the Ethiopian wolf (Fig. 5c). them to move to lower elevations, connecting and exchanging 4 COMMUNICATIONS BIOLOGY | (2023)6:530 | https://doi.org/10.1038/s42003-023-04908-w | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-04908-w ARTICLE Fig. 3 Buccal morphology of different Plio-Pleistocene and extant African canid mandibles. The fossil specimens are arranged in three lines in geographic areas (i.e., North, East, and South African sites). Site numeration: (1) Ahl al Oughlam (Morocco); (2) Grotte de Rhinoceros (Morocco); (3) Aralee Issie (Ethiopia); (4) Koobi Fora, Upper Burgi-KBS Members (Kenya); (5) Sterkfontein Member 4 (South Africa); (6) the red star, Melka Wakena (Ethiopia). a right mandible MNHN.F.Th1-10728 of Lupulella mohibi from Grotte de Rhinoceros36. b, c right mandible fragments of Lupulella sp. from Koobi Fora, Upper Burgi-KBS Members37: b KNM-ER 332, and c KNM-ER 895. d Canis brevirostris from Sterkfontein Member 439: 3D model of left hemimandible fragment DNMNH.STS137 (courtesy of Dr. M. Tawane and J. Adams, see Adams et al.40). e, f buccal views of right mandible fragments of Lupulella paralius from Ahl al Oughlam36: e MNHN.F.AaO-4119, and f MNHN.F.AaO-3499. (g, h) left mandible fragments of Eucyon kuta from Aralee Issie38: g left hemimandible fragment MSD-VP-4/438, and h left hemimandible ARI-VP-1/64038. i–k buccal views of mandibles of Lupulella adusta: i right mandible AMNH-216344, j left mandible, reversed, MZUF 8496, and k right mandible AMNH-52049. l–n buccal views of mandibles of Lupulella mesomelas: l right mandible MZUF 2834, m right mandible MZUF 1707, n right mandible MZUF 1843. o–q buccal views of mandibles of Canis lupaster: o right mandible NHM 7434, p left mandible, reversed, MZUF 1851, q right mandible MZUF 2110. r–t buccal views of mandibles of Canis simensis: r right mandible CE 818, s right mandible CE 17801, t right mandible CE 17800. Map and Silhouettes of canids made by authors (S.B.L.). genes with adjacent populations. The overall population size C. simensis, with a succession of several bottlenecks during its would recover when the different populations connected evolutionary history31. physically. In summary, we show that the Ethiopian wolf’s By modeling its favorable climatic conditions, we infer suitable favorable conditions and geographic range as predicted from our territories at present and forecast the fate of C. simensis in the hind-sight models, expanded during cold phases of the Ethiopian highlands over the next 80 years. Our models imply that Pleistocene (Supp. Note 5); the Last Glacial Maximum (LGM present and future global warming constitutes a detrimental factor about 26–18 ka) could have been one of the most favorable phases for the Ethiopian wolf’s survival. In the past 2 Ma, the Ethiopian (Figs. 5b and 6e). These analyses show that C. simensis was always highlands never experienced the temperature increase forecasted an endemic species, with its maximum potential geographic range for the future (Fig. 5a). This increase implies a significant shrinkage restricted to the Ethiopian highlands (Fig. 6b–e), represented by of the Ethiopian wolf’s habitat, which may drastically reduce its small populations, mostly fragmented during the warmer phases geographic range and bring populations to the brink of a of the Quaternary. This is coherent with the hypothesis of the bottleneck. There is also evidence that the current global climate Long-term Small Populations Size and small genetic diversity of change has already induced phenological shifts in top predators, COMMUNICATIONS BIOLOGY | (2023)6:530 | https://doi.org/10.1038/s42003-023-04908-w | www.nature.com/commsbio 5 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-04908-w Fig. 4 Morphometric comparison analyses of the mandible MW5-B208 from Melka Wakena as opposed to extant canids. a Principal component analysis on dentognathic variables. The graph shows a clear separation between the mesocarnivorous canids (on the left) and the three species of pack hunting hypercarnivores (on the right). Yellow star represents the position of MW5-B208 from Melka Wakena. b Bootstrap cluster analyses on PCA scores, confirming the similarity of MW5-B208 to extant Canis simensis. Numbers represent the percentage p-values supporting the node: values b above 95% are reported in bold red. c Linear discriminant function between the three species of jackals and the Ethiopian wolf. This function reclassifies correctly 100% of the specimens of the living canid species and assigns the Melka Wakena fossil jaw to the group of C. simensis. However, the score of the specimen in the discriminant function is closer to the limit between the groups compared (R0) than to the group centroid of Ethiopian wolves (R1), which suggests that the fossil jaw does not show the fully evolved condition of a mandible of C. simensis. Silhouettes of canids made by authors (S.B.L.). with drastic consequences for their reproductive success32. Current physical connections (MIROC6 model, GLM regression, Supp. climatic conditions already determine very restricted contiguous Note 5 and Supp. Table 5), although different future climatic areas suitable to the species (Fig. 6f and Supp. Table 3). For the projections may impact differently on the species’ habitat (Supp. future, even optimistic projections, taking into account adequate Note 5 and Supp. Table 6). climate policies (Models: CNRM-CM6-1 and CNRM-ESM2-1, ssp: 126, with mean global temperature increase lower than 2 °C in 2100, Supp. Note 5), indicate that the most favorable habitats will Discussion be reduced up to about 65% of its current extent by 2080 (Supp. Our comparative analysis of the extant specimens of African Note 5) and that they will totally disappear afterward (Supp. canids shows that MW5-B208 belongs to a fossil form of the Note 5, Supp. Tab. 4). Indeed, more pessimistic projections of Canis simensis lineage (Figs. 3 and 4; Supp. Tab. 1, Supp. Note 4). faster increases of temperature and CO2 emissions predict a habitat Both anatomical and morphometric analyses confirmed similarity loss of ~95% (Model: BCC-CSM2-MR, ssp: 585, year: 2100, Supp. between the specimen from Melka Wakena and the extant Note 5) and up to 100% (Model: CNRM-CM6-1, ssp: 370 and 585, Ethiopian wolf. It is worth noting that the fossil takes a dis- year: 2080; Model: CNRM-ESM2-1, ssp: 585, year: 2100) of the criminant score that situates it in the lower limit of the range of most suitable territories (see Supp. Note 5, Supp. Table 4, for values shown by Ethiopian wolves in the discriminant function. intermediate scenarios). By considering an intermediate- This may be tentatively interpreted as indicating that while pessimistic increase of temperature (BCC-CSM2-MR, ssp: 370), MW5-B208 shows the metric features of the lower dentition and we forecasted a reduction of the species’ most suitable conditions of jaw of a specimen of C. simensis, it does not represent the fully about 28% by 2040 (Fig. 6g), ~63% by 2060 (Fig. 6h), ~74% by evolved condition of this species. Given that the Early Pleistocene 2080 (Fig. 6i), and ~87% by 2100 (Fig. 6j). In general, the preferred altitude of Melka Wakena was similar to the current situation conditions for the Ethiopian wolf will be confined to the mountain (~2300 m above sea level)26, the specimen MW5-B208 constitutes tops but will be highly reduced geographically (OLS regression, the first-ever evidence that C. simensis was present on the Supp. Note 5) and fragmented into many patches with poor Ethiopian highlands from at least 1.6–1.4 Ma. 6 COMMUNICATIONS BIOLOGY | (2023)6:530 | https://doi.org/10.1038/s42003-023-04908-w | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-04908-w ARTICLE Fig. 5 Temporal trend of ecological and environmental parameters of C. simensis in Ethiopia and Melka Wakena. a Temporal trend of Habitat Suitability Index (HSI) values computed for the last 2 Ma at the site of Melka Wakena (violet line), as mean value in the whole Ethiopia (blue line). At the base of this panel (a) a sequence of green points indicates the temporal intervals with environmental conditions in Ethiopia better than the present time interval; below, the sequence of orange points indicates temporal intervals with conditions similar to the current climate in Ethiopia. In some cases, these two different sequences of green and orange points apparently overlap in time, thus generating time bins that seem to include both better and similar conditions to present Ethiopia. This is only a graphical issue due to limit of minimum pixel size in the image. In order to make present and past time bins comparable, current climatic conditions were determined by using the Raia et al.55 raster layer for the present. b Temporal trend of the species’ geographic range’s Log10 surface area. For species’ geographic range computation, see “Methods” section. c Temporal trend of the Log10 mean altitude sampled by the species in the reconstructed geographic range. Orange vertical dashed lines indicate time intervals important to the fossil localitity of Melka Wakena: 1690 ka, the maximum estimated radiometric age of the site; 1657 ka, age with the highest HSI value reconstructed for the site in the estimated age range; 1437 ka, age with the lowest HSI value reconstructed for the site in the estimated age range; 1324 ka, the minimum estimated radiometric age of the site; LIG, Last Interglacial (~123 ka); LGM, Last Glacial Maximum (~22 ka). The geographic position of Melka Wakena coincides with the likelihood was high that the species’ fossilized remains would territory where our model identifies some of the highest values of have been found in this area34. the HSI (i.e., indicating highly favorable conditions for the spe- Our study renders the absence (until the current discovery) of cies) throughout Ethiopia over the last 2 Ma (these values were C. simensis from the fossil record a likely result of the scarcity of higher than 99% of the area of the whole region). Since habitat documented archaeological and paleontological sites on the suitability is typically correlated with population abundance33, the Ethiopian highlands. The bulk of the Quaternary fossil record in COMMUNICATIONS BIOLOGY | (2023)6:530 | https://doi.org/10.1038/s42003-023-04908-w | www.nature.com/commsbio 7 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-04908-w Fig. 6 Temporal trend of climate in Ethiopia and of habitat suitable to C. simensis. a Reconstructed Annual Mean Temperature (BIO 1) averaged for the whole Ethiopia territory during the last 2 Ma, based on Raia et al.55, the present and future forecasts in 2040, 2060, 2080, and 2100 (based on the model BCC-CSM2-MR with ssp: 370, see Supp. Note 5 for further information). Before averaging temperature values, present and future temperature raster were upscaled to a resolution of 50 × 50 km, as for the past bioclimatic rasters. Blue dots indicate the temperature values for the present, for 2040, 2060, 2080, and 2100 (see Supp. Note 5 for other future models and scenarios). The Last Interglacial (LIG) and the Last Glacial Maximum (LGM) are also indicated by the dashed-gray vertical lines at 123 and 22 ka; b–e maps of the reconstructed Ethiopian wolf’s Habitat Suitability Index in Ethiopia for 1657 Ka (b), 1437 Ka (c), 123 ka (LIG, d), and for 22 ka (LGM, e); and f–j Habitat Suitability Index map of the present (f) and forecasted maps for the 2040 (g), 2060 (h), 2080 (i), and 2100 (j) according to the CMIP6 model BCC-CSM2-MR with the shared socio-economic pathway 370 (see “Methods” for further information about CMIP6 models). Map of Africa modified by the authors (F.C.) from https://neo.gsfc.nasa.gov/view.php?datasetId=BlueMarbleNG, credits to NEO (Nasa Earth Observations, https://neo.gsfc.nasa.gov/). Maps showing the C. simensis spatial distribution through time generated modified by authors (F.C.) by applying Maxent models to the bioclimatic variables’ rasters. eastern Africa derives from the rich paleoanthropological and since their arrival to the Ethiopian highlands. Additional informa- paleontological sites in the rift valley, where C. simensis most tion, with the application of new methods to reconstruct the geo- probably would not be present, given its adaptation for highland graphic distributions of past species35, is needed to understand how environments. The discovery of MW5-B208 thus underscores the C. simensis overcame past periods of habitat loss, which might help importance of paleoanthropological research outside the East in building conservation programs for this iconic animal. African Rift System. The Ethiopian wolf is one of the most iconic species on the African continent. While our study is not aimed to outline con- Methods Comparative sample, programs, and software used for the morphological and servation strategies, it does underscore the urgent need for such morphometric analyses. The present study describes and analyses a canid programs in order to safeguard the remaining habitats of this mandible (MW5-B208) recovered from the Early Pleistocene site Melka Wakena ecologically specialized, highly endangered species and to prevent (Ethiopia), comparing it taxonomically to other African Plio-Pleistocene medium- its reproductive depression. Since its adaptation to the Afroalpine sized canids and to selected extant species of the genera Lupulella and Canis. The subject fossil is housed at the Ethiopian Heritage Authority (EHA) in Addis Ababa. environments, at least 1.6–1.4 Ma ago, C. simensis most probably The Villafranchian and Epivillafranchian canids from Africa housed at Musée has been at the brink of extinction several times due to climate National d’Histoire Naturelle (MNHN, Paris) as well others from the literature changes, yet it has always recovered. This species is a true bio- were utilized as comparative fossil material. This fossil comparative sample logical legacy that belongs to the world’s nature heritage. It is an includes specimens of Lupulella paralius Geraads, 2011 from Ahl al Oughlam, ethical as well as a scientific duty to preserve and draw the Morocco36; Lupulella mohibi Geraads, 2011 from Grotte de Rhinoceros, Morocco36; Lupulella sp. from Koobi Fora, Upper Burgi-KBS Members37; Eucyon Ethiopian wolf from the brink of extinction and turn it into a kuta Werdelin et al., 2015 from Aralee Issie, Ethiopia38 and Canis brevirostris Ewer, persistent animal. 1956 from Sterkfontein Member 4, South Africa39 (digital data obtained thanks to The discovery of the ~1.6–1.4 Ma specimen MW5-B208, coupled the courtesy of Dr. M. Tawane and Dr. J. Adams from Morphosource40). Extant with the genetic data that suggest its speciation at an even earlier specimens housed at the American Museum of Natural History (AMNH, New York), “La Specola” Museo di Zoologia di Firenze (MZUF, Florence), Nat- date, suggests that populations of the Ethiopian wolf likely survived urhistorisches Museum Basel (NMB, Basel), Musée National d’Histoire Naturelle unfavorable environments throughout a number of climatic cycles (MNHN, Paris), and EHA (Addis Ababa) were also used for morphological and 8 COMMUNICATIONS BIOLOGY | (2023)6:530 | https://doi.org/10.1038/s42003-023-04908-w | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-04908-w ARTICLE metrical comparisons. We examined specimens of Lupulella adusta (Sundevall, at the posterior border of the canine; JDc: the height of the mandible at the 1847), Lupulella mesomelas (Schreber, 1775), Canis lupaster Hemprich and posterior border of the canine; JBp3p4: breadth of the mandible measured at the Ehrenberg, 1832, and Canis simensis Rüppell, 1840. Moreover, a wide dataset of interalveolar space between the third and fourth premolars; JDp3p4: the height of craniodental measurements of modern canids, taken by Prof. Blaire Van Valk- the mandible measured at the interalveolar space between third and fourth enburgh, was used also in some statistical comparisons, including a discriminant premolars; JBp4m1: breadth of the mandible measured at the interalveolar space analysis between omnivorous (i.e., meso- and hypo-carnivorous) and hypercarni- between the fourth premolar and the carnassial; JDp4m1: the height of the vorous canids, in order to reach paleoecological inferences for the Melka Wakena mandible measured at the interalveolar space between the fourth premolar and the canid. Cranial and dental measurements were taken with a digital caliper to the carnassial; JBm1m2: breadth of the mandible measured at the interalveolar space nearest 0.1 mm41. Photographs of the specimens were taken with a camera Canon between the carnassial and the second molar; JDm1m2: the height of the mandible EOS 200D and with a Canon 60D and subsequently processed with Photoshop CC measured at the interalveolar space between the carnassial and the second molar; L: 2017 (release 2017.0.1). Inkscape ver 0.92.4 was used to compose and finalize the maximum mesiodistal length of the tooth measured at the base of the tooth crown; figures. Three-dimensional data of specimen MW5-B208 were acquired using a Lcm2: the length between the anterior tip of the canine and the posterior border of handheld high-resolution 3D scanner EinScan Pro HD (www.einscan.com) located the second molar; m1tal: trigonid portion of the lower carnassial; m1trig: trigonid at EHA. Scanning was carried out by one of the authors (T.G.) with the help of portion of the lower carnassial. EHA curator Sahlesellasie Melaku. This hand-held scanner is compact and highly suitable for in situ scans of small to medium-sized objects. Its field-of-view (FOV) ranges between a maximum of FOV = 31.0 × 24.0 cm and a minimum of Modeling the bioclimatic niche through time. The reconstructions of past and 20.9 × 16.0 cm with an adjustable scan speed between 10 to 30 frame per seconds. future spatial distribution of the Ethiopian wolf and related analyses require Its accuracy is around 0.1–0.3 mm depending on the chosen preset (respectively mathematical modeling of the species’ bioclimatic niche followed by the projection “scan” or “rapid scan”), making it suitable for small objects. In order to obtain the of this niche onto the geographic space. This procedure is known as Species Dis- 3D model of the specimen, raw data were acquired with multiple passes of the tribution Modeling (SDM)50 and is widely employed for both living51 and fossil scanner on the mandible, which was oriented in different views and positions to species35,52. The output of this projection is a map of the species’ Habitat Suitability allow a full coverage of the surface of the specimen. Acquisition and subsequent Index (HSI map), which is the reconstruction of the territories suitable to the elaboration of the data were performed in the native software EinScan-Pro series species’ habitat in the geographic space and that can be used to derive the species’ ver. 3.1.0.2. This program enabled cleaning, alignment, and refinement of the data geographic range and metrics of the landscape connectivity structure. The proce- prior to the creation of the 3D mesh of the specimen. The rapidity and quality of dures used to model the bioclimatic niche of C. simensis, to forecast its geographic the scanner allowed obtaining a valuable digital reproduction of the fossil in a range expansion or contraction in the next 80 years and to project these conditions short time. to the past 2 Ma are described in the following section. We present future pro- Morphometric analyses42–47 of the metric affinities of the Melka Wakena jections for C. simensis’ under different scenarios of climate change. Using models specimen were carried out using nineteen metric variables, measured on the teeth of past temporal intervals, we explored whether the habitat of Ethiopian wolf and on the portions of the mandibular ramus preserved in the fossil jaw. The latter already existed before and during the time span suggested for human occupations include: the anteroposterior length and buccolingual breadth of the lower canine at the Melka Wakena site-complex and to study the species’ temporal evolution in (Lc and Bc, respectively) measured at the base of the tooth crown, the length and the area of present-day Ethiopia; for the future, we aimed at estimating the degree breadth of the fourth premolar (Lp4 and Bp4, respectively), the length of the of habitat loss caused by global warming. Furthermore, the past and future habitat trigonid blade in the lower carnassial (Lm1trig), the maximum breadth of the reconstructions were built upon different environmental variables and at different carnassial (Bm1), the length and breadth of the talonid basin in the carnassial spatio-temporal resolutions. This means that past and future C. simensis’ habitat (Lm1tal and Bm1tal, respectively), the length and breadth of the second molar reconstructions cannot be compared to each other although being accurately (Lm2 and Bm2, respectively), the length between the anterior tip of the canine and reconstructed for their scopes. We named the models for projecting the bioclimatic the posterior border of the second molar (Lcm2), the height and breadth of the niche onto past temporal intervals the “Hindcasting Model” (HM), whereas the mandible at the posterior border of the canine (JDc and JBc, respectively), the models for projections onto the future the “Forecasting Model” (FM). height and breadth of the mandible at the contact between the third and fourth premolars (JDp3p4 and JBp3p4, respectively), the height and breadth of the mandible at the contact between the fourth premolar and the carnassial (JDp4m1 Species’ occurrences. Information about the present occurrences of C. simensis is and JBp4m1, respectively), and the height and breadth of the mandible at the unsystematic and sporadic. GBIF public repository (www.gbif.org) registered 96 contact between the carnassial and the second molar (JDm1m2 and JBm1m2, occurrences, but the error of the geographical positions was too high to allow an respectively). Abbreviations are also reported in a separate section. accurate spatial modeling. Instead, we used the geographic range polygons pro- A principal components analysis (PCA) was performed after logarithmic vided by the IUCN (International Union for the Conservation of Nature, https:// transformation of these metric variables on a sample of seven species of living canids, www.iucnredlist.org/) for C. simensis as a framework for simulating species which includes four mesocarnivorous species [the Ethiopian wolf (Canis simensis, occurrences. According to the most recent available estimates, there are 197 mature n = 20), the golden jackal (C. aureus, n = 19), the side-striped jackal (Lupulella adusta, reproductive individuals within the total census population4; therefore, we simu- n = 21), the black-backed jackal (Lup. mesomelas, n = 21) and the coyote (C. latrans, lated a total of 200 random points falling inside the polygons of the species’ n = 21)]. The sample also includes three hypercarnivorous, pack hunting species [the geographic range. This procedure is widely used for modeling species distribution gray wolf (C. lupus, n = 22), the African painted dog (Lycaon pictus, n = 19), and the in the absence of accurate true presence data53. Because the species is narrowly dhole (Cuon alpinus, n = 21)]. Two other tests were performed to assess the similarity distributed and the IUCN polygons are drawn from empirical data and the species of MW5-B208 to the other fossil species in our sample. Firstly, we performed a linear actual occurrences and movements, simulating presence points inside the polygons discriminant analysis between the three species of jackals (n = 61) and the Ethiopian approximates the populations’ general distribution. In our simulation, only one wolf (n = 20), using the Wilks lambda criterium for selection of variables (λwilks = presence point could occur in any single grid cell of the bioclimatic variables raster 0.132, c2 = 153.912, p < 0.0001). The linear function obtained incorporates six variables used to model the specie’s niche. This procedure is helpful to reduce overfitting (Φ = −0.577Lp4 + 1.545Bp4 + 0.277Lm1trig − 0.839Bm1tal + 0.256Lcm2 − 0.623 during the modeling process. Since we used bioclimatic variables with different JBp3p4 − 14.966) attempting to reclassify each specimen via cross-validations with the spatial resolutions for the Hindcasting and Forecasting models we used (see below), leave-one-out method. A second classification algorithm performed on the resulting we produced two sets of simulated occurrences, one each for the past and for the data of the PCA is a permutation cluster analysis: bootstrapping cluster analysis future models, with spatial grids of 4 × 4 and 1 × 1 km cell resolution, respectively. (hereafter BCA) is a partitioning method, helpful to test the robustness of the clustering results. We briefly outline the way the methodology works: a preliminary clustering of observations is performed via an UPGMA algorithm (unweighted pair Bioclimatic variables. We used different sets of bioclimatic variables to first model group method with arithmetic mean)48 to be used as a reference cluster in which the the species’ bioclimatic niche on, and then to project it onto future and past time fossil species are grouped according to similarity in their values. Subsequently, the intervals. The Forecasting models and the Hindcasting models derive from the algorithm assigns a probability value to each partitioning level in the reference same set of current bioclimatic conditions, differing in the number of the variables clustering. Then a sample of the original data (i.e., rotated a) is randomly chosen and considered in the models. Following this, the Forecasting and Hindcasting models used to feed a new UPGMA cluster analysis, which yields a new cluster. The were projected onto different temporal intervals, in the next 80 years and within comparison between the reference and the sample clusters yields a similarity index (G) the past 2 Ma, respectively. that ranges between 0 (if the two clusters are totally different) and 1 (if the clusters For calibrating the Hindcasting models, we used four of the nineteen generated with the original and sampled data coincides). The G* is then compared bioclimatic variables of the present time interval as provided by www.worldclim. with the expected similarity value between reference and sample-based clustering (G°) 54 org , with a spatial resolution of 4 × 4 km grid cell combined with a Net Primary under the null hypothesis that the sampled dataset is a truly random sample of the Productivity (NPP) layer (see the next paragraph for details). For past projections original data. The resampling the original dataset was reiterated 1000 times. In the end, of the species’ niche models, we used the improved version55 of the PLASIM- if the probability that G* is higher or equal to G° is higher than the significant level GENIE paleoclimatic emulator data56 including four bioclimatic and NPP layers at (P(G° ≤ G*); α = 0.05), the partitioning levels of the reference cluster analysis are sharp a spatial resolution of 50 × 50 km grid cell and a temporal step of 1000 years. The (for a detailed explanation of the algorithm, see DePatta Pillar49). PLASIM-GENIE bioclimatic dataset spans the last 5 Ma, including the Abbreviations for the morphometric analyses: B: maximum buccolingual breadth reconstruction of the current climatic conditions at the same spatial resolution, a of the tooth measured at the base of the tooth crown; JBc: breadth of the mandible fact that makes past habitat reconstructions comparable with present habitat COMMUNICATIONS BIOLOGY | (2023)6:530 | https://doi.org/10.1038/s42003-023-04908-w | www.nature.com/commsbio 9 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-04908-w estimations (see below). We used these climate projections to reconstruct the Hindcasting models, we only used five bioclimatic variables of past intervals. While species’ spatial distribution from the present to 2 Ma. employing PC axes would have reduced variable’ collinearity, a small number of For calibrating the Forecasting models for projection onto the future, we used predictors is not advisable because it could determine model’s overfitting. Hence, all the nineteen bioclimatic variables at a spatial resolution of 1 km side grid cell, as the strategy of Breiner et al.62 was prioritized. We trained five small models, each provided at http://www.worldclim.org, in order to reconstruct the most detailed built upon a combination of four out of the five available PC variables. For fore- current realized species’ niche. We used the bioclimatic variables provided by casting, we built 28 models, each using a combination of six out of the eight http://www.worldclim.org at a 4 × 4 km grid cell resolution for the four future time bioclimatic variables. intervals: 2021 to 2040, from 2041 to 2060, from 2061 to 2080, and from 2081 to 2100. Future climates are simulated by different research institutions that are The Maxent model. The SDM we used to predict the species’ habitat suitability coordinated by the Coupled Model Intercomparison Projects, now at its stage 6 was the Maxent algorithm63. Maxent belongs the family of presence-only methods (CMIP6). Each institution uses a specific emulator to predict future climate. in which only species occurrences are required other than the environmental Nevertheless, all these future projections are modeled by taking into account variables. Being currently one of the most popular algorithms for presence-only different shared socio-economic pathways scenarios (ssp). These scenarios datasets, Maxent relies on the assumption that the observed species occurrence hypothesize future socio-economic projections and political environments which locations represent a random sample of its true geographic distribution. The aim of would drive to different CO2 emissions’ rates in future years57. The CO2 emission the algorithm is to find the true species probability distribution in the study area, rate is a parameter used to forecast the future temperature variation. Each which is initially unknown. Further, since many theoretical assumptions are vio- institution’s emulator yields a climate forecast that has a specific response to the lated in the real world, the Maxent output is better interpreted as measure of the above-mentioned CO2 emission scenarios. This means that each emulator species habitat suitability in the study area than its occurrence probability (for a produces a climate forecast by considering a specific “climate sensitivity”, i.e., the broad discussion about the Maxent output interpretation, see Merow et al.64 and rate of temperature increase with the increase of CO2 in the atmosphere. The Elith et al.65). In Maxent algorithm the true (unknown) species’ occurrence CMIP6 project produced future climates with a sensitivity ranging from 1.8 and probability in the study area (the landscape, sensu Elith et al.65 is a function of the 5.6. The IPCC AR5 report58 considered as plausible a range from 1.5 and 4.5. For Probability Density Function (PDF) of the environmental variables sampled at the our future projections of the Ethiopian wolf climatic niche, we consider four out of presence locations and of a prior hypothesis. Maxent aims to find an approx- all the available CMIP6 emulators: the MIROC6 with a sensitivity of 2.6; the BCC- imation of this unknown true probability distribution. The theoretical framework CSM2-MR emulator with a climate sensitivity of 3; CNRM-CM6-1 with a of this algorithm is given by the maximum entropy principle66, according to which sensitivity of 4.3; and CNRM-ESM2-1 with a more extreme sensitivity of 4.8. For the best approach to approximate the true PDF of an event is to find the dis- each CMIP6 model, we considered the four most important ssps with an increasing tribution with the maximum entropy, while still satisfying any known constraint63. CO2 emission rate through time (for further details, Meinshausen et al.57): ssp 126, In other words, to assess the possible state of an event at an unknown location, we ssp 245, ssp 370, and ssp 585. Present, past and future bioclimatic variable rasters have to maximize the number of possible choices of the state of the event. Max- were projected into a Lambert Azimuthal Equal Area projection with projection imizing the number of choices means maximizing the entropy of the system67. The centre at the African continent before performing Species Distribution Models Maxent algorithm starts from the most maximized and unconditioned PDF, which (SDM). We used all the bioclimatic data as provided at the downloading date, is represented by the uniform distribution63,64, and applies some constraints 13 June 2021. determined by the characteristics of the PDFs of the variables at presence points. In thermodynamic terms, finding the PDF with the maximum entropy equals to Variables’ pre-processing. The model used to project the Ethiopian wolf climatic expanding a gas into a larger volume after providing some energy to its molecules niche onto past time intervals could not use all the present bioclimatic variables (thus increasing the molecules’ chaotic movements) but controlling for the since the raster dataset54,55 includes only four out of the nineteen bioclimatic volume’s shape. With a uniform distribution a species has the same chance to be variables available for the present time interval, with an additional NPP layer. We present in any location of the study area. In ecological terms, this means that the included in the Hindcasting Models the following variables: BIO 10 (Mean Tem- probability of the species to be present in a location is a function of the relative perature of Warmest Quarter), BIO 11 (Mean Temperature of Coldest Quarter), abundance in the landscape of the environmental conditions found there64. Hence, BIO 18 (Precipitation of Warmest Quarter), BIO 19 (Precipitation of Coldest the probability density function of environmental variable values in the study area Quarter) and NPP. To reduce the effect of multicollinearity between these vari- represents our prior hypothesis about the distribution with maximum entropy. If Z ables, we performed a Principal Component Analysis and considered the related is the set of predictor variables geographically distributed over the study, this prior five PC scores for all the modeling processes. Santini et al.59 recommended to check is defined as Q(z) and the Maxent algorithm can be explained by calculating the for exportability through time of the models, since models calibrated with present probability of finding the variables’ values favouring the species’ presence. This variables might produce extrapolations when projecting onto climatic conditions of turns that Maxent algorithm has to find the distribution of predictor values that is different time intervals. We performed a Multivariate Environmental Similarity the closest to Q(z). The probability density function of these variables at the Surface analysis (MESS)60 on all the variables of all the time intervals to detect locations where the species is observed is P(z), whereas P*(z) is the distribution we those raster’s pixels with values’ combinations not included in the present climate. want Maxent calculates by forcing it to be as closed as possible to Q(z). If we Before performing the MESS analyses, we first turned past variables into PCA remember the definition of maximum entropy described above, Maxent has to find scores by applying on them the PCA model calibrated with the variables of the a P*(z) that also satisfies some constrains. These constrains are determined by the present time interval. After the MESS analyses, we ruled out a specific past bio- statical characteristics of the environmental values at the presence points. Maxent climatic variable if, considering all the temporal intervals, its mean proportion of applies these constraints to the raw environmental variables and turns them into the out-of-the-range values in the total Africa raster would have exceeded 20%. We “features”, i.e., functions of original predictors. The details about the features are anticipate here that none of the five PCA raster variables would be ruled out from explained in Phillips et al.63 and Merow et al.64. the modeling procedures. The mathematical formulation of Maxent model is the following, from Merow For future projections of the climate, we performed the same MESS analyses et al.64, Equation 1. but, because future climate projections shared the same bioclimatic variables with Qðxi Þ expðzðxi ÞIÞ current conditions, we did not perform the PCA. Instead, we reduced the possible P* ðzðxi ÞÞ ¼ ð1Þ effects of multicollinearity by removing from the set of present bioclimatic variables ∑ i Qðxi Þ expðzðxi ÞIÞ those with a Variance Inflation Factor higher than 361. After this procedure, eight where z(xi) is the vector of predictor/environmental values at location xi and I is out of nineteen bioclimatic variables proved useful for the modeling process. For the vector of regressor coefficients. The aim of Maxent is to maximize the similarity the future projections, too, none of the remaining bioclimatic variables showed between the first term of the equation (the distribution of environmental vales critical issues of extrapolation after the MESS analyses. favouring the species’ presence/the prediction) and the second term, i.e., the We performed the described analyses and pre-processing procedures with all constrained prior distribution. The prior distribution of the environmental the CMIP6 future projections and ssp scenarios. variables is taken from a sample of the entire study area called background points. The way we choose these points can strongly affect the Maxent predictions64. Since the second term of the equation has as denominator the PDF of variables in the The hindcasting and forecasting modeling procedures. As mentioned above, whole study area, the sum of all the P*(z(xi)) is 164. This means that the final scores Hindcasting Models (HMs) are trained for past niche projections, whereas the yielded by the algorithm are very small and difficult to interpret. Maxent turns Forecasting Models (FMs) were employed for projecting species’ niche onto future these values into a logistic output that can be interpreted as Habitat Suitability ones. Both sets of models were calibrated by combining the simulated species’ Index scores64,65. occurrences, as described above, with the climate of the present. For the calibration of both Hindcasting and Forecasting models we employed the procedure called “ensemble of small models”62. Rather than creating a single model with all the Model cross-validation and parameter’ tuning. Classical calibration-testing available variables, this procedure combines a number of models built upon dif- protocols require the splitting of the occurrence dataset into training and testing ferent combinations of a subset of the predictor variables. The authors suggest to data, the former used to build the mathematical model and the latter to measure average the predictions of all the models by weighting them according to their the model’s prediction accuracy. This procedure is known as cross-validation and is accuracy. This procedure is a strategy to model the bioclimatic niche of small- particularly important when such models have to be transferred to different ter- ranged, endemic species, like C. simensis, in order to avoid model overfitting, an ritories for conservation strategies or for reconstructing the temporal evolution of a issue that may result in an underestimation of the species’ niche breadth. For the species’ habitat, as in the current study. The calibrated model is then applied to 10 COMMUNICATIONS BIOLOGY | (2023)6:530 | https://doi.org/10.1038/s42003-023-04908-w | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-04908-w ARTICLE reconstruct the known testing data and the related predictive performance is threshold value as a species geographic range. There are many strategies to find the measured by specific indices, like AUC, Boyce, and TSS. The cross-validation is an most appropriate threshold value71, such as using the suitability value maximizing iterative process, and its algorithm can be built according to the purposes of the the difference between correctly identified true presences and true absences. research. In machine learning protocols, the cross-validation procedure has the Threshold values are computed by using the data for training the models and then important role of choosing the best combination of model parameters, thus applied to the habitat suitability maps. In our very case, since the bioclimatic allowing control for the model’s accuracy and complexity, whose balance has a variables of the present that were used for training the model and those of the past prominent role for the model’s transferability59,68. The Maxent algorithm needs the time intervals have different spatial resolutions, computing the threshold value by tuning of different parameters, such as the kind of features to combine and to apply using the training data is not a valid strategy. Coarse-grained habitat suitability to the predictor environmental variables (i.e., the constraints discussed above)65 maps of the past might not include all the HSI values computed for the present and a regularization multiplier, the function of which is to avoid model overfitting. temporal interval. Indeed, when we used threshold values computed with training We applied the Maxent algorithm of Phillips et al.63, as provided in the package data, many past temporal intervals yielded no (i.e., zero) geographic range of the “ENMeval” ver 2.0.369 of the R software environment42. This package allows to set species. different cross-validation strategies for the tuning of all the model parameters Here we used a special iterative strategy to compute a threshold value for (features and the regularization multiplier). We followed the strategy suggested by successfully binarizing coarse-grained HSI maps computed with the bioclimatic Helmstetter et al.68 and Santini et al.59, which resulted in the best way to tune a rasters55. This strategy, which we dubbed Minimum Useful Threshold (MUT), Maxent parameters in order to ensure transferability of the model through space prevented us from computing a null species geographic range for the past time and time. This procedure consists of dividing the occurrence data into different, intervals: (1) we created a vector of different thresholds ranging from 0.01 to 1, spatially separated folds, using subsets of the data to calibrate the model. The with a step of 0.01; (2) we applied each one of these threshold values to binarize the remaining subsets of data were used to validate the model’s predictive performance suitability maps computed from 0 to 2 Ma; (3) we counted the number of raster cell by means of the AUC metric. The algorithm builds many Maxent models itera- selected as species geographic range for each time bin. Then we chose as a valid tively by testing different combinations of parameters and then provides the related threshold the lowest suitability value higher than zero geographic range in the performance metric. In the end, the user can choose the best performance model, complete considered time bin. We then applied this threshold to the continuous i.e., that which most correctly predicts the occurrences in the testing dataset. This is habitat suitability maps of all the time bins, to derive the species’ geographic range a well-known “leave one out” strategy. Specifically, we divided the simulated and to provide statistics about its temporal evolution. After computing the species’ occurrence dataset into twelve spatially separated groups, roughly coinciding with geographic range, we tested if in the past time bins there were climatic conditions the polygons of the species’ spatial distribution provided by https://www. as suitable—or more beneficial—to the species as those in the current territory of iucnredlist.org/ (the same used to simulate presences, as described above). Since Ethiopia. To this aim, we computed for each of the past 0.01 Ma-time interval a these polygons have a north-south geographic distribution, we are sure that this Wilcoxon test between the habitat suitability values sampled by the species in its spatial block strategy is the one that can ensure temporal model’s transferability. estimated geographic range at the specific time bin and the same values sampled in This strategy was used for both Hindcasting and Forecasting models. the current geographic distribution of the species. In addition, we counted the time intervals for which we reconstructed statistically significant habitat suitability values similar to, or higher than, the computed ones for the present in the territory Background record simulation. Maxent algorithm is a presence-only method and of Ethiopia. Since the amplitude of glacial cycles switched from a ~41 ka to ~100 ka only requires species’ occurrences and environmental data to work. As explained in at about 1 Ma, we hypothesized that the variance of the species’ geographic ranges Eq. 1, the species occurrence PDF needs to be regularized by dividing it by the PDF computed for the time intervals younger than 1 Ma was higher than the variance of environmental variables of the Landscape. Maxent uses a sample of the land- during intervals older than 1 Ma. To test this hypothesis, we first averaged HSI scape represented by what are known as background point data. These points can values within the species’ geographic range in each time bin and then performed a be simulated over the study area following different strategies and should be Brown-Forsythe test for equality of variance, to compare the variances between the sampled over the territory that could possibly host the species, although not found bins younger and older than 1 Ma. We performed the Brown-Forsythe test by using there. The background points could be sampled randomly over the study area, but the R package “onewaytests” ver 2.629. this strategy implies that we know a-priori that the species could inhabit any area of Since Santini et al.59 warned about computing species’ geographic range by the territory being studied. Since this is far from a true species’ natural binarizing the habitat suitability maps according to threshold, we also computed distribution63, different strategies could be more suited to the ecology of the species for each time bin the surface area of cells in the maps with a specific range of being studied. We followed a strategy that simulates a random distribution of the suitability values. Specifically, we divided the past Ethiopian territories into two background points but taking into account two constraints: (1) the occurrence complementary spatial classes: one included cells with good-to-best habitat probability of background points is inversely proportional to the distance to pre- suitability values (GB, HSI values higher than 0.5) and another included cells with sence points; (2) the territory where we simulated the background points was not good-to-worst habitat conditions (GW, HSI values equals or lower than 0.5). Then, the whole of the area of present-day Ethiopia, but rather the territory enclosed by we counted the number of cells in each of the two spatial classes per time bin. Since the biomes actually inhabited by the Ethiopian wolf plus the geographic dis- our aim was to understand how many times in the last 2 Ma the species was in tribution of its main prey, the big-headed African mole rat Tachyoryctes macro- better and in worse conditions than at present, when a time bin showed a cephalus. To build this territory we used the world’s ecoregion polygons70. The percentage of GB class higher than the GW we interpreted this information as a range polygon of Tachyoryctes macrocephalus was downloaded from https://www. species geographic range expansion, whereas the reverse result indicated a stage of iucnredlist.org/. Once we have identified the polygons of the ecoregions inter- its range contraction. secting the Ethiopian wolf spatial polygons, we combined them with the polygon of This habitat suitability discrimination is important in order to understand the big-headed African mole rat. Then we simulated several background points under which bioclimatic conditions the species was able to survive. The procedure inside this unique combined polygon, verifying that there was only one background of HSI classification described above was applied to the portion of the territory of point per cell of the 4 × 4 km side resolution bioclimatic raster. We simulated Ethiopia where we reconstructed the species’ geographic range by binarizing the 10,000 background points59. habitat suitability map using the MUT strategy (described above). This geographic limitation was necessary because, in all of the time bins considered in this study, Statistical analyses. We used the predictions of the hindcasting and forecasting nearly the whole of the Ethiopia territory rendered very small suitability values, this models to reconstruct the geographic range of the Ethiopian wolf through time and biasing the relative abundance of the two HSI classes. to estimate how much of the conditions suitable for the species habitat was present in the past and will be available to the species in the future scenarios. The relationship among the species’ geographic range evolution, the temperature changes, and the species’ altitude shift in the past. Since the Ethiopian wolf is known Species geographic range and habitat suitability reconstructions in the past. At first, to be adapted to the Afroalpine habitat72 we tested the relationships between the after computing the continuous habitat suitability maps for each time interval from changes in the species’ geographic range through time with the temperature var- 0 to 2 Ma, we extracted the HSI values corresponding to the geographic coordinates iation in the last 2 Ma and the altitudes at which the species dwelt as response to of the Melka Wakena fossil site and to the territory of the whole Ethiopia for climate change. Indeed, since the species currently inhabits high-altitude territories, comparisons. Then, we reported the time bins for which we recorded the highest we hypothesized that any size variation of past geographic ranges was determined and lowest HSI values for Melka Wakena site during its estimated age range. by the species movement to higher elevations during unfavorable conditions during For estimating the species geographic range during past time bins, it is warm climates, and towards lower altitudes during periods of colder climates, when important to remember that past bioclimatic rasters were provided by Raia et al.55 C. simensis expanded its territories. In order to obtain a value for average at a spatial resolution of 50 km side square cell—a rather coarse resolution temperature for each time bin, we extracted from the territory of Ethiopia the compared to the 4 × 4 km cell raster used to train the models (see above). Hence, raster values of the annual maximum and minimum temperatures provided in when comparing past and present habitat reconstructions, we used the projections Raia et al.55. Then, for each time bin, we averaged these two measures of tem- we made of the models onto the present bioclimatic layers as provided by Raia perature and get an annual mean temperature in Ethiopia related to time bins. For et al.55. This allowed us to perform an informative comparison between the present the geographic range of the species, we used the surface area of the HSI maps and past suitable territories available to the species. At first, we reconstructed the binarized by the MUT strategy, whereas for the altitude, we recorded the mean temporal evolution of the species’ geographic range by binarizing the habitat altitude value obtained by intersecting the species’ geographic range with the suitability map obtained for each temporal interval. The binarization is obtained by altitude raster map of the Ethiopia territory as provided by http://www.worldclim. considering all the cells in the raster map with a suitability index higher than a org and rescaled to the same spatial resolution of Raia et al.55 bioclimatic raster COMMUNICATIONS BIOLOGY | (2023)6:530 | https://doi.org/10.1038/s42003-023-04908-w | www.nature.com/commsbio 11 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-04908-w maps. Then, we regressed the log10 of species geographic range surface area versus The following analyses were performed for each of the considered CMIP6 both the mean temperature values and log10 mean altitude per time bin. future climatic scenario projections, separately, since these models come from As ours are time series data, we had to account for issues related to temporal different starting hypotheses and physical processes the laboratories considered for auto-correlation in simple linear regressions’ residuals in order to have an accurate their modeling experiments. Then we used the landscape metrics to analyze the calculation of the regression coefficients73. This was done by applying an Auto temporal variation (including the metrics of the present time interval) of the mean Regressive Integrated Moving Average model (ARIMA)73. This model includes patch surface area (“area_mn”), while accounting for the differences between the three different components to properly compute the regression’s coefficients. The HSI classes. More specifically, we performed four simple linear regression models, auto regressive (AR) component allows to reduce the effect of auto-correlation in one for each of the CMIP6 future projections, in which we considered the Log10 of the series of observations, thus complying with the simple linear regression premise the mean patch surface area as a response variable, the temporal step of the models’ of independence of observations. Specifically, we calculated the order of the AR projections (from the present to the future) and the HSI classes as explanatory component that tells the ARIMA model how many consecutive values in the time variables. series affect any single observation. For instance, an AR component of order 3 Further, we tested the relationships between the four HSI classes and the means that any observation can be predicted by considering the values of the 3 landscape metrics to investigate if patches of species habitat suitability classes are preceding observations in the time series73. The Integrated component of the actually exploitable by the species in the present, and whether they will be so in the model (I) allows detrending the time series in order to have stationarity in mean future. We set four Generalized Linear Model (GLM) regressions in which the HSI and variance of the measurements, whereas the moving average (MA) component classes were the response variables, and all the landscape metrics were the allows smoothing the series values to reduce the outliers’ effect. At first, we explanatory variables. In these regression models, in order to take into account the performed an ordinary least square (OLS) regression between response (Log10 temporal variation of the variables, all the explanatory variables were considered in geographic range) and predictor variables (mean temperature and Log10 altitude). interaction with the time bins of the climatic projections, from the present to the Then, we measured the degree of auto-correlation in model residuals. In the case of future. We emphasize that, since C. simensis exists in a very restricted territory of significant autocorrelations between these residuals, we used ARIMA models to the Ethiopian highlands, most of the geographic surface considered in this study properly compute the relationship between the response and predictor variables provides unfavorable bioclimatic conditions. Most of the area of Ethiopia shows described above. We anticipate here that the residuals of the OLS showed low HSI values (falling within HSI class 1; see above), i.e., are barely suitable for C. significant autocorrelations. simensis. This overabundance of raster pixels of low HSI values occurred in the We carried out two different strategies to set the parameters of ARIMA model map of present conditions and in any of the considered future climatic scenarios components. In the first strategy, we used different combinations of parameters and might have masked any between-CMIP6 model differences and statistical based on the OLS residuals time series. Only the AR component order of the OLS relationships. Therefore, we excluded the HSI class 1 values from PerMANOVA residuals was estimated automatically by means of maximum likelihood method. and regression analyses. Moreover, since some HSI classes totally disappeared As regards the I component, since the time series of OLS residuals showed a clear (area_mn= 0) for two CMIP6 future projections, it was impossible to compute auto-regressive component, we also used the first order (1) difference between landscape metrics for these cases and data for the analyses were unbalanced. models’ residuals to confer stationarity to the residuals’ distribution73. As regards Therefore, we first standardized all data in a way such that all the CMIP6 model- the MA component, we performed six ARIMA models, each with fixed AR and I related projections had the same HSI classes computed and then performed all the components computed as described above, but with different MA component statistical analyses. We did not apply this correction to the four regressions values ranging from 0 (no moving average window) to 6 (i.e., by averaging up to 6 describing the relationship between mean patch area and time and HSI classes, consecutive values in the time series). In the second strategy we estimated all the since these analyses do not require the computation of the other landscape metrics. three ARIMA components by means of maximum likelihood estimation as It is important to note that, for all the same analyses, we did not consider the implemented in auto.arima function in the package “forecast” ver 8.1674 of the landscape metric “division” since it showed zero variance. software R42. After computing all the ARIMA models described in the two We emphasize that our present and future habitat suitability projections as well strategies, we chose the model with the lowest AIC. as the computation of the landscape metrics were not conducted with the aim of providing any strategy for the species conservation plans. They did turn out to Future suitable climatic conditions estimation. As mentioned above, since Santini highlight the acute threat that the Ethiopian wolf is already subjected to and the et al.59 discouraged the use of a unique binarized suitability map to estimate a species’ growing risks it would be exposed to in the near future. range size, we considered intervals of different habitat suitability values separately to analyze the temporal evolution of some landscape metrics characterizing the terri- Statistics and reproducibility tories suitable to C. simensis within the study area. For the present and for each of the Morphometric study of the fossil. Morphometric analyses were performed with SPSS future time bin maps, we considered four classes of habitat suitability value pixels: the statistics program ver. 25, and with Rv. 3.6.042 in RStudio ver. 2021.09.1 + 372 first class, with the worst HSI values (HSI < 0.25); the second, with low values “Ghost Orchid” (Release 8b9ced188245155642d024aa3630363df611088a, 2021-11- (HSI ≥ 0.25 and HSI < 0.5); a third one with good values (HSI ≥ 0.5 and HSI < 0.75); 08) for Windows using the packages “ggfortify” (ver. 0.4.013)43,44, “ggplot2” (ver. and a fourth class with the best habitat suitability conditions (HSI ≥ 0.75). For each of 3.3.5)45, “pvclust” (ver. 2.2-0)46,47, and “stats” (ver. 4.1.2.)42. the habitat suitability classes, we computed the following landscape structural con- nectivity indices. We did not compute any measure of distance between patches since, Modeling the bioclimatic niche through time given that C. simensis occupies mountain tops, differences of HSI values between near Simensis polygon: We simulated 200 Canis simensis occurrences over the grid of patches may be strongly affected by the steep altitudinal variation. The indices used in the bioclimatic variables by considering the geographic range polygon downloaded this analysis were: the number of patches (“np”); mean patch surface area from IUCN (https://www.iucnredlist.org/species/3748/10051312) as spatial (“area_mn”); Patch aggregation (“ai”) is the ratio between the number of the adjacent boundaries. For the present and future time intervals, we used the bioclimatic patches of the same class and the maximum possible number of like adjacencies variables at the spatial resolution of 1 × 1 and 4 × 4 km wide cell side (www. involving the corresponding class, which is achieved when the class is maximally wrldclim.org) combined with a rescaled raster of the Net Primary Productivity clumped into a single, compact patch. Its values range from 0 to 100; Clumpiness (https://cmr.earthdata.nasa.gov/search/concepts/C1631984056-LPDAAC_ECS. (“clumpy”) indicates how many equal-class patches are spatially adjacent if compared html), whereas for the past time bins, we used the variables as provided by Raia to a totally random spatial distribution of all the available different classes and its et al.55 with the spatial resolution of 50 × 50 km wide cell side. To train and values can range from −1 to 1; Cohesion (“cohesion”) measures the physical con- calibrate Maxent algorithm for the Species Distribution Modeling, we used the nectedness of the corresponding patch class (values ranging from 0 to 100); the “ENMeval” v. 2.0.369. Division index (“division”, ranging from 0 to 1) is the probability that two randomly The ARIMA regressions used for investigating the relationship among the chosen places in the landscape are not situated in the same undissected patch75. We species’ geographic range evolution, the temperature changes, and the species’ computed landscape metrics by using the R package “landscapemetrics” ver 1.5.476. altitude shift in the past were performed by using the “forecast” package v. 8.1674. We then used these landscape metrics to test the statistical differences between Habitat fragmentation metrics were performed by using the “landscapemetrics” the CMIP6 future climatic projections while considering the related four shared package v. 1.5.476. socio-economic pathways (ssp) as repeated experiments. For this, we set a two-way We performed two-way PerMANOVA analysis by using the package “vegan” v. permutational MANOVA (two-way PerMANOVA) in which the landscape metrics 2.6-2 and the related pairwise version with the “pairwiseAdonis” package ver 0.478. were the explanatory variables, the four CMIP6 models and the related future time All the statistical analyses were performed with the R statistical software v. bins were the factors and the four shared socio-economic pathways (ssp) were 4.1.042. experiment’s stratifications, i.e., repeated measures. PerMANOVA allows testing differences between groups with small sample size and without the assumption of normal distribution and equality of variance by relying on the computation of Reporting summary. Further information on research design is available in the Nature between-groups distance matrix and permutations. We considered the Log10 of all Portfolio Reporting Summary linked to this article. the used landscape metrics and computed the between-groups distance matrices by using the Euclidean distance. p-values were computed after setting 9999 permutations. Where significant differences were detected, a pairwise version of the two-way PerMANOVA with stratification was applied. We performed two-way Data availability PerMANOVA by using the R package “vegan” ver 2.6-277 and the related pairwise The codes and all the information to run the analyses with ecological data and landscape version by means of the “pairwiseAdonis” package ver 0.478. metrics for present, past and future time intervals, are freely available and stored at 12 COMMUNICATIONS BIOLOGY | (2023)6:530 | https://doi.org/10.1038/s42003-023-04908-w | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-04908-w ARTICLE https://zenodo.org/badge/latestdoi/618339598. The codes run under the free R statistical 28. Bonnefille, R., Melis, R. & Mussi, M. in The Emergence of the Acheulean in software and include the lines to install the necessary packages to perform the analyses. East Africa and Beyond. Vertebrate Paleobiology and Paleoanthropology (eds. Gallotti, R. & Mussi, M.) (Springer, 2018). Received: 1 December 2022; Accepted: 3 May 2023; 29. Dag, O., Dolgun, A. & Konar, N. M. Onewaytests: an R package for one-way tests in independent groups designs. R. J. 10, 1 (2018). 30. Mondanaro, A. et al. The role of habitat fragmentation in Pleistocene megafauna extinction in Eurasia. Ecography 44, 1619–1630 (2021). 31. Mooney, J. A., Marsden, C. D., Yohannes, A., Wayne, R. K. & Lohmueller, K. E. Long-term small population size, deleterious variation, and altitude adaptation in the ethiopian wolf, a severely endangered canid. Mol. Biol. Evol. References 40, msac277 (2022). 1. Yalden, D. W. The extent of high ground in Ethiopia compared to the rest of 32. Abrahms, B., Rafiq, K., Jordan, N. R. & McNutt, J. W. Long-term, climate- Africa. Sinet: Ethiop. J. Sci. 6, 35–38 (1983). driven phenological shift in a tropical large carnivore. Proc. Natl Acad. Sci. 2. Malcolm, J. R. & Ashenafi, Z. T. In The Ethiopian Wolf: Status Survey And USA 119, e2121667119 (2022). Conservation Action Plan (eds. Sillero-Zubiri, C. & Macdonald, D. W.) 61–63 33. Weber, M. M., Stevens, R. D., Diniz‐Filho, J. A. F. & Grelle, C. E. V. Is there a (IUCN/SSC Canid Specialist Group, 1997). correlation between abundance and environmental suitability derived from 3. Sillero-Zubiri, C., Marino, J., Gottelli, D. & Macdonald, D. W. In The Biology ecological niche modelling? A meta‐analysis. Ecography 40, 817–828 (2017). and Conservation of Wild Canids (eds. Sillero-Zubiri, C. & Macdonald, D. W.) 34. Noto, C. R. In Taphonomy 287–336 (Springer, 2011). (Oxford University Press, 2004). 35. Carotenuto, F. et al. MInOSSE: a new method to reconstruct geographic 4. Marino, J. & Sillero-Zubiri, C. Canis simensis. The IUCN Red List of ranges of fossil species. Meth. Ecol. Evol. 11, 1121–1132 (2020). Threatened Species Version 20112. https://doi.org/10.2305/IUCN.UK.2011-1. 36. Geraads, D. A revision of the fossil Canidae (Mammalia) of North-western RLTS.T3748A10051312.en (2011). Africa. Palaeontology 54, 429–446 (2011). 5. Gottelli, D., Sillero-Zubiri, C., Marino, J., Funk, S. M. & Wang, J. Genetic 37. Werdelin, L. & Lewis, M. E. Koobi Fora Research Project, Volume 7, the structure and patterns of gene flow among populations of the endangered Carnivora 333p (California Academy of Sciences, 2013). Ethiopian wolf. Anim. Conserv 16, 234–247 (2013). 38. Werdelin, L., Lewis, M. E. & Haile-Selassie, Y. New species of Eucyon 6. Castelló, J. R. Canids of the World: Wolves, Wild Dogs, Foxes, Jackals, Coyotes, (Mammalia; Carnivora; Canidae) from the Pliocene of the Woranso-Mille and Their Relatives Vol. 135 (Princeton University Press, 2018). Area, Afar Region, Ethiopia, and a critical review of African species of Eucyon. 7. EWCP Annual Report 2020. www.wildcru.org/news/ewcp-annual-report- Pap. Palaeontol. 1, 33–40 (2015). 2020-now-out (2020). 39. Ewer, R. F. The fossil carnivores of the Transvaal caves: Canidae. Proc. Zool. 8. Atickem, A. & Stenseth, N. C. The role of rodents in the conservation of Soc. London 126, 97–120 (1956). endangered species in the Ethiopian highlands. Therya 13, 73–77 (2022). 40. Adams, J. A., Olah, A., McCurry, M. R. & Potze, S. Surface model and 9. Traill, L. W., Bradshaw, C. J. & Brook, B. W. Minimum viable population size: tomographic archive of fossil primate and other mammal holotype and a meta-analysis of 30 years of published estimates. Biol. Conserv. 139, 159–166 paratype specimens of the Ditsong National Museum of Natural History, (2007). Pretoria, South Africa. PLoS ONE 10, e0139800 (2015). 10. Ferràndez-Cañadell, C., Ribot, F. & Gibert, L. New fossil teeth of 41. A. von den Driesch, A guide to the measurement of animal bone from Theropithecus oswaldi (Cercopithecoidea) from the Early Pleistocene at Cueva archaeological sites. Peabody Museum Bulletin 1 (Peabody Museum of Victoria (SE Spain). J. Hum. Evol. 74, 55–66 (2014). Archaeology and Ethnology, 1976). 11. Gottelli, D. et al. Molecular genetics of the most endangered canid: the 42. R. Core Team R: A language and environment for statistical computing. R Ethiopian wolf Canis simensis. Mol. Ecol. 3, 301–312 (1994). Foundation for Statistical Computing, Vienna, Austria. https://www.R-project. 12. Gottelli, D., Marino, J., Sillero‐Zubiri, C. & Funk, S. M. The effect of the last org/ (2021). glacial age on speciation and population genetic structure of the endangered 43. Horikoshi, M. & Tang, Y. ggfortify: Data visualization tools for statistical Ethiopian wolf (Canis simensis). Mol. Ecol. 13, 2275–2286 (2004). analysis results. https://CRAN.R-project.org/package=ggfortify (2016). 13. Hovers, E. et al. The expansion of the Acheulian to the Southeastern Ethiopian 44. Tang, Y., Horikoshi, M. & Wenxuan, L. ggfortify: unified interface to visualize Highlands: insights from the new early Pleistocene site-complex of Melka statistical result of popular R packages. R. J. 8, 478–489 (2016). Wakena. Quat. Sci. Rev. 253, 106763 (2021). 45. Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer-Verlag 14. Lindblad‐Toh, K. et al. Genome sequence, comparative analysis, and New York, 2016). haplotype structure of the domestic dog. Nature 438, 803–819 (2005). 46. Suzuki, R. & Shimodaira, H. Pvclust: an R package for assessing the 15. Gopalakrishnan, S. et al. Interspecific gene flow shaped the evolution of the uncertainty in hierarchical clustering. Bioinformatics 22, 1540–1542 (2006). genus. Canis. Curr. Biol. 28, 3441–3449 (2018). 47. Suzuki, R., Terada, Y. & Shimodaira, H. pvclust: hierarchical clustering with 16. Ciucani, M. M. et al. Evolutionary history of the extinct Sardinian dhole. Curr. p-values via multiscale bootstrap resampling. R package version 2.2-0. https:// Biol. 31, 5571–5579.e6 (2021). CRAN.R-project.org/package=pvclust (2019). 17. Perri, A. R. et al. Dire wolves were the last of an ancient NewWorld canid 48. Sokal, R. R. & Michener, C. D. A statistical method for evaluating systematic lineage. Nature 591, 87–91 (2021). relationships. Univ. Kans. Sci. Bull. 28, 1409–1438 (1958). 18. Azzaroli, A. Quaternary mammals and the end-villafranchian dispersal event 49. DePatta Pillar, V. How sharp are classifications? Ecology 80, 2508–2516 (1999). —a turning point in the history of Eurasia. Palaeogeo. Palaeoclim. Palaeoecol. 50. Elith, J. & Leathwick, J. R. Species distribution models: ecological explanation 44, 117–139 (1983). and prediction across space and time. Ann. Rev. Ecol. Evol. Syst. 40, 677–697 19. Tedford, R. H., Wang, X. & Taylor, B. E. Phylogenetic systematics of the (2009). North American fossil caninae (Carnivora: Canidae). Bull. Am. Mus. Nat. 51. Araújo, M. B. et al. Standards for distribution models in biodiversity Hist. 2009, 1–218 (2009). assessments. Sci. Adv. 5, eaat4858 (2019). 20. Bartolini-Lucenti, S., Bukhsianidze, M., Martínez-Navarro, B. & 52. Chiarenza, A. A., Mannion, P. D., Farnsworth, A., Carrano, M. T. & Varela, S. Lordkipanidze, D. The wolf from Dmanisi and augmented reality: review, Climatic constraints on the biogeographic history of Mesozoic dinosaurs. implications, and opportunities. Front. Earth Sci. Paleont. 8, 131 (2020). Curr. Biol. 32, 570–585 (2022). 21. Martinez-Navarro, B. et al. A new species of dog from the Early Pleistocenesite 53. Cunze, S. & Klimpel, S. From the Balkan towards Western Europe: Range of Venta Micena (Orce, Baza Basin, Spain). C. R. Palevol. 20, 297–314 (2021). expansion of the golden jackal (Canis aureus). A climatic niche modeling 22. Bartolini-Lucenti, S. et al. The early hunting dog from Dmanisi with approach. Ecol. Evol. 12, e9141 (2022). comments on the social behaviour in Canidae and hominins. Sci. Rep. 11, 54. Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1km spatial resolution climate 13501 (2021). surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017). 23. Martínez-Navarro, B. & Rook, L. Gradual evolution in the African hunting 55. Raia, P. et al. Past extinctions of Homo species coincided with increased dog lineage systematic implications. C. R. Palevol. 2, 695–702 (2003). vulnerability to climatic change. One Earth 3, 480–490 (2020). 24. Lordkipanidze, D. et al. Postcranial evidence from early Homo from Dmanisi, 56. Holden, P. B. et al. PALEO-PGEM v1. 0: a statistical emulator of Georgia. Nature 449, 305–310 (2007). Pliocene–Pleistocene climate. Geosci. Model Dev. 12, 5137–5155 (2019). 25. Lordkipanidze, D. et al. A complete skull from Dmanisi, Georgia, and the 57. Meinshausen, M. et al. The shared socio-economic pathway (SSP) greenhouse evolutionary biology of early. Homo. Sci. 342, 326–331 (2013). gas concentrations and their extensions to 2500. Geosci. Model Dev. 13, 26. Xue, L., Alemu, T., Gani, N. D. & Abdelsalam, M. G. Spatial and temporal 3571–3605 (2020). variation of tectonic uplift in the southeastern Ethiopian Plateau from 58. IPCC. Climate Change 2014: Synthesis Report. Contribution of Working morphotectonic analysis. Geomorphology 309, 98–111 (2018). Groups I, II and III to the Fifth Assessment Report of the Intergovernmental 27. Marino, J. Threatened Ethiopian wolves persist in small isolated Afroalpine Panel on Climate Change (eds. Core Writing Team, R. K. Pachauri, R. K. & enclaves. Oryx 37, 62–71 (2003). Meyer, L. A.) 151pp (IPCC, 2014). COMMUNICATIONS BIOLOGY | (2023)6:530 | https://doi.org/10.1038/s42003-023-04908-w | www.nature.com/commsbio 13 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-04908-w 59. Santini, L., Benítez‐López, A., Maiorano, L., Čengić, M. & Huijbregts, M. A. California, Berkeley (T.G.); this contribution is framed within a wider project on late Assessing the reliability of species distribution projections in climate change Neogene vertebrate evolution granted with “Fondi di Ateneo” by the Earth Sciences research. Divers. Distrib. 27, 1035–1050 (2021). Department of the University of Florence (L.R., S.B.L.); the field geological and strati- 60. Elith, J., Kearney, M. & Phillips, S. The art of modelling range‐shifting species. graphic work was supported by the National Geographic Society (Grant #9801-15) and Meth. Ecol. Evol. 1, 330–342 (2010). the Fritz Thyssen Foundation (E.H), and the School of Earth Sciences, Addis Ababa 61. Zuur, A. F., Leno, E. N. & Elphick, C. S. A protocol for data exploration to University (A.A.); the Spanish Ministry of Science (grants CGL2017-92166-EXP and avoid common statistical problems. Meth Ecol. Evol. 1, 3–14 (2010). PID2019-111185GB-I00) and The Junta de Andalucía (grants P18-FR-3193 and 62. Breiner, F. T., Guisan, A., Bergamini, A. & Nobis, M. P. Overcoming UMA18-FEDERJA-188) (B.F., P.P.); the Ann and Gordon Getty and John Templeton limitations of modelling rare species by using ensembles of small models. Foundations, the National Science Foundation (Grant BCS-2020044) (P.P.R., E.M.N.); Meth. Ecol. Evol. 6, 1210–1218 (2015). and the Braginsky Center for the Interface between Science and the Humanities at the 63. Phillips, S. J., Anderson, R. P. & Schapire, R. E. Maximum entropy modeling Weizmann Institute of Science, the Israel Academy of Science and Humanities COVID- of species geographic distributions. Ecol. Mod. 190, 231–259 (2006). 19 emergency postdoctoral fellowship, and the Israel Science Foundation research grant 64. Merow, C., Smith, M. J. & Silander, J. A. Jr A practical guide to MaxEnt for 1745/21 awarded to U. Smilansky (G.H.). modeling species’ distributions: what it does, and why inputs and settings matter. Ecography 36, 1058–1069 (2013). 65. Elith, J. et al. A statistical explanation of MaxEnt for ecologists. Divers. Distrib. Author contributions B.M.N., E.H., and T.G. designed the study; E.H. and T.G. supervised the field study at the 17, 43–57 (2011). MW site-complex; E.H, T.G., and G.H. conducted the archaeological excavation; T.G. 66. Jaynes, E. T. Information theory and statistical mechanics. Phys. Rev. 106, 620 found the fossil specimen; A.A. carried out the field work and analyzed the stratigraphic (1957). and sedimentological data; B.M.N. cleaned the specimen and carried out the anatomical 67. Shannon, C. E. A mathematical theory of communication. Bell Syst. Tech. J. description of the fossil hemimandible; S.B.L., P.P., B.F., and L.R. carried out the com- 27, 379–423 (1948). parative morphological and morphometrical study of the specimen; P.R.R. and E.M.N. 68. Helmstetter, N. A., Conway, C. J., Stevens, B. S. & Goldberg, A. R. Balancing carried out the geochronological study; F.C. conducted the diachronic modeling of the transferability and complexity of species distribution models for rare species Canis simensis bioclimatic niche; B.M.N., E.H., T.G., F.C., S.B.L., A.A., and P.P. wrote the conservation. Divers. Distrib. 27, 95–108 (2021). manuscript. All the authors participated in the discussion and agreed to the final version 69. Kass, J. M. et al. ENMeval 2.0: redesigned for customizable and reproducible of the manuscript. modeling of species’ niches and distributions. Meth. Ecol. Evol. 1602–1608, https://doi.org/10.1111/2041-210X.13628 (2021). 70. Dinerstein, E. et al. An ecoregion-based approach to protecting half the Competing interests terrestrial realm. BioScience 67, 534–545 (2017). The authors declare no competing interests. Borja Figueirido is an Editorial Board 71. Scherrer, D., D’Amen, M., Fernandes, R. F., Mateo, R. G. & Guisan, A. How to Member for Communications Biology since April 1st, 2023, but was not involved in the best threshold and validate stacked species assemblages? Community editorial review of, nor the decision to publish this article. optimisation might hold the answer. Meth. Ecol. Evol. 9, 2155–2166 (2018). 72. Marino, J., Sillero‐Zubiri, C. & Macdonald, D. W. Trends, dynamics and resilience of an Ethiopian wolf population. Anim. Conserv 9, 49–58 (2006). Additional information 73. Cowpertwait, P. S. & Metcalfe, A. V. Introductory Time Series with R (Springer Supplementary information The online version contains supplementary material Science & Business Media, 2009). available at https://doi.org/10.1038/s42003-023-04908-w. 74. Hyndman, R. J. & Khandakar, Y. Automatic time series forecasting: the forecast package for R. J. Stat. Softw. 27, 1–22 (2008). Correspondence and requests for materials should be addressed to Bienvenido Martínez- 75. McGarigal, K. FRAGSTATS: Spatial Pattern Analysis Program for Categorical Navarro, Tegenu Gossa or Erella Hovers. Maps. Computer software program produced by the authors at the University of Massachusetts, Amherst. http://www.umass.edu/landeco/research/fragstats/ Peer review information Communications Biology thanks Pasquale Raia, Lars Werdelin, fragstats.Html (2002). and Ignacio Lazagabaster for their contribution to the peer review of this work. Primary 76. Hesselbarth, M. H., Sciaini, M., With, K. A., Wiegand, K. & Nowosad, J. Handling Editor: Luke R. Grinham. Landscapemetrics: an open‐source R tool to calculate landscape metrics. Ecography 42, 1648–1657 (2019). Reprints and permission information is available at http://www.nature.com/reprints 77. Oksanen, J. et al. vegan: Community Ecology Package. R package version 2.5- Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in 6. 2019 (2020). published maps and institutional affiliations. 78. Martinez Arbizu, P. pairwiseAdonis: Pairwise multilevel comparison using adonis. R package version 0.0, 1 (2017). Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, Acknowledgements adaptation, distribution and reproduction in any medium or format, as long as you give We thank the Ethiopian Heritage Authority (EHA), the Culture and Tourism Bureaus of appropriate credit to the original author(s) and the source, provide a link to the Creative the Oromia Regional State, the West Arsi zone, and the Gadeb-Hasasa woreda for their Commons license, and indicate if changes were made. The images or other third party support during the field and post-excavation phases of the study. We are grateful to the material in this article are included in the article’s Creative Commons license, unless people of Aluba village for their hospitality and help in the field. We thank B. Van indicated otherwise in a credit line to the material. If material is not included in the Valkenburgh for facilitating us the measurements of several canid mandibles housed in article’s Creative Commons license and your intended use is not permitted by statutory different institutions. Finally, we thank P. Raia, L. Werdelin, and I. A. Lazagabaster for regulation or exceeds the permitted use, you will need to obtain permission directly from their insightful comments that helped to improve the manuscript. This research was the copyright holder. To view a copy of this license, visit http://creativecommons.org/ supported by the Palarq Foundation, the Spanish Ministry of Science and Innovation licenses/by/4.0/. through the “María de Maeztu” excellence accreditation (CEX2019-000945-M) and the Generalitat de Catalunya, 2021SGR 01238 (AGAUR) (B.M.N.); the Gerda Henkel Foundation and Glynn Isaac post-doctoral Fellowship of HERC at the University of © The Author(s) 2023 14 COMMUNICATIONS BIOLOGY | (2023)6:530 | https://doi.org/10.1038/s42003-023-04908-w | www.nature.com/commsbio