Isolation by distance and past climate resistance shaped the distribution of genealogical lineages of a neotropical lizard

Organisms adapted to open environments in South America have recently been used to understand the origins of the high Neotropical biodiversity. In the Caatinga, the largest continuous block of Seasonally Dry Tropical Forest in South America, phylogeographic studies have uncovered the role of historical climate changes and rivers (i.e., the São Francisco River, the largest perennial river in Caatinga), in promoting genetic differentiation and speciation of lizards and amphibians. We used mitochondrial data, demographic analyses, paleodistribution models, and landscape genetic methods to test the effects of spatial distances, historical climate fluctuations, and landscape heterogeneity on the genetic variation of the generalist lizard Tropidurus hispidus in the semi-arid Caatinga in northeastern Brazil. Four haplogroups with moderate geographical structure diverged in the Pleistocene, and exhibited high genetic diversity. Ecological niche models revealed large suitable climatic areas for T. hispidus in the past 790 thousand years, connecting the Caatinga and other regions via a narrow corridor. Part of the genetic differentiation in T. hispidus resulted from spatial distances among populations and isolation by resistance through climatic unsuitability areas in the Last Glacial Maximum (LGM), which probably reduced population connectivity and gene flow. Our findings highlight the role of the historical factors of the Caatinga, through LGM climate, and the generalist condition of species in shaping the genealogical histories of populations. Although the results are based on a single-locus approach, our study is a first step to shed light on the main drivers of the evolutionary history of T. hispidus, in a highly diverse and still poorly studied region.


Introduction
The interplay between historical climatic changes, landscape evolution, and organismal traits, such as physiological tolerance and dispersion ability, shapes species' geographic distributions through time. Through geomorphological and climate dynamics, landscape history and heterogeneity can also modulate gene flow between metapopulations and the spatial structure of genetic variation (Albert et al., 2017;Wang & Bradburd, 2014). For example, physical attributes (e.g., topographic variation and the presence and specific characteristics of rivers) across the landscape may lead to populational genetic differentiation, even if not wholly halting gene flow (Fonseca et al., 2021;Rivera et al., 2020;S anchez-Ram ırez et al., 2018). Likewise, present-day climatic gradients and severe climatic oscillations in the recent past, such as Pleistocene climatic changes, should have influenced populations' connectivity and structure (Sexton et al., 2014;Vasconcellos et al., 2019). Combining population genetics and landscape heterogeneity may help us unravel the putative roles of environmental features on genetic variation and gene flow (Manel et al., 2003).
The development of new theoretical and analytical approaches that account for the potential association between landscape heterogeneity and genetic differentiation allows us to study microevolutionary processes that shape the current spatial distribution of genetic variation (Manel et al., 2003;McRae, 2006). The null model in landscape genetics, in contrast, considers geographical distances as the sole driver of genetic differentiation. In these cases, the inter-population restricted dispersal and genetic drift can produce patterns of isolation by distance (IBD), in which genetic and geographic distances are positively correlated (Slatkin, 1993;Wright, 1943). For several organisms, spatial distances are a significant predictor of genetic structure and IBD is present in most species (Jenkins et al., 2010;Perez et al., 2018). Besides, IBD is more likely for ectotherms, suggesting that natural history features, such as metabolic basis, can influence gene flow rates and genetic structure (Jenkins et al., 2010). Conversely, in natural systems, the landscape is more heterogeneous than straight-line distances can capture, presenting current or past environmental differences and potential barriers to dispersal that may impact directly genetic differentiation of populations (McRae, 2006;Sexton et al., 2014;Thom e, Carstens, Rodrigues, Galetti Jr, et al., Thom e et al., Thom e et al., 2021).
The degree of complexity of a given landscape can generate, among others, patterns of isolation by the environment (IBE) (Sexton et al., 2014;Wang & Bradburd, 2014) and isolation by resistance (IBR) (McRae, 2006;McRae & Beier, 2007). In the former, environmental differences mediate connectivity among populations, affecting species dispersion and, consequently, the patterns of spatial genetic divergence, apart from spatial distances (Sexton et al., 2014;Wang & Bradburd, 2014). IBE is the main predictor of genetic divergence for several tropical lizard species (Fenker et al., 2021). Complementarily, IBR predicts that genetic divergence is mediated by landscape features that impose physical restrictions on dispersal and connectivity among populations (McRae, 2006), acting as a biogeographic barrier (e.g., unsuitable habitats, complex topography, and rivers). However, species idiosyncrasies modulate the permeability of landscape features to gene flow (Zamudio et al., 2016). Whether IBD, IBE, or IBR (or their combinations) is the primary factor responsible for genetic structuring may depend on the species' lifehistory traits, such as dispersal ability (Paz et al., 2015).
Environmental variability and historical landscape changes shaped the spatial patterns of diversity in different scales in the Caatinga (e.g., from communities to population genetics), making it an excellent environment to study the roles of landscape features on diversification. For harsh environments such as the arid Caatinga, water availability should modulate spatial variation in species richness (Hawkins et al., 2003). Temperature and precipitation fluctuations, for example, influenced the distribution of amphibians and reptiles (de Oliveira & Diniz-Filho, 2010). However, at a finer scale, several other contemporary and historical attributes may have contributed to the current spatial distribution of the genetic diversity and structure. Some studies documented IBD, IBE, and IBR patterns among Caatinga bees, frogs, and lizards (Jaff e et al., 2019;Thom e, Carstens, Rodrigues, Galetti Jr, et al., 2021;. Annual temperature variation, forest cover, and topography variables explained most of the gene flow in the stingless bee Melipona subnitida (Jaff e et al., 2019). Also, the asynchrony of seasons explained genetic differentiation and mediated gene flow of populations in the toad Rhinella granulosa (Thom e, Carstens, Rodrigues, Galetti Jr, et al., 2021;. Conversely, both IBD and IBR predicted the genetic differentiation of the whiptail lizard Ameivula ocellifera, with the resistance of main rivers explaining most of the variation . Topographically, highlands such as the Espinhaço Mountain Range (Chapada Diamantina) and plateaus of the Caatinga could have influenced the evolutionary history of species in at least two non-excluding ways: a milder climate gradient generated by higher elevations (higher than 500 m, reaching up approximately 2,000 m), and a higher topographic complexity that can reduce the movement and connectivity among populations (Jablonski et al., 2016;Rodr ıguez et al., 2015). In addition, the Caatinga is highly threatened, being the target of intense conversions of its fragments through chronic anthropogenic modifications (Antongiovanni et al., 2020). Thus, species widely distributed across the biome offer an excellent opportunity to test the main drivers of diversification of a highly threatened and still poorly studied dry environment.
The Caatinga is a highly diverse region and an area of endemism for the genus Tropidurus, harbouring 13 of the 28 species described (Carvalho et al., 2013;Carvalho et al., 2016;Uetz et al., 2020). Tropidurus species use several habitat types (e.g., sandy soils, rocky outcrops, semi-arboreal, and open areas) and can have from very narrow (e.g., paleodunes from the São Francisco River, rock fields from Espinhaço Mountain Range) to broad geographic distributions (Mesquita et al., 2017;Rodrigues, 2005). In contrast to its endemic congeners, T. hispidus has a wide geographic distribution, covering the Caatinga and extending into coastal areas in northeastern Brazil and Amazonian savannas enclaves (Carvalho, 2013;Rodrigues, 1987). Furthermore, it is a medium-size and generalist species found even in cities (Andrade, 2020;Carvalho, 2013;Vitt, 1995) but with a low dispersal rate (Vitt, 1995;Vitt et al., 1996). Given these traits, genetic structure would be expected to be absent or reduced and genetic variation should follow an isolation-by-distance model, as generalist species potentially are less affected by the patchiness of the landscape. In contrast, T. hispidus is distributed along a broad climatic gradient and heterogeneous landscape, suggesting that IBR and IBE could explain genetic variation and structure.
The genetic structure in Caatinga lizards has been related to vicariant barriers such as rivers, climatic gradients, differences in current and historical climate suitability, and geographical distances (Lanna et al., 2020;Oliveira et al., 2015;Werneck et al., 2015). Given its wide distribution in the Caatinga and generalist condition, Tropidurus hispidus is an excellent model to test these diversification processes upon the distribution of genealogical lineages and investigate the effects of landscape complexity on the genetic variation of Caatinga organisms. Here, we use a Caatinga-wide sampling of T. hispidus to assess three non-mutually exclusive hypotheses on the structuring of the species' genetic variation in the biome: (1) genetic differentiation is associated with spatial distance, (2) environmental gradients across the Caatinga are the primary drivers of genetic differentiation, and (3) landscape resistance predicts genetic divergence through physical barriers and climatically unsuitable areas. Hypotheses, their associated mechanisms, and expectations are summarized in Table 1.

Study area
The Caatinga is the largest continuous block of Seasonally Dry Tropical Forest in South America and is an exclusive biome from Brazil that extends for ca. 900,000 km 2 , equivalent to 11% of the Brazilian territory (Silva et al., 2017). The climate is semi-arid with high air temperature (annual average from 25 to 30 C), sparse and irregular rainfall, and seven to 10 months of an intense dry season (Silva et al., 2017;Prado, 2003). The semi-arid nature of this region and the low precipitation rates result from the low penetration of the moisture coming from the Atlantic Ocean, which reaches only the forest zone on the coast (Ab'Saber, 2003). The Caatinga topography, comprising interplateau depressions (300-500 m above sea level) and elevations up to 1,000 m (Ab 'Saber, 1974;Andrade-Lima, 1981), increases landscape complexity. There are few main perennial rivers in the Caatinga (see Fig. 1), and they have been intensely disturbed by human activities (da Silva et al., 2017).

Sampling and DNA extraction, amplification, and sequencing
We collected 59 specimens of Tropidurus hispidus from 27 localities throughout the Caatinga and adjacent areas in northeastern Brazil (Fig. 1). We complemented our dataset by downloading available sequences of T. hispidus from GenBank (http://www.ncbi.nlm.nih.gov/genbank), totaling 66 specimens from 31 localities (Fig. 1). These sampling points include areas from both margins of the São Francisco River and plateaus and depressions. Detailed information on geographic localities and voucher specimens is available in the Supplemental Material online.

Phylogenetic analyses
We downloaded sequences from other species of Tropiduridae from the GenBank to assess the monophyly of Tropidurus hispidus and check species identifications. We also used this dataset to estimate a timecalibrated gene tree. Detailed information on all specimens used in this analysis is available in the Supplemental Material Table S1.
First, we estimated the best substitution model and partition schemes based on the Bayesian Information Criterion using PartitionFinder2.1.1 (Lanfear et al., 2012), with a linked branch lengths model and a 'greedy' search. The best models were TrNEF þ I þ G for COI 1 st position, HKY þ I for COI 2nd position, and TrN þ G for COI 3rd position. We then estimated a calibrated, partitioned gene tree using the best substitution models, a birth-death tree model, and a lognormal uncorrelated relaxed clock as priors. The birth-death prior accounts for both speciation and extinction processes and is robust when including inter and intraspecific sampling (Ritchie et al., 2017). To avoid relevant impacts of the substitution rate heterogeneity among lineages, as our data also contains distant-related lineages, we assumed the lognormal uncorrelated relaxed clock (Sarver et al., 2019). We used a mean substitution rate of 1.07 Â 10 À2 site/million years (SD ¼ 0.002) available from Bernardo et al. (2019) and a normal prior distribution of divergence times to calibrate the tree based on the COI gene, following the callibration used for phrynosomatid lizards. Next, we implemented three independent runs in BEAST 1.10.4 (Suchard et al., 2018), each consisting of 5 Â 10 7 generations and sampled at every 5 Â 10 3 generations. We discarded each run's initial 25% generations as burn-in and combined the results with LogCombiner v1.10.4 (Suchard et al., 2018). Finally, we assessed the convergence of runs based on effective sample sizes (ESS ! 200) using Tracer v1.6 (Rambaut & Drummond, 2007).

Delimitation of mtDNA haplogroups and population genetics
We implemented a Bayesian Analysis of Population Structure (BAPS) to assess the presence and nature of spatial genetic structure in Tropidurus hispidus. Because IBD was present in our data (see results), we estimated the number of genetic clusters (k) by applying the spatially explicit BAPS model (see Perez et al., 2018), using the option 'spatial clustering of individuals'. We implemented this analysis in BAPS v6.0 (Corander et al., 2008) with a population mixture model and k ranging from 1 to 10. Then, we conducted population admixture analysis with the number of iterations to estimate the admixture coefficients set to 100, 200 reference individuals from each population, and 20 iterations per individual. We applied the tree-based species Table 1. Hypotheses of genetic differentiation in the generalist lizard Tropidurus hispidus within the Caatinga. Each mechanism is followed by a general pattern expected in the data.

Hypothesis
Mechanism Expectation Isolation by distance (IBD) Gene flow is more likely among closest populations because of dispersal constraints.  Kapli et al., 2017) to identify the number of independent evolving units in mtDNA data. Among the tree-based single-locus species delimitation methods, mPTP was shown to be more effective with unbalanced sampling (Blair & Bryson, 2017). Because this analysis uses the number of substitutions to estimate branching processes, and requires a phylogram, we transformed the ultrametric tree generated in BEAST (scaled by time) with the "phangorn" v.2.8.1 R package (Schliep, 2011). We used the resulting optimized phylogram as input on the mPTP website (https://mptp.h-its.org/). Because mtDNA haplogroups delimited by BAPS and mPTP were partially discordant (see Results), we made comparisons among groups considering the most conservative results from BAPS. We estimated the following statistics for each haplogroup identified by BAPS using DnaSP v5.10 (Librado & Rozas, 2009): number of haplotypes, number of polymorphic sites, haplotype diversity, nucleotide diversity, Tajima's D, Fu's F, and Ramos-Onsins and Rozas' R 2 (using 1,000 coalescent simulations). Next, we used a Hierarchical Analysis of Molecular Variance (AMOVA) to assess the departure of population genetic structure from panmictic expectations with Arlequin v.3.5.2 (Excoffier & Lischer, 2010). We also calculated average genetic distances (uncorrected p-distance) among haplogroups in MEGA X (Kumar et al., 2018). Finally, we built a haplotype network using the median-joining method with PopArt v1.7 (Leigh & Bryant, 2015).
To assess temporal demographic dynamics, we built Bayesian Skyline Plots (Drummond et al., 2005) for each haplogroup with more than ten samples with BEAST v1.10.4 (Suchard et al., 2018). For all haplogroups we used the HKY þ G substitution model. We used a substitution rate of 1.07 Â 10 À8 /site/year (SD ¼ 2 Â 10 À9 ) (Bernardo et al., 2019) with a normal prior distribution and a strict clock model. Each MCMC run was performed with 2 Â 10 7 generations, sampling every 2 Â 10 3 generations and a 10% burn-in. We assessed the convergence of runs and generated plots in Tracer v1.6 (Rambaut & Drummond, 2007). We did not build a skyline plot for the Southeast haplogroup because we had less than ten samples.

Environmental suitability models
We constructed ecological niche models (ENM) to predict environmentally suitable regions for Tropidurus hispidus using species occurrences and environmental data. ENMs results were used to analyze whether current or past climatic conditions explain patterns of genetic differentiation in T. hispidus (see below). First, we gathered T.
hispidus occurrence data in scientific collections and literature (Carvalho, 2013;Gorzula & Señaris, 1998), totaling 185 unique geographic records. We filtered records less than 10 km apart to reduce sampling bias using the "spThin" R package v0.2.0 (Aiello-Lammens et al., 2015), resulting in 175 records. Second, we downloaded environmental data at a spatial resolution of 2.5 arc-min from CHELSA (Karger et al., 2017). We opted for this resolution as it was the finest for some paleoclimate projections and took into account localities' potential uncertainties with approximate geo-referencing projections. We obtained 19 bioclimatic variables and removed those high collinear (Pearson's r! 0.8) while keeping biologically relevant bioclimatic variables, as follows: temperature seasonality (BIO4), mean temperature of the driest quarter (BIO9), precipitation of wettest month (BIO13), precipitation seasonality (BIO15), precipitation of driest quarter (BIO17), precipitation of warmest quarter (BIO18), and precipitation of coldest quarter (BIO19).
To build the ENMs, we used a broad spatial extent as background, including areas for the potential distribution of Tropidurus hispidus in South America. We used the MaxEnt algorithm (Phillips et al., 2006) to generate models and the "ENMeval" R package v2.0.3 (Muscarella et al., 2014) for tuning and evaluating models, with the following settings: 10,000 background points, checkerboard1 method, and L, H, LQ, LQH, LQHP, and LQHPT feature class combinations. We conducted model selection based on the Akaike information criterion corrected for small samples (AICc) and used the area under the curve (AUC) to assess model performance. AUC ranges from 0 to 1, in which 0.5 means random prediction, values between 0.80 and 0.90 suggest that model performance is good, while higher values are considered excellent. We also assessed model performance with the Boyce index (values range from À1 to 1) in the "ecospat" R package v3.2.1 (Di Cola et al., 2017), in which values close to zero are no better than a random model and values closer to one means better models (Hirzel et al., 2006). Based on the seven bioclimatic variables used in modeling, we then projected environmental suitability into four periods during the past 787,000 years: Marine Isotope Stage 19 in the Pleistocene (MIS19; 787 kya) (Brown et al., 2018), Last Interglacial (LIG; 120 kya) (Otto-Bliesner et al., 2006), Last Glacial Maximum (LGM; 21 kya) (Karger et al., 2017), and Mid-Holocene (Holocene; 8.326-4.2 kya) (Fordham et al., 2017). We downloaded the paleoclimate data from Paleoclim (Brown et al., 2018; http://www. paleoclim.org/).

Landscape heterogeneity and population connectivity
To evaluate the effects of landscape heterogeneity on gene flow and movement of populations, we assessed the relative contributions of IBD, IBE, and IBR on the genetic differentiation in Tropidurus hispidus. To this end, we tested nine potential promoters of genetic differentiation (more details below; see also Table 1), and we used 19 localities with at least two samples to build a genetic distance matrix with G ST -values (Nei, 1973), considering the averaging individuals within sites, as the response variable, using R package "mmod" v.1.3.3 (Winter, 2012). We chose G ST because it is a non-negative measure of genetic differentiation.
IBD predicts that genetic differentiation increases with geographical distance simply due to restricted gene flow between populations far apart. We used geographical coordinates of each locality directly in Generalized Dissimilarity Modelling (see below).
To test the scenario where environmental differences across the landscape can predict genetic differentiation (i.e., IBE), we used 19 present-day bioclimatic variables from CHELSA (see the above section). We extracted these variables for each locality and transformed them via principal component analysis with the 'prcomp' function in R. We retained the first three PCs, which explained 77%, 15%, and 5% of the variance, respectively, and calculated pairwise environmental (Euclidean) distances between populations for each PC, with the 'dist' function in R.
To assess the influence of resistance surfaces upon genetic differentiation, i.e., the IBR hypothesis, we used topography features (terrain slope and roughness-a measure of surface complexity), the presence of large rivers, and climatic suitability (current and LGM maps) as predictors. We derived slope and roughness rasters from elevation (NASA Jet Propulsion Laboratory; available at https://landscape.jpl.nasa.gov/). We downloaded a shapefile of major rivers (at 1:10 m scale) from Natural Earth Dataset (https://www.naturalearthdata.com) and converted it to ASCII format in GDAL/OGR. Rasters were at a 2.5 arc-min spatial resolution. We used circuit theory to generate resistance models for movement between sample sites and produced pairwise resistance matrices with Circuitscape 5 (Anantharaman et al., 2020). Based on habitat heterogeneity, Circuitscape calculates the dispersal cost between two points across a landscape, considering all possible paths between them (McRae, 2006;McRae & Beier, 2007). To use environmental suitability maps (current and LGM models) as resistance models, we inverted them such that higher values indicate environmental resistance. We also produced resistance models from rivers in which the main perennial rivers (grid cells set to 1) were considered barriers, implying higher dispersal costs. Similarly, we considered areas with steeper slopes and higher roughness as offering higher resistance to dispersal across the landscape. We replaced all zeros by 0.0001 because Circuitscape reads 0 as hard barriers. Therefore, we produced resistance matrices based on: (i) current climatic unsuitability, (ii) LGM climatic unsuitability, (iii) rivers, (iv) slope, and (v) roughness.
To investigate the effects of IBD, IBE, and IBR on the genetic differentiation of Tropidurus hispidus, we implemented Generalized Dissimilarity Modelling (GDM) with the "gdm" R package v1.5.0 (Manion et al., 2018). This analysis is a matrix regression procedure that accounts for non-linear rather than linear relationships using I-spline basis functions (Ferrier et al., 2007). Initially, we fitted the response variable (G ST matrix) to the nine predictor variables, as follows: pairwise matrices of geographical distances (Euclidean distances matrix), climate (PC clim1 , PC clim2 , and PC clim3 ), and resistance (current and LGM unsuitability, rivers, slope, and roughness). The relative importance of each predictor can be viewed through the maximum height of response curves while keeping all other variables constant (Fitzpatrick & Lisk, 2016). To estimate predictor importance and significance (P 0.05), we used the 'gdm.varImp' function of the "gdm" R package, with 1,000 permutations. This procedure consists of matrix permutation and backward elimination, in which one predictor is removed at a time, and then the importance and significance of the remaining predictors are reassessed. The backward elimination continues until only significant predictors with unique contributions remain in the final model. We conducted all analyses in the R environment v4.1.2 (R Development Core Team, 2020).

DNA polymorphism, delimitation of mtDNA haplogroups, and population genetics
We obtained a final alignment of 540 bp for 66 individuals of Tropidurus hispidus, and we found no stop codons. We found high levels of genetic diversity for T. hispidus, with 70 polymorphic sites, 46 haplotypes, and nucleotide diversity of 0.025 ( Table 2). The mPTP and BAPS analyses were overall congruent, considering significant splits. However, mPTP showed additional internal breaks, delimiting ten haplogroups, while BAPS recovered four (Fig. 2). The haplotype distributions of both analyses can be viewed in Fig. 1. The BAPS mtDNA haplogroups are arranged as follows ( Fig. 1): one distributed mainly from eastern Caatinga (hereafter East haplogroup), one predominantly from western Caatinga (West haplogroup), and two from southern Caatinga (South and Southeast haplogroups). The South and Southeast haplogroups occur on the right margin of the São Francisco River (SFR), while the West haplogroup is distributed mainly on the left margin of the SFR. The East and West haplogroups partially overlap. The South and Southeast haplogroups occur in sympatry in Morro do Chap eu, Bahia state (Chapada Diamantina); and probably in the ecotone region between Atlantic Forest and Caatinga in Bahia, and haplogroups East and West cooccurred in Monte Alegre, Sergipe state. Although the haplogroups generally had a small geographic overlap, the haplotype genealogy showed no shared haplotypes (Fig. 1B).
The four haplogroups showed high levels of genetic diversity, with haplotype diversity varying from 0.89 to 0.97 (Table 2). Population size change (R 2 ) tests had significant values for all haplogroups, while significant values for the neutrality (Fu's Fs) were detected only for South and East ( Table 2). The hierarchical AMOVAs revealed that most genetic variation (68%, P < 0.0001) is among BAPS haplogroups. The uncorrected p-distance among haplogroups showed genetic differences of approximately 3.0%, as follows the pairwise comparisons: 3.5% between South and East, 3.6% for South/Southeast, 3.1% for South/West, 2.1% for East/Southeast, 3.1% for East/West, and 2.8% between Southeast and West.

Phylogenetic analysis and dating
One sample of Tropidurus torquatus (GenBank number KM588063) grouped with T. hispidus (Fig. 2) and, therefore, we could not ensure the monophyly of T. hispidus based on the COI Bayesian tree. However, this might result from misidentifying the sample because it is from a contact zone between the two species in southeastern Bahia (h11 from Fig. 1) and presents a low genetic distance to some T. hispidus samples. In such a case, the monophyly of T. hispidus based on the COI Bayesian gene tree is assured.

Environmental suitability models
The ENMs had high accuracy, with an AUC value of 0.85 and Boyce index of 0.98 (Supplemental Material  Table S2. However, past models over-predicted areas for southeast and central Brazil, Paraguay, Bolivia, and Peru, in which there are no records for Tropidurus hispidus. The models projected suitable climatic areas connecting the Caatinga and reaching Venezuela, Guiana, and northern Brazil, via a narrow corridor (Fig. 4). This corridor vanished during specific periods, such as the Holocene and LIG. During the past 790 kya, northeastern Brazil had large and interconnected climatically suitable areas for T. hispidus, with a more pronounced shrinkage in LGM. Our projections for the LIG and MIS19 indicated suitable areas in the northern Atlantic Forest, suggesting that hotter and wetter periods favor range expansion in T. hispidus (Fig. 4).

Landscape heterogeneity and population connectivity
The GDM indicated that isolation by distance (IBD) and landscape resistance (IBR) predicted most genetic differentiation in Tropidurus hispidus. Five predictors accounted for 41.24% of the null model deviance, presenting I-splines above zero: geographic distances (I- Fig. 2. Tropidurus hispidus maximum clade credibility Bayesian mtDNA tree (COI) highlighting the divergence time of haplogroups (below the tree topology) and posterior probabilities above 0.7 (in each node). Bars represent the 95% Highest Posterior Density (HPD) interval for divergence dates. Haplotype numbers are also given in the terminals. Different colours in mtDNA haplogroups are based on Bayesian Analysis of Population Structure (BAPS) and mPTP delimitations (see also Fig. 1). spline ¼ 0.79), resistance through LGM unsuitability (Ispline ¼ 0.68), resistance through rivers (I-spline ¼ 0.69), resistance through differences in roughness (Ispline ¼ 0.55), and environmental differences in PC clim1 and PC clim3 (I-spline ¼ 0.08 and 0.04, respectively) (Fig.  5). However, only resistance through LGM unsuitability and geographic distance contributed significantly (P 0.05), with a joint contribution of 37.89% of the genetic variation. In this case, genetic divergence increased with the geographic distance and higher unsuitability surfaces during the LGM period. The variance partitioning analysis showed that geographical distance and resistance through LGM unsuitability accounted for 76% and 24% of the explained variation, respectively.

Discussion
The generalist lizard Tropidurus hispidus is characterized by high genetic diversity and moderate geographic structure. We found that isolation by landscape Fig. 3. Bayesian skyline plots for each mtDNA haplogroup of the generalist lizard Tropidurus hispidus. The x-axis represents years from present to past, while the y-axis represents the effective population size (Ne). The blue line indicates the median of Ne, while the black lines indicate a 95% higher posterior probability. resistance through historical climatic unsuitability and spatial distances explained much of the genetic variability, with a more significant influence of the latter. Physical barriers, such as complex topography and rivers, had no significant effect on structuring the species' genetic variation. These results seem to fit the expectations for a widely distributed and generalist lizard species, in which IBD is the predominant factor explaining intraspecific genetic differentiation. Furthermore, the relationship between genetic differentiation and LGM unsuitability suggests that past climates could have mediated gene flow among populations by decreasing connectivity through unsuitable climatic areas. Still, we found highly divergent haplotypes in higher elevation regions, such as Morro do Chap eu, Bahia (Chapada Diamantina), highlighting topographic complexity in driving genetic differentiation (Rodr ıguez et al., 2015).

Effects of geography, environment differences, and landscape resistance
Dispersal ability and spatial/environmental features of the landscape can mediate population connectivity. Only IBR and IBD significantly explained the genetic divergence between sampled demes. Regardless of the potential effects of environmental gradients and landscape complexity, isolation by distance likely reinforced genetic differentiation among haplogroups by favouring gene flow among closer populations. This finding agrees with previous studies on fauna diversification in the Caatinga, in which IBD explained part of the differentiation among several taxonomic groups (Magalhaes et al., 2014;. Conversely, landscape resistance can be more critical to explaining gene flow (Jaff e et al., 2019;. Indeed, we recovered climatically unsuitable areas during the LGM as crucial barriers to gene flow in T. hispidus. Rivers and terrain roughness, conversely, had no significant effects on genetic differentiation, likely reflecting the generalist nature of the species and its ubiquity in the Caatinga, from pristine to disturbed habitats. Thus, evolutionary studies should consider natural history aspects, as they can influence how barriers and environmental gradients affect gene flow (Zamudio et al., 2016).
In a diversification scenario with no evident current or past geographic barriers between lineages, ecological features along environmental gradients (e.g., climate, habitat type, elevation, seasonality) may decrease gene flow, promoting differentiation among lineages, and ultimately lead to speciation (Endler, 1973;Nosil, 2008Nosil, , 2012. Intuitively, we expect that gene flow in species widely distributed to be mediated by the discontinuities or connections across environmental gradients. Conversely, Ito et al. (2020) detected no genetic structure in Tropidurus hispidus suggesting high connectivity among demes at least at small geographic scales. Our results based on a single-locus approach, but covering a broader spatial distribution in the Caatinga, indicate a moderate geographic structure for the species in this region. However, the strong genetic structure, inferred by the lack of shared haplotypes among mtDNA clusters and the presence of IBD suggests that other scenarios not evaluated in our study can explain such inconsistencies between geography and genetic structures. The first one relies on the haploid and usually matrilineal inheritance of the mtDNA data. Given the smaller N e , high mutation rates and faster coalescence time than to single nuclear genes, mtDNA has proven to effectively estimate the divergence at an inter-populational level (Hung et al., 2016). However, mito-nuclear discordances can emerge, among other factors, like sex-biased dispersion (Latta, 2006), and the use of mtDNA alone can limit our inferences of the genetic/geographic structure patterns. Males of T. hispidus seem to be highly dispersive as shown by bi-parental markers at a finer geographical scale (Ito et al., 2020); thus, we cannot rule out the effects of solely using mtDNA data on the phylogeographical structure we found.
Although this first scenario can partially explain the genetic structure found, it does not clarify how haplotypes of different clusters are present in nearby localities, or even distant from their core distribution (see cooccurrence of haplotypes h42 and h29 in Fig. 1). The second and not exclusive scenario explaining this pattern is human-mediated dispersion. Tropidurus hispidus occurs in high densities in urban environments and can be easily found on concrete walls and other artificial structures used as shelter or egg-laying sites (Andrade, 2020). In northeastern Brazil, Cear a and Bahia states are responsible for most of the production of ceramic materials sold to other states in the region (Mello et al., 2011). The commerce of these materials may promote lizard dispersal. Besides, a denser sampling could explore a finer resolution of the clusters' boundaries, as several unsampled or extinct haplotypes were inferred in the haplotype network. Further studies using diffusion/ dispersal analyses coupled with more variable nuclear markers such as SNPs could help track these dispersal events more precisely. The combination of IBD and resistance due to historical climatic unsuitability LGM promoted a substantial genetic differentiation ($38%) in Tropidurus hispidus. However, the use of mtDNA markers in studies at fine spatial and temporal scales must be weighted (Anderson et al., 2010). Compared to other more variable markers, as stated above, mtDNA can produce delayed responses to fast landscape changes (Anderson et al., 2010;Hall & Beissinger, 2014). As LGM unsuitability was the most critical predictor of genetic distances, it is reasonable to be cautious since the current genetic differentiation may result from the mtDNA's slower mutation rate. For example, genetic diversity in the co-distributed lizard Ameivula ocellifera was positively associated with LGM suitability, while genetic differentiation was mainly predicted by rivers resistance . The role of microevolutionary processes in the diversification of Caatinga vertebrates is still underappreciated. To our knowledge, this is the second study that considers the putative effects of landscape features on the lizard's genetic differentiation in this region.

Spatial genetic structure of Tropidurus hispidus and the Late Pleistocene climatic shifts in Caatinga
During the Pleistocene (2.58 Mya À11,700 Mya), significant paleoenvironmental changes impacted the Caatinga landscape  and were crucial for the intraspecific genetic divergence of species (Lanna et al., 2020;Werneck et al., 2012;Werneck et al., 2015). Geomorphological events in this period changed the course of the SFR, and in the wettest periods, rainforests expanded through the Caatinga, leaving genetic signatures Carnaval & Bates, 2007;Thom e et al., 2016;Werneck et al., 2015). Most intraspecific divergence of Tropidurus hispidus happened in the Quaternary (ca 2 Mya), suggesting that environmental changes in this period drove the structuring of the species' genetic variation. Likewise, the lineages of the T. semitaeniatus species group, a saxicolous congener endemic to the Caatinga and sympatric with T. hispidus, also diverged in the same period (Werneck et al., 2015).
The impact of recurrent climatic changes during Pleistocene on the "diagonal" of open formation's biota has been studied with environmental niche modeling and molecular methods that try to reconstruct the diversification history of lineages through time (Collevatti et al., 2015;Costa et al., 2018;Prado et al., 2012;Werneck, Costa, Colli, Prado, & Sites Jr, Werneck et al., Werneck et al., 2011). Predicted biomes' distribution models for South America suggest a contraction and fragmentation of savannas and Seasonally Dry Tropical Forests during the LGM Silveira et al., 2019), when many species with different life-history traits experienced a demographic expansion in the Caatinga (Gehara et al., 2017). We expected that populations of Tropidurus hispidus were subject to varying degrees of demographic reduction during periods of rainforest expansion, while recolonizations and demographic expansions would follow in hotter/drier periods. Our results partially corroborate these expectations. Indeed, paleoclimatic models suggest large suitability areas along the Late Pleistocene, comprising most of the current distribution in Caatinga. Besides that, different climatic stability areas in the past 790 thousand years may have kept viable populations, reducing the impacts of extinctions during the expansion of rain forests. During the expansion of open and dry environments, gene flow among populations of T. hispidus from Caatinga with those from westernmost distributions (e.g., Amazonia and dry habitats of Venezuela) may have been possible, mainly during LGM (Fig. 4). Conversely, demographic history varied among haplogroups (Fig. 3). While the South haplogroup showed signals of demographic expansion after the LGM, the East showed signs of an older population expansion during LIG, while the West haplogroup remained stable in the Late Pleistocene. Previous studies on the herpetofauna of the Caatinga region suggest that most anurans, lizards, and snakes show synchronous population expansions in the Late Pleistocene (Gehara et al., 2017).
Many questions on the patterns and processes of diversification in open arid and semi-arid environments remain open. Lizards are highly diverse in arid and semi-arid environments (Pie et al., 2017;Powney et al., 2010;Wiens et al., 2013). Given their biology (e.g., ectothermic and low water dependence), tropidurid lizards are a promising group for testing the main drivers of diversity in these environments. Because Tropidurus is a genus incredibly diverse in open and dry environments of South America, the hypothesis of niche conservatism should be considered when testing the leading promoters of speciation in the group. Studies using several loci and a representative geographic sampling of species, testing explicit biogeographic hypotheses, should elucidate this group's tempo and mode of diversification. Given the large contact zones, the number of evolutionary lineages, and speciation processes inferring whether the T. torquatus group can be a case of radiation in the open environments of South America and if hybridization plays an essential role in the evolutionary history of these animals. Future studies using morphological and genetic data should clarify the taxonomic status between T. torquatus species group populations in contact zones and whether these species can hybridize.
Our results are a first step to shedding light on the main drivers of the evolutionary history of Tropidurus hispidus. We showed that isolation by resistance and isolation by distance shaped the genetic structure of this generalist lizard. Climatically unsuitable areas during the LGM reduced gene flow among populations, being a significant predictor of the species' genetic variation. Thus, territorial behaviour and climatic differences along the Caatinga can partially explain the genetic structure observed in T. hispidus. However, besides a denser sampling along the Caatinga and other environments, given the limitations associated with using only one marker, future studies should focus on using more variable markers such as SNPs and microsatellites. This should allow assessing evolutionary scenarios accounting for mito-nuclear discordance and intersexual differences in dispersion rates. We are starting to understand how the Caatinga's landscape influenced the genetic differentiation in lizards. Future evolutionary studies using climate-associated candidate genes are necessary to evaluate whether local adaptation occurs in this species and test possible connecting routes and divergence times of Amazon populations.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Supplemental material
Supplemental material for this article may be accessed here: https://doi.org/10.1080/14772000.2022.2084470.