Niche diversification of Mediterranean and southwestern Asian tortoises

Background Tortoises of the genus Testudo are widely distributed throughout the Mediterranean region and southwestern Asia. However, the evolutionary mechanisms of diversification in this genus are still poorly understood. Methods In this study, we assessed the evolutionary patterns in the climate niches of five species and 11 subspecies of the genus Testudo using ecological niche models and evaluated the niche overlap based on species phylogenetic distances. Results The ecological models indicated that most species differ in their climate niches, but show overlap, with gradual transitions at range boundaries. As expected, the ecological divergence among subspecies was lower than that among species. Evaluation of the phylogenetic signal indicated that climate niches have been weakly conserved, but sister species also show high evolutionary divergence.


INTRODUCTION
Closely related species tend to respond similarly to environmental variation because they share similar life histories and physiological traits (Cruz et al., 2011;Dahlke et al., 2020). This lineage-specific evolutionary inertia in environmental response is known as phylogenetic niche conservatism (Crisp & Cook, 2012). Niche conservatism is common during the speciation process of a lineage, and it is produced by the fragmentation of ancestral subpopulations in environmentally analogous regions (Peterson & Nyari, 2008). Climatic oscillations, including glacial phases, promote niche diversification, as some subpopulations are subject to substantial contractions in their ranges during unfavourable periods, reducing or preventing gene flow (Leroy et al., 2017;Andersen et al., 2019).
A vicariant mode of speciation is evident in several lineages of Mediterranean reptiles (Kindler et al., 2017;Senczuk et al., 2017). Most reptiles are ectothermic, having entire reliance on external heating sources to control body temperature and complete embryonic development (Brattstrom, 1965;Booth, 2006). Therefore, the temperature is a major environmental factor influencing species richness among reptile assemblages in temperate regions (Rodríguez, Belmontes & Hawkins, 2005). The oscillatory phases throughout the Pleistocene, which involved a general drop in temperature and increased aridification (Magri & Parra, 2002), had profound impacts on the process of reptile diversification and distribution (Araújo et al., 2008;Fitze et al., 2011;Anadón et al., 2015). However, the influence of recent climatic history on the diversification process of Mediterranean tortoises is little known.
In this study, we tested for phylogenetic climatic niche conservatism and modelled the ranges of the tortoise species of the genus Testudo. This genus includes species that are morphologically similar (size ranges between 20 and 35 cm) and are generalist herbivores (Bonin, Devaux & Dupré, 2006). The genus Testudo comprises five species and 15 subspecies, distributed throughout Eurasia and North Africa in open and relatively warm environments (Calzolai & Chelazzi, 1991;Fritz et al., 2009;Escoriza et al., 2022). Most Testudo species are either parapatrically or allopatrically distributed, although some species are sympatric in a partial range of their distribution in the southern Balkans, such as T. graeca-T. hermanni and T. hermanni-T. marginata (Valakos et al., 2008). This pattern of occurrence suggests: (i) that the climatic niche of the species is conserved, being isolated by competitive exclusion or environmental barriers (Pasch, Bolker & Phelps, 2013;Jankowski et al., 2013;Li, Shao & Li, 2020); or (ii) that the species are ecologically divergent but contact in narrow transitional areas within the margins of environmental tolerance of the two parapatric species (Jones et al., 2020). We tested these hypotheses using ecological niche models (ENMs).

Study region and species
The study area included the western Palaearctic biogeographic realm (southern Europe and northern Africa) and southwestern Asia between Pakistan and Turkey ( Fig. 1). Data on the occurrence of all the five species of the genus Testudo were obtained during several samplings in the region (Escoriza & Ben Hassine, 2017;Escoriza & Pascual, 2021), from bibliographic sources, the GBIF (Global Biodiversity Information Facility, https://www.gbif.org) and iNaturalist (only Testudo marginata; https://www.inaturalist.org/) databases (Appendix S1 and S2). These databases are to be reliable data sources for biodiversity research (García-Roselló et al., 2015;Hochmair et al., 2020). We selected from the GBIF those records added since 2000, with a fine spatial resolution (error < 1,000 m), far from cities and within the known ranges of the species, as defined by Bonin, Devaux & Dupré (2006). To reduce spatial autocorrelation bias (Boria et al., 2014), record points that are closer than 10 km from each other were removed from the analysis using spThin (Aiello-Lammens et al., 2015) in R (R Core Team, 2022). This process reduced the number of locations by 54% (from 964 to 519) (Appendix S1): T. graeca (274 records), T. hermanni (103), T. horsfieldii (77), T. kleinmanni (4), and T. marginata (61) (Fig. 1).

Environmental data
Data on 19 bioclimatic variables, available in World Clim2 (Fick & Hijmans, 2017), were obtained from the species occurrences. Model complexity was reduced by estimating the increase in the variance inflation factor (VIF) (Craney & Surles, 2002). To estimate the VIF, we defined buffer polygons of 200 km around the location of each species and generated 1000 random points within these polygons. The climate data obtained for the random points and the location of each species were included in a logistic regression and used to estimate the variations in VIF for each new variable added to the model. We started with a simple model including only bio 10 (mean temperature of the warmest quarter), bio 11 (mean temperature of the coldest quarter), and bio 12 (annual precipitation). These variables were initially chosen because of their relevance in the life cycles of tortoises, and they included summer and winter temperatures (specifically associated with embryonic development and hibernation) and moisture availability (Lambert, 1983;Anadón et al., 2012). Variables with VIF > 10 were excluded (Schroeder, Lander & Levine-Silverman, 1990).

Data analyses
We visualized the realised niche of each species using principal component analysis (PCA) of the normalised variables. The climate niche of each species was generated using Maxent, as this method is suitable for models based only on presence data (Elith et al., 2011). The best Maxent model was selected by comparing several candidates built using various features (L: linear; Q: quadratic; H: hinge; P: product; T: threshold) and regularisation multipliers (RMs). This iterative selection was conducted using the Akaike criterion corrected for finite samples (AICc), using the ENMeval functions (Muscarella et al., 2014) in R. The features and the RM for the optimal model were used to generate 30 replicates in Maxent 3.4.4 (Phillips, Anderson & Schapire, 2006), calibrated using 75% of the species occurrences. The reliability of these models was assessed using the AUC statistic and variable contributions based on permutations of the presence-background data (Phillips, Anderson & Schapire, 2006). The best ENMs were used to infer the location and extent of potential palaeoclimate refugia. To do this, we projected the current relationships between climate and species to the paleoclimatic conditions, using the default clamping function in Maxent (Zhao et al., 2021). We used several palaeoclimatic layers covering the period from the early Pleistocene (ca. 787 ka), and a complete glacial cycle that included the last interglacial (ca. 130 ka), the last glacial maximum (ca. 21 ka), and the mid-Holocene (ca. 6 ka) (Otto-Bliesner et al., 2006;Brown et al., 2018). The logistic result of these models was binarised using the fuzzify function in Quantum GIS (mid point = 0.5). The four layers generated (along with one generated for present conditions) were summed, producing a new layer in which values close to five identified those areas that were climatically suitable as refugia throughout most of the Pleistocene-Holocene periods. The niche of T. kleinmanni was not modelled, because the number of localities was low. We also evaluated the niche divergence between pairs of species and subspecies of the genus Testudo. These comparisons were conducted for all species and between parapatric conspecific subspecies. Those taxa having fewer than 19 occurrences were excluded from the niche tests (Jenkins & Quintana-Ascencio, 2020). Niche divergence was estimated using a suite of tests (overlap, linear and blob range-breaking, and identity and background tests) that estimated the D index, which varies between 0 (no overlap) and 1 (full overlap) (Schoener, 1970). The overlap test quantified the similarity between ENMs using Latin hypercube sampling (Warren et al., 2021). The range-breaking tests assess the occurrence of sharp changes in environmental conditions at the range boundaries of a species. The linear test assumes that this limit is linear, whereas the blob test assumes that the limit is geometrically irregular; the latter is more robust if the species occupy geographical areas differing in extent, and the sample sizes differ (Glor & Warren, 2011). Identity tests assessed whether the models generated for the occurrence of two species differed statistically from those generated from random subsamples of their occurrences (Warren, Glor & Turelli, 2008). Background tests evaluated the environmental conditions that surrounded a species range, enabling assessment of whether the niches of two species were more similar to each other than to the available conditions (Warren, Glor & Turelli, 2008). We defined background regions as buffers of 500 km for the species and 250 km for the subspecies and natural micro-endemic T. marginata. These background regions were defined according to the distribution of each species or subspecies, maximizing the inclusion of environmental variation around the occurrence areas while minimizing the effect of non-informative regions (such as temperate-boreal biomes or deserts, not occupied by any species of the genus) (Acevedo et al., 2012).
The identity and background tests were conducted based on ENMs (Warren, Glor & Turelli, 2008) and an ordination method (ecospat) (Broennimann et al., 2012). The ecospat method has the advantage that it reduces the chances of model overfitting (Broennimann et al., 2012). We also estimated the Spearman rank correlation between environmental variation and predicted suitability (Warren et al., 2021). This correlation varies from −1 to 1, with 0 representing random correlation (Warren et al., 2021). These tests were implemented using generalized linear models and their statistical significance was assessed after 500 replicates (Warren et al., 2021). These analyses were conducted using the ENMTools package (Warren, Glor & Turelli, 2010) in R.
The evolutionary relationships among several species and subspecies of the genus Testudo were taken from a published molecular phylogeny (Graciá et al., 2017) (Fig.  2). The phylogenetic signal was determined only for a set of taxa in which the extent of divergence was quantified for the same genes (Fig. 2). The test of phylogenetic influence over environmental niche occupancy (defined by the same subset of low-correlated variables) was determined using Blomberg's K (Blomberg, Garland Jr & Ives, 2003). Blomberg's K values vary between 0 and ∞, where values of K < 1 represent less phylogenetic signal than that expected under Brownian motion (Blomberg, Garland Jr & Ives, 2003). Blomberg's K provides an acceptable estimate of the phylogenetic signal even when sample sizes are limited, although is susceptible to false positives (Blomberg, Garland Jr & Ives, 2003;Münkemüller et al., 2012). The value of the phylogenetic metric was calculated after 10,000 resamplings of the climate matrix, thus avoiding the error associated with estimations based only on the mean value per species (Harmon & Losos, 2005). These calculations were performed using the phytools package (Revell, 2012) in R.

RESULTS
The scatter plot (Fig. 3) showed that the species occupy well-differentiated positions along the environmental gradient, from a mesophilic species (with high scores in the variables that describe precipitations), such as T. hermanni, other species that occupied an intermediate position (T. graeca, T. marginata), to xerophilic species such as T. horsfieldii and T. kleinmanni, that occupy one of the climate extremes (Fig. 3A). The subspecies show a tendency to overlap or occupy contiguous positions in the niche space (Fig. 3B), although the western group of T. graeca occupies a more limited niche than the eastern group of subspecies (Fig. 3B).
All Maxent models showed good performance (Appendix S2). The explanatory capacity of the climatic variables determining the species distributions differed significantly (Appendix S2). For T. graeca the variable having the greatest importance was annual precipitation (39.1%), for T. hermanni it was annual precipitation (77.9%), for T. horsfieldii it was the temperature of the wettest quarter and precipitation in the warmest quarter (29.1%), and for T. marginata it was precipitation in the warmest quarter (32.3%). The mapped results of the ENMs are shown in Fig. 4. The ENM projections identified several palaeoclimate refugia, which differed in most of the species (Fig. 5).
The results of the niche tests did not indicate that these species are niche conservative, both in the group of parapatric species and between sister species (T. hermanni-T. horsfieldii). The lack of statistical significance for most of the comparisons of the range breaking tests allows us to reject the hypothesis of sharp environmental variations at the species' range boundaries (Table 1). In general, climatic niche overlap between species was lower than between subspecies, and species tend to occur in distinct climate niches (Tables  1 and 2). Similarly, species showed negative values or close to 0 in their correlations between climate variation and predicted suitability, while in some pairs of subspecies (e.g., T. g. graeca and T. g. marokkensis) this correlation was positive, indicating similar environmental responses (Tables 1 and 2). The phylogenetic signal estimation indicated lability during the evolutionary processes for the genus Testudo, with absent or weak phylogenetic signal for the entire subset of bioclimatic variables included in the ENMs (i.e., BIO8-BIO12, BIO15, and BIO18) (K ≤ 0.475; Table 3).

DISCUSSION
The combined results indicate niche partitioning among the species of the genus Testudo. The temperate and humid region is largely occupied by T. hermanni and as the conditions become more arid, the species progressively replace one another, appearing two species entirely adapted to extreme climates, cold steppes (T. horsfieldii) or specialized in sand semi-deserts (T. kleinmanni) (Bringsøe & Buskirk, 1998;Valakos et al., 2008;Bondarenko & Peregontsev, 2017;Arakelyan et al., 2018). Consistently, we found that climate niches are not phylogenetically conserved, as the values for Blomberg's K in our analysis range from 0.263 to 0.475. Mapped projections of the ENMs showed marked differences in the climate niches of these species of tortoise, and niche overlap was generally low. The analyses also suggest that climate (and particularly annual precipitation, and precipitation and temperature of the warmest quarter) plays a key role in determining tortoise distribution in the study area, given the high predictive performance of the ENMs, as already indicated by previous studies involving the eastern clades of T. graeca (Turkozan, Karacaoğlu & Parham, 2021).
The analyses revealed that, in most cases, the geographical boundaries between these tortoise species and subspecies were not determined by rapid climate transitions. This implies that the species can partially differ in their niches, but also share zones of environmental tolerance, as shown by the ENM results. These regions (e.g., the southern Balkans or Sardinia) are environmentally suitable for several species, and natural populations (or populations resulting from historical introductions; Vamberger et al., 2011) of two or three species coexist. In transitional regions, the species turnover may be determined by negative interspecific interactions or habitat properties, but they can also occur sympatrically (Wright, Steer & Hailey, 1988;Valakos et al., 2008).
The ENMs and ecospat identity tests did not show complete consistency. The results of the ecospat method were more conservative (i.e., they rejected the null hypothesis of niche equivalency less frequently). We considered that the species are ecologically divergent if the results of the two identity tests were consistent. These tests showed a pattern of niche divergence for most of the species pairs, except for T. graeca-T. marginata. However, in this case, the environmental response was not convergent (i.e., the values of the rank correlation coefficient were close to 0), and they were no more conservative than expected given their available conditions (i.e., the background tests were not statistically significant). The subspecies showed greater similarities in their climate niches, with pairs occurring in equivalent climate niches and with convergent environmental responses (for example, the T. g. whitei-T. g. nabeulensis pair). Lower niche divergence was expected between subspecies since their genetic separation spans smaller timescales (Peterson, 2011;Meynard et al., 2017).
Niche divergence was especially marked for the eastern and western groups of T. graeca, as indicated by significant differences in the identity tests and negative values for the rank correlation coefficient. Indeed, western populations of T. graeca are confined to warm and dry environments in southern Europe, whereas the species extends to colder regions in eastern Europe (Anadón et al., 2006;Buică, Iosif & Cogălniceanu, 2013). This high niche divergence observed between both clades would initially support their separation as species, as proposed by Van der Kuyl, Ballasina & Zorgdrager (2005). The North African subclade has evolved isolated in the southern Mediterranean, in a hot and arid region on the fringes of the Sahara, which could have reduced its cold tolerance, compared to Asian subspecies. However, it should be evaluated whether these differences in the realized niches between the western and eastern subclades of T. graeca also extend to the fundamental niches, through experimental mechanistic studies (Jiménez et al., 2019;Bujan et al., 2022).
Complementary to the niche tests, evaluation of the phylogenetic signal suggested a poorly conserved climate niche. This result contrasts with data for other reptile groups from the region, including lacertid lizards, which show a strong phylogenetic signal in species climate niches (Escoriza, Pascual & Mestre, 2021). However the power of Blomberg's K tests is limited when using small phylogenies (Münkemüller et al., 2012), and future studies have to evaluate these evolutionary patterns with a larger number of tortoise species.

CONCLUSIONS
The speciation process of tortoises from the genus Testudo is closely linked to regional climate history. The distribution of subspecies can be explained by the presence of stable climate refugia throughout the Pleistocene. Most species occupy well-defined and separate climate niches, reducing the chances of contact and interspecific interactions. In the case of occupation of equivalent niches, the species may occur in separate geographic regions (e.g., the T. graeca-T. marginata).