Niche Divergence at Intraspecific Level in the Hyrcanian Wood Frog, Rana pseudodalmatina: A Phylogenetic, Climatic, and Environmental Survey

The role of ecological niche divergence in lineage speciation has recently stimulated the interest of evolutionary biologists and ecologists. Phylogenetic analysis has revealed that the Hyrcanian wood frog, Rana pseudodalmatina, has diverged into two western and eastern regional clades (WRC and ERC) within the Hyrcanian forest. The goal of this study was to investigate whether the ecological niches of WRC and ERC are conserved or diverged, as well as to figure out what variables promote niche conservatism or divergence. For this purpose, the maximum entropy model was employed to assess environmental niche modeling in geographical (G) space utilizing climatic and macro-environmental data. The niche overlap, equivalency, and similarity tests based on PCAenv analyses were used to assess niche divergence or conservatism in environmental (E) space. The findings strongly support the hypothesis that WRC and ERC have undergone substantial niche divergence and are constrained by a unique set of climatic and macro-environmental conditions. This study by ecological niche comparisons based on phylogenetic data provides new insights into the exploration of species diversification processes in the Hyrcanian forests.

The evidence supporting niche divergence or conservatism appears to be conflicting at first glance (Hu et al., 2015). However, there is a considerable structure when the patterns are time-structured. The temporal frames of the study systems and niche tests (i.e. niche similarity vs niche equivalency tests) influence results of niche conservatism or divergence (Peterson, 2011;Broennimann et al., 2012). The equivalency tests are more important to determine whether niche models are transferable across time and space, whereas niche similarity tests are more appropriate to test evolutionary and biogeographic assumptions (Peterson, 2011). Niches generally tend to conserve over time in terms of distribution patterns and speciation, whereas niche divergence considers being the driving force behind species diversity along ecological gradients (Graham et al., 2004;Raxworthy et al., 2007). Short-term and recent events, such as distributional shifts over short periods or species invasions, demonstrate significant tendencies toward conservatism. Longerterm events, such as phylogenetic differentiation, however, demonstrate a tendency for conservatism to break down (Wiens and Graham, 2005;Peterson, 2011). According to evidence for niche divergence, divergent natural selection promotes diversity by allowing individuals to adapt to new environments (Schluter, 2009).
The Hyrcanian wood (or brown) frogs, Rana pseudodalmatina (Anura, Ranidae), have been recorded from southeastern Azerbaijan to the Golestan National Park in northeastern Iran (Najibzadeh et al., 2018;Kidov and Litvinchuk, 2021). This study focuses on Iranian frog populations in the Hyrcanian (or Caspian, temperate) forests at altitudes ranging from −22 to 2850 meters above sea level along the northern slopes of the Alborz Mountains and the southern edge of the Caspian Sea. Based on two partial mitochondrial markers (cytochrome b and 16S rRNA), Najibzadeh et al. (2018) clearly showed the existence of two clades in the west and east of the Hyrcanian forest for this species (Najibzadeh et al., 2018). In another study using similar markers, Amiri et al. (2021) demonstrated that there existed two regional patterns within the species distribution range in the west and east, which diverged in the Pleistocene (1.6 Mya) (Amiri et al., 2021). Using species distribution modeling (SDMs hereafter), this study proposed the hypothesis of "refugia-withinrefugia" in this area (Amiri et al., 2021). Based on biogeographic analyses, they suggested that dispersal and vicariance (rise of Caspian Sea water levels) processes may have influenced this species' genetic structure (Amiri et al., 2021).
For this purpose, the goal of this study was to investigate whether the R. pseudodalmatina clades occupied the equivalent or similar niche, and if not, what role niche divergence or conservatism had in their evolution. Specifically, we used a maximum entropy approach (MaxEnt) to create ENMs for each clade based on large-scale environmental factors and occurrence data. Following that, using environmental and occurrence data, we investigated whether these two clades occupy more similar environments than estimated based on background (Hyrcanian regions) environmental differentiation using niche overlap, niche equivalency, and niche similarity tests based on principal component analysis (PCA env ). The results of these studies, when combined with the phylogeographic relationships (Najibzadeh et al., 2018(Najibzadeh et al., , 2021Amiri et al., 2021), give a comprehensive and multidimensional perspective of niche variation and differentiation in the Hyrcanian wood frog.  Table 1). We excluded the Azerbaijan data from our analysis due to a lack of genetic information to assess its occurrence in a particular clade in this distribution region. Therefore, we classified occurrences into two western and eastern regional clades (WRC and ERC hereafter) based on the previous phylogenetic analyses (Figure 1; Najibzadeh et al., 2018Najibzadeh et al., , 2021Amiri et al., 2021). We evaluated all occurrences equally, regardless of population size, and double-checked them using spreadsheets and GIS to search for duplication and possible georeferencing errors. In parallel, we checked occurrences spatially at a spatial resolution of 30 arcseconds to ensure that only one record per grid cell remained for each clade. The final dataset included 71 georeferenced occurrences, with 33 for WRC and 38 for ERC (Figure 1 and Supplementary

Explanatory Variables
We gathered 32 environmental variables to quantify environmental heterogeneity within the distribution range of Hyrcanian wood frogs. These comprised 19 bio-climatic FIGURE 1 | Study area. Hyrcanian forests along the northern slopes of the Alborz mountain range and the southern coastline of the Caspian Sea. The green circles represent the western regional clades (WRC) and the red circles represent the eastern regional clades (ERC) of the Hyrcanian wood frog, Rana pseudodalmatina within its range. A topographic map serves as the background.
(BCV hereafter) and 13 macro-environmental (MEV hereafter) variables, including climate, soil, hydrology, land cover, topography, sunlight, and human impact, all of which identified as significant factors potentially affecting species distribution limits (Supplementary Table 2; Hu et al., 2016). Because high colinearity between environmental variables may inflate SDMs' model accuracy, dimension-reduction methods (e.g., clustering algorithms and/or correlation analysis) can be employed to minimize correlations among variables (Hu and Jiang, 2010;Synes and Osborne, 2011). Therefore, we reduced the number of predictor variables by combining the results of Pearson's correlation tests (certain BCV variables eliminated due to high correlations with other BCV with | r| > 0.75) (Goudarzi et al., 2019;Vaissi, 2021a,b) and a jackknife analysis (retaining MEV with the higher value when used in isolation) (Phillips et al., 2006;Hu et al., 2016 (Karger et al., 2017a,b); annual actual evapotranspiration (AET hereafter), annual aridity index (AI hereafter), annual potential evapotranspiration (PET hereafter) from the consortium for spatial information 3 ; land-cover (LandCov hereafter) from the Google Earth Engine (GEE) 4 (Ghorbanian et al., 2020); altitude (Alt hereafter) from the worldclim database 5 and slope and aspect, which produced from the digital elevation model using the surface tools in ArcMap (v 10.8) (Supplementary Table 2). The spatial resolution of all variables was 30 arc-seconds (∼1 km). Variables measured in this study reflected environmental factors that frogs encountered and that affect amphibian physiology and survival (Araújo et al., 2006;Buckley and Jetz, 2007;Wells, 2010;Hof et al., 2011;Hu et al., 2016).

Environmental Niche Modeling
To create ENMs, we used the maximum entropy model (MaxEnt v 3.4.4; 6 a presence-background approach that has a high prediction power even with a small sample of presence sites (Phillips et al., 2006) and widely used in spatial ecology (Banerjee et al., 2019;Cardoza-Martínez et al., 2019;Martínez-Freiría et al., 2020;Rodríguez-Rodríguez et al., 2020;Goudarzi et al., 2021). We constructed nine distinct models by MaxEnt using (1) BCV at species-level (SL hereafter) and two WRC and ERC, (2) MEV at SL and two WRC and ERC, and (3) integrated BCV + MEV at SL and two WRC and ERC. According to recent studies, using the default setting is not always the best option, particularly when the sample size is small (Merow et al., 2013;Radosavljevic and Anderson, 2014;Morales et al., 2017). Furthermore, the number of model parameters has an impact on model complexity, which can lead to overfitting (Shcheglovitova and Anderson, 2013). Therefore, we calculated the beta multiplier and reduced the number of included features and variables to reduce model complexity and over-parameterization (see more details in the section of explanatory variables) (Moreno-Amat et al., 2015). We evaluated the most parsimonious WRC, ERC, and SL models using different regularization beta-multiplier values (0.5, 1, and 2) and Akaike's Information Criterion compensated for a small sample size (AICc). As a result, the final models used to have a beta multiplier of 0.5. Furthermore, we used the bootstrap method to generate 30 replicates, with 75% of the occurrences for model training and 25% for testing (Goudarzi et al., 2021). We used the percentage of contribution and jackknife analyses to determine the importance of each variable. We calculated the average value of replicates using the AUC metric (Area under the ROC curve), which is a threshold-independent measure of discriminating capacity (Swets, 1988). We chose the logistic output format for displaying the anticipated environmental suitability because of its simplicity of comprehension, with values ranging from 0 (lowest) to 1 (highest) (Phillips and Dudík, 2008).

Ecological Niche Metrics
We compared and visualized the niches between the WRC and ERC ranges using the niche overlap test, principal component analysis (PCA env hereafter), niche equivalency (or identity) test, and niche similarity (or background) test from the eight BCV and seven MEV. Niche overlap is the intersection of two niches in a niche space. We used Schoener's D-metric (Schoener, 1970;Warren et al., 2008), which ranges from 0 (no overlap) to 1 (complete overlap), to calculate the niche overlap between WRC and ERC. We utilized the PCA env method to assess equivalency and similarity between the realized niches of WRC and ERC. This method estimates the available environmental space defined by the first two axes of the PCA by comparing the environmental circumstances available for a species within a specific study extent (background) with its observed occurrences. A smooth kernel density function is used to adjust for sampling bias in this method .
The niche equivalency and niche similarity tests frequently use the niche conservatism or divergence hypothesis (Glennon et al., 2014). We utilized the niche equivalency based on Brown and Carnaval's (2019) ordination technique in environmentalspace (E-space hereafter), which is the outcome of extending the methods of Broennimann et al. (2012), Petitpierre et al. (2012), and Qiao et al. (2017). We also used Brown and Carnaval's (2019) method of expanding Beale et al. (2008), Warren et al. (2010), and Nunes and Pearson (2017) techniques to measure niche similarity tests. The niche equivalency test is conservative since it only considers whether the two entities' niches (here ERC and ERC) are equal in their niche spaces, considering only the occurrences and leaving the backdrop out (Warren et al., 2008;Aguirre-Gutiérrez et al., 2015). The niche similarity test determines whether one range's environmental niche is more similar to the other range's niche than predicted by chance . The null hypothesis is that specific niches are equivalent when comparing a non-significant niche equivalency test with a significant niche similarity test. The null hypothesis of niche equivalency is rejected when a niche equivalency test is statistically significant, regardless of the significance (or non-significance) of the niche similarity test (Brown and Carnaval, 2019).
We used the ecospat package in R 4.1.1 7 to perform niche overlap test analyses (Di Cola et al., 2017). The Humboldt package was to perform PCA env , niche equivalency test, and niche similarity test analyses (Brown and Carnaval, 2019).

Environmental Niche Modeling
The ENMs displayed high predictive accuracy as measured by the AUC metric for the WRC (0.97 ± 0.008 based on BCV; 0.95 ± 0.02 based on MEV; 0.95 ± 0.03 based on BCV + MEV), ERC (0.97 ± 0.008 based on BCV; 0.97 ± 0.006 based on MEV; 0.98 ± 0.005 based on BCV + MEV), as well as at SL (0.94 ± 0.008 based on BCV; 0.96 ± 0.007 based on MEV; 0.97 ± 0.005 based on BCV + MEV) ( Supplementary  Figures 1-3). Potentially suitable habitats based on BCV, MEV, and integrated BCV + MEV for WRC and ERC, as well as at SL, are illustrated in Figures 2-4I-K. The WRC and ERC clearly illustrate distinct environmental suitability regions. Based on BCV, MEV, and BCV + MEV maps, it appears that for the western population, the west of the distribution range has high suitability compared to the east (Figures 2-4I). Similar results also appear for dispersed populations in the east (Figures 2-4J).
The percent contribution and permutation importance of BCV, MEV, and BCV + MEV for the WRC, ERC, and SL are given in Table 1. In the BCV Jackknife evaluation, the major climatic factors influencing WRC, ERC, and SL distribution were found to be Bio18, Bio18, and Bio17, respectively. According to the MEV jackknife test for WRC and ERC and at SL, the major macroenvironmental factors affecting the distribution of Hyrcanian wood frogs were LandCov, Alt, and LanCov, respectively. Based on the integrated BCV + MEV jackknife test for WRC and ERC and at SL, the environmental variables with the highest gain when utilized in isolation were Bio17, Bio18, and Bio17, respectively (Supplementary Figures 4-6).
According to Bios 12, 13, 14, 17, and 18, the greatest probability of the presence of Hyrcanian wood frogs in the west is in the precipitation of about 1100 (mm/year), 170 (mm/month), 50 (mm/quarter), 150 (mm/quarter), and 180 (mm/quarter), respectively. While in the east of the distribution, due to the same variables, the most likely presence of the species is in the lower precipitation of about 600 (mm/year), 80 (mm/month), 30 (mm/quarter), 70 (mm/quarter), and 110 (mm/quarter), respectively. Although there was no significant difference in the highest probability of species presence between WRC and ERC due to temperature (Bios 1, 2, and 8), the curves show that western populations are more likely to be present at lower temperatures than eastern populations (Supplementary Figures 7, 10).
Based on MEV, the probability of species presence increases with increasing AET, AI, and slope for populations in the West. While for the population distributed in the east, although the probability of the species' presence first increases with increasing AET, AI, and slope, but then decreases. In the west, increasing PET decreases the possibility of a species' presence, but in the east, it first increases and then decreases. According to Alt curves, western populations are more likely to be found at higher elevations than eastern populations (Supplementary Figures 8, 11).

Ecological Niche Metrics
The species niche of Hyrcanian wood frogs in the western and eastern distribution ranges along the two first axes of the PCA env based on BCV, MEV, and BCV + MEV are illustrated in Figures 2-4A,B. Figures 2-4C,D indicate the niche overlaps between WRC and ERC based on BCV, MEV, and BCV + MEV, which have little overlap. The first two components of PCA env based on BCV explained 88.05% of the total variation along climatic gradients (PC1 = 71.22%, PC2 = 16.83%) (Figure 2H). Based on MEV, the first two components of PCA env explained 70.64% of the total variation along micro-environmental gradients (PC1 = 51.25%, PC2 = 19.39%) (Figure 3H). According to both climatic and macro-environmental scales, PCA env indicated that 75.81% of the environmental variation was explained by the first two PCs (PC1 = 57.29%, PC2 = 18.52%) (Figure 4H).
The occupied niches by the two WRC and ERC were significantly non-equivalent based on BCV (p = 0.01), MEV (p = 0.03), and BCV + MEV (p = 0.009) (Figures 2-4E). According to BCV, the realized niches of WRC and ERC had a similarity test value of D = 0.03 and were not significant when comparing WRC to ERC (p = 0.18) ( Figure 2F) and ERC to WRC (p = 0.76) (Figure 2G). Based on MEV, the realized niches of WRC and ERC had a similarity test value of D = 0.03 and were not significant when WRC was compared to ERC (p = 0.94) (Figure 3F), or when ERC was compared to WRC (p = 0.76) ( Figure 3G). Based on BCV + MEV, the realized niches of WRC and ERC had a similarity test value of D = 0.002 and were not significant when WRC was compared to ERC (p = 0.28) ( Figure 4F) or ERC was compared to WRC (p = 0.84) ( Figure 4G).

DISCUSSION
Several studies have focused on the genetic structure, phylogeny, and habitat suitability of the Hyrcanian wood frog, R. pseudodalmatina (Najibzadeh et al., 2017(Najibzadeh et al., , 2018(Najibzadeh et al., , 2021Amiri et al., 2021;Kidov and Litvinchuk, 2021). However, the present study, using ENMs and multivariate niche analyses, for the first time, investigated environmental differentiation between the WRC and ERC of the Hyrcanian wood frog in the Hyrcanian forests. The hypothesis focused on whether R. pseudodalmatina lineages preserved similar niches after divergent, which might reveal new details about Hyrcanian wood frog diversification. For this purpose, niche divergence or conservatism and its relationship to lineage diversification were explored based on recent phylogeographic understanding (Najibzadeh et al., 2018(Najibzadeh et al., , 2021Amiri et al., 2021). The results of this study found evidence of niche divergences between WRC and ERC that may have accelerated the evolution of clades of R. pseudodalmatina. On the other hand, this research by ecological niche comparisons may provide new insights into the exploration of species diversification processes in the Hyrcanian forests. Furthermore, understanding the major ecological constraints to species distribution is critical for current conservation efforts as well as future research into the effects of climate change on biodiversity (Aguirre-Gutiérrez et al., 2015).
The factors that impact species distributions must be addressed and recognized when evaluating differences in species' environmental niches (Guisan and Thuiller, 2005;Colwell and Rangel, 2009). Therefore, this study focused on the gradients of abiotic factors such as climate, soil, hydrology, land cover, topography, sunlight, and human impact (Supplementary Table 2) that can constrain species ranges (Gaston, 2003). In Frontiers in Ecology and Evolution | www.frontiersin.org the Hyrcanian regions, the annual precipitation decreases from west to east as well as from lowlands to highlands (Tohidifar et al., 2016). The mean annual temperature also decreases from the lowlands to the mountain ranges. According to the CHELSA database, 8 annual precipitation in the Hyrcanian area is around 1430 mm in the west (Gilan province) and around 8 http://chelsa-climate.org 280 mm in the east (north of Khorasan province) (Karger et al., 2017a,b). The mean annual temperature in the lowland belt is around 19.2 • C (at an altitude of about −20 m a.s.l. in Gilan and Mazandaran provinces), whereas it is about 5.2 • C in the upper-montane belt (at an altitude of around 2,500 m a.s.l. in Mazandaran province). According to climate station databases in lowland cities, an increase in the mean annual temperature and a decrease in annual precipitation have been recorded from FIGURE 3 | The macro-environmental niche of the Hyrcanian wood frog, Rana pseudodalmatina across Hyrcanian forests, Iran. Panels (A,B) illustrate the species' niche in the western and eastern distribution ranges, respectively, along the two first axes of the PCA. The density of the species' occurrences in the cell is shown in the gray shade. The solid and dashed contour lines represent 100 and 50% of the available (background) environment, respectively. Panels (C,D) represent the niche overlap values for western (WRC) and eastern regional clades (ERC), respectively. Histograms (E-G) show the observed niche overlap D between the two ranges (bars with a diamond) and simulated niche overlaps (gray bars) on which tests of niche equivalency (E), niche similarity of WRC to ERC (F), and niche similarity of ERC to WRC (A). (H) The contribution of the macro-environmental variables on the two axes of the PCA env and the percentage of inertia explained by the two axes. Potentially suitable areas are based on macro-environmental variables according to WRC (I), ERC (J), and species-level (K).
west to east (Gholizadeh et al., 2020). Furthermore, according to a new bioclimatic classification, the climate in the eastern part of these forests differs from that in the western part and is classified as Mediterranean pluviseasonal oceanic and temperate oceanic, respectively (Djamali et al., 2011). The vegetation in these forests is diverse, as are the geological substrates, which include acidic rock (granite) in the west, dolomite in the center, and sandstone in the east (Siadati et al., 2010;Gholizadeh et al., 2020). Brown soil covers almost 90% of these forests, which can be followed by alluvial, colluvial, rendzina, and ranker soils (Sagheb-Talebi et al., 2014).
According to the findings, climatic and macro-environmental factors have a substantial impact on the distribution of the Hyrcanian wood frog. In more detail, based on percent FIGURE 4 | The integrated climatic and macro-environmental niche of the Hyrcanian wood frog, Rana pseudodalmatina across Hyrcanian forests, Iran. Panels (A,B) illustrate the species' niche in the western and eastern distribution ranges, respectively, along the two first axes of the PCA. The density of the species' occurrences in the cell is shown in the gray shade. The solid and dashed contour lines represent 100 and 50% of the available (background) environment, respectively. Panels (C,D) represent the niche overlap values for western (WRC) and eastern regional clades (ERC), respectively. Histograms (E-G) show the observed niche overlap D between the two ranges (bars with a diamond) and simulated niche overlaps (gray bars) on which tests of niche equivalency (E), niche similarity of WRC to ERC (F), and niche similarity of ERC to WRC (A). (H) The contribution of the bio-climatic and macro-environmental variables on the two axes of the PCA env and the percentage of inertia explained by the two axes. Potentially suitable areas are based on integrated bio-climatic and macro-environmental variables according to WRC (I), ERC (J), and species-level (K).
contribution and Jackknife analysis ( Table 1 and Supplementary  Figures 4-6), climatic variables, especially precipitation, have a greater impact on the distribution range of western and eastern populations than macro-environmental variables when both BCV and MEV are taken into account. As a result, Bio17 and then Bio18 and Bio14 are the most significant variables influencing the distribution of western populations, whereas Bio18 and then Bio12 and Alt are the most important variables affecting the eastern populations' range (Supplementary Figure 6). Based on contribution percentages, these factors were Bio18 and then 1 | Percent contribution (PC) and permutation importance (PI) of eight bio-climatic (BCV) and seven macro-environmental (MEV) variables for the two western and eastern regional clades (WRC and ERC) and at species-level (SL).

Variables
Based LandCov and slope for western populations, and Bio17 and then Bio1 and slope for eastern populations ( Table 1, the section of BCV + MEV). However, it should be noted that biotic interactions such as predators, competitors, and parasites all have an important role in restricting species distribution, which needs attention (Gaston, 2003;Gotelli et al., 2010;Sunday et al., 2011). Differences in distribution and environmental suitability between the two lineages might be explained by niche differences and divergent reactions to the environment (Peñalver-Alcázar et al., 2021). Evaluating whether ecological features of one lineage can effectively forecast the geographic distribution of another lineage (or itself), as well as vice versa, is a reliable way to investigate this prediction (Peterson et al., 1999;Peterson and Holt, 2003). In this study, ENMs were carried out separately for WRC and ERC as well as at the SL. The habitat suitability predictions for each population were significantly different between the WRC and the ERC. The results clearly reveal that the east of the distribution region is unsuitable for the western populations in terms of BCV, MEV, and BCV + MEV. Eastern populations exhibit similar results, with the exception that in the west, restricted distribution areas for these populations have been found, although they are not highly suitable (Figures 2-4). The results of the probability of species presence show that western populations have evolved in wetter and low-temperature environments that have high precipitation, AET, AI, Alt, and slope, as well as low PET, whereas eastern populations have adapted to drier and warmer conditions with less precipitation, AET, AI, Alt, and slope, as well as high PET ( Supplementary  Figures 7-15). As a result, these data indicate the potential for intraspecific niche differences in R. pseudodalmatina in geographic (G) space. Therefore, these two regional clades may respond differently to future climate change.
The ecological niche metrics analysis, similar to ENMs results, indicates that WRC and ERC have undergone considerable niche divergence in both climatic and macro-environmental niches (Figures 2-4). As a consequence, the niche overlap observed revealed that the WRC and ERC have little overlap (Schoener's D-metric), and the rejection of the niche equivalence hypothesis confirmed that the two regional clades exist in distinct environmental niche spaces. The results of the similarity test also support these predictions (Figures 2-4). Furthermore, the PCA env results in E-space appear to be credible for the Hyrcanian wood frog, as the two regional clades occupy regions with substantial climate and environmental changes in the Hyrcanian forests (Figures 2-4). These findings indicate that niche divergence may accelerate the evolution of western and eastern R. pseudodalmatina populations within Hyrcanian forests. These results follow Brown and Carnaval (2019) conclusions that niches diverged considerably when niche overlap was low, the equivalent test was significant, and the similarity test was significant or not. However, it should be highlighted that without comprehensive physiological investigation, it may be impossible to determine whether two species' or clades' fundamental niches differ due to divergence, or whether their fundamental habitats are the same but their prospective niches differ due to biotic variables (Brown and Carnaval, 2019). Therefore, further research is needed on these subjects.
The niche overlap results indicated that the WRC and the ERC share a limited hybrid (or sympatric) region. This finding may confirm the results of the previous study (Amiri et al., 2021), which showed that using Circuitscape analysis, the highest connectivity between the two populations is in the central part of the distribution range (i.e., the westernmost margin of the eastern distribution range and the easternmost margin of the western distribution). According to Amiri et al. (2021), this connectivity might be due to a secondary post-glacial connection. However, one possible reason is that the WRC and ERC have similar environmental requirements in areas where admixture between the regional clades is possible (Hu and Jiang, 2012). This approach could be beneficial for identifying and designating sample zones in other species where allopatric lineages meet and exchange genes, especially in those with a deep phylogenic structure (Peñalver-Alcázar et al., 2021). Furthermore, it emphasizes that connecting multiple sources of knowledge (e.g., ecological and historical biogeography, phylogeographic and phylogenetic techniques) aids in improving the species' evolutionary and ecological understanding (Wiens and Donoghue, 2004).

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because this study is modeling and does not require sampling.

AUTHOR CONTRIBUTIONS
SV conceived and designed the study, prepared the data, interpreted the data, made the figures, and wrote, finalized, and revised the manuscript. SV and SR analyzed the data. All authors contributed to the article and approved the submitted version.