Niche models at inter- and intraspecific levels reveal hierarchical niche differentiation in midwife toads

Variation and population structure play key roles in the speciation process, but adaptive intraspecific genetic variation is commonly ignored when forecasting species niches. Amphibians serve as excellent models for testing how climate and local adaptations shape species distributions due to physiological and dispersal constraints and long generational times. In this study, we analysed the climatic factors driving the evolution of the genus Alytes at inter- and intraspecific levels that may limit realized niches. We tested for both differences among the five recognized species and among intraspecific clades for three of the species (Alytes obstetricans, A. cisternasii, and A. dickhilleni). We employed ecological niche models with an ordination approach to perform niche overlap analyses and test hypotheses of niche conservatism or divergence. Our results showed strong differences in the environmental variables affecting species climatic requirements. At the interspecific level, tests of equivalence and similarity revealed that sister species were non-identical in their environmental niches, although they neither were entirely dissimilar. This pattern was also consistent at the intraspecific level, with the exception of A. cisternasii, whose clades appeared to have experienced a lower degree of niche divergence than clades of the other species. In conclusion, our results support that Alytes toads, examined at both the intra- and interspecific levels, tend to occupy similar, if not identical, climatic environments.


Methods
Organism presence records and environmental data. The  The habitats of these five species show a large degree of differentiation. A. obstetricans requires areas with high amounts of precipitation and occupy a wide range of habitats, from mountain ranges to crops 21 . A. muletensis is only known in a few localities on northern Mallorca Island 22 , whereas A. maurus is restricted to a few localities in the Rif and Middle Atlas Mountains of Morocco and typically occupies humid sites in karst and steep areas 23 . Regarding the two Iberian endemics, A. dickhilleni is restricted to the Betic region in southeastern Spain 24 . A. cisternasii occupies the central-south-western sector of the Iberian Peninsula, usually between 0 and 700 m 25 , and, in comparison to the other species of this genus, it is associated with a hotter and drier climate 26 .
We used 676 localities to build ecological niche models. A total of 319 A. dickhilleni presence points were surveyed by the authors within all their Andalusian distribution areas. We identified 170 local population locations of A. cisternasii (including our own collection data and data from Amphibiaweb 27 ), 162 of A. obstetricans (our data, Amphibiaweb, and in addition, revised data from www.obser vatio n.org), 14 of A. maurus (our data, Amphibiaweb and Donaire et al. 23 ), and 11 of A. muletensis (our data and Amphibiaweb). All populations were separated by at least 200 m. Figure 1 shows the distribution of points selected for all five species.
To establish intraspecific comparisons among A. cisternasii, A. dickhilleni and A. obstetricans, we used presence data for the genetic lineages reported for these species. Presence data were assigned to lineages according to Gonçalves 9 (the central-eastern lineage of A. obstetricans was excluded because of the low number of presence points). Some local populations could not be assigned to any of the reported lineages, as these locations are in contact zones our outside the areas surveyed in the cited studies, and therefore, these were removed from the data set. We excluded 42 locations of A. dickhilleni and 6 of A. cisternasii. For the geographically restricted A. muletensis and A. maurus, the low number of records and lack of adequate genetic lineage information precluded intraspecific analysis. Even in a well-designed data survey across all of the range of a species, model outputs are sensitive to sampling bias 28  www.nature.com/scientificreports/ and lineages included in the study (random selection by MaxEnt). In addition, intraspecific figures show the presence points that were selected for each lineage and then used in the intraspecific correlative model. We considered climatic and topographic factors to explain the realized distributions and niches of the five Alytes species. Climatic variables were obtained from WorldClim version 2 database at 30 s 29 , and topographic data were obtained from SRTM (https ://www2.jpl.nasa.gov/srtm/datap relim descr iptio ns.html). The study area was fixed to the distribution limits of the Alytes sp. Genus, and in order to avoid a high correlation and redundancy among the predictors, we performed pairwise Pearson correlations, and for r > 0.6, the variable with lower biological relevance was excluded. The final data set included five climatic variables: mean diurnal range, isothermality, minimum temperature of the coldest month, mean temperature of the driest quarter, and precipitation of the wettest month. niche model analyses. We performed two different modelling approaches: ecological niche models (ENMs) and an ordination technique 30 . First, of the different ENMs types, we selected the widely used machine learning method MaxEnt (version 3.4.1) to build the SDMs 31 . The model was fitted using hinge, product, linear, and quadratic features with a maximum of 10,000 background points, 1,000 replicates, and clamping. Models were fitted by using the area under the ROC curve (AUC and ROC represents "receiver operating characteristic" 31 ). We used the Cloglog output format. Although this measure has been extensively used to fit models, its usefulness has been criticized, especially for presence/background models such as MaxEnt. Thus, in addition to AUC, we present the values of its components, sensitivity and/or specificity 32 . To obtain these two indices, the continuous MaxEnt output was transformed into a categorical variable (predicted presence/absence). For this transformation, we applied the threshold of the minimum training presence (the lowest suitability scores associated with the populations of each lineage/species) given in the MaxEnt output sheet. We calculated the specificity from the confusion matrix 33 . This is a conservative and realistic threshold since it may include even small suitability scores whenever the lineages/species are present 34 . In this case, a calculation of sensitivity was not necessary since the applied threshold was fixed at the maximum. In addition, we have followed the Raes and Ter Steege validation method 57 , including 95% I.C AUC values of null models created with random points of the same size of the presences included in our models. ENM AUC values that are higher than their corresponding 95% CI AUC value of the fitted null model, significantly deviate from what would be expected by chance (p < 0.05).
The ordination technique approach was applied to perform niche overlap analyses, either between pairs of the five species or between pairs of the lineages within A. cisternasii, A. dickhilleni and A. obstetricans. We used the tool Ecospat, which incorporates null hypotheses 35 . For these analyses, we performed the following tests: niche equivalency tests (are niches identical?), similarity tests (are niches more similar than expected by chance?), and niche principal component analysis (PCA). As a measurement of the realized niche overlap, we calculated a Schoener's D index through the niche-PCA. This index ranges from 0 to 1 to reflect no overlap to total overlap, respectively 36 . For both the niche equivalency and similarity tests, we used the argument = "greater" (overlap greater than expected by chance) to test the conservatism hypothesis and the argument = "lower" (overlap lower than expected by chance) to test the divergence hypothesis 35 . We performed 1,000 permutations for each analysis. Additionally, we integrated phylogenies (inter-and intraspecific levels) using the age-range correlation function www.nature.com/scientificreports/ of the Phyloclim package 38 . This function is used to test for phylogenetic signals in patterns of niche overlap. Slopes and intercepts derived from a linear model can be used to characterize speciation mode (allopatric versus sympatric) or niche evolution (conservatism versus flexibility) in the clade 39 .
Regarding the choice of the geographical extent (and as a consequence, environmental background), we used the software QGIS 40 to compile and process environmental data using the extension point sampling tool 41 . Analyses were conducted using R 42 . The study area for the interspecific analysis with MaxEnt and ECOSPAT was adapted to the total distribution of the genus to allow comparisons. At the intraspecific level, we used the same area in each cluster of lineages (the species ranges).

Results
At the interspecific level, the Maxent outputs are shown in Fig. 2. The variables with a relatively higher percentage of contribution in the Maxent models were isothermality (62.5%; A. obstetricans), mean temperature of the driest quarter (56.  Table 1. We found minimal niche overlap with significant p values for the equivalency test in the divergence hypothesis, except for the overlap between A. obstetricans and A. maurus. The magnitude and sign of the variables in the principal component plots of the niche are shown in Fig. 3. Maxent outputs for intraspecific A. cisternasii, A. dickhilleni and A. obstetricans indicate differences in climatic suitability among lineages. The model outputs are shown in Fig. 4, and the ENM AUC, 95% I.C. AUC of  The results of niche overlap (Schoener's D) and similarity and equivalency analyses showed low to medium values depending on the pair of lineages compared and the species (see Table 2 for A. dickhilleni and A. cisternasii and Table 3 for A. obstetricans). The south-western lineage of A. dickhilleni has lower overlap values than those of the other lineages, and in the equivalency test, we obtained significant results for the divergence hypothesis compared with the niches of the other two lineages (equivalency test, divergence hypothesis, Table 2). There were no significant divergences or convergences between the niches of A. cisternasii lineages ( Table 2). For A. obstetricans, we found significant values in the equivalency test between most lineage comparisons, but we also obtained several significant results for the niche conservatism hypothesis ( Table 3)   www.nature.com/scientificreports/

Discussion
Our results showed a clear spatial niche segregation when we examined interspecific niche variation in Alytes toads. On the Iberian Peninsula, the distribution outputs of the MaxEnt models fit neatly in some places, although they exhibited a slight overlap in certain areas that coincides with the present sympatric distributional area between A. cisternasii and A. obstetricans. The pattern derived from niche similarity and equivalence tests revealed that each Alytes species occupies a non-identical environmental niche since no significant p values were found for the hypothesis of complete niche overlap (equivalence test), but instead, highly significant distinctions occurred for the divergence hypothesis (equivalence test) for any pair-species comparisons with the exception of the A. obstetricans versus A. maurus contrast. This suggests an evolutionary scenario where niches are less equivalent (identical) than expected by chance in relation to different non-exclusive processes, including local adaptation. In addition, we did not find significant results in the similarity test (conservatism hypothesis), rejecting the hypothesis of retained similarity. We also did not find significance in the case of the divergence hypothesis for the similarity test, indicating no greater than expected divergence. However, rejecting the retained similarity, determining that the divergence was not greater than expected was already relatively divergent. A. cisternasii has been described as the most phylogenetically distant group within the genus Alytes 17 . This is compatible with an evolutionary scenario where the complex formed by A. obstetricans, A. maurus, A. muletensis, and A. dickhilleni share a more recent natural history and, consequently, they could also share similar environmental/ climatic requirements to a high degree. However, our results did not support this prediction, thus suggesting that the niches of these species evolved in a complex scenario, creating a wide diversity of adaptations. A. maurus and A. obstericans presented the widest climatic range (see Fig. 2). In comparison to the other species (with the exception of A. maurus), A. obstetricans inhabits colder areas with higher precipitation. A. cisternasii occupies warmer areas with relatively high precipitation (at least in the wettest period); A. maurus faces a wider range of temperatures and a high rainfall amount but with the widest variable range. It is important to remark that our models shows two separate areas in the ecological niche models (Ecospat) of A.maurus, this being possibly this an artefact due to the lack of knowledge about the distribution or alternatively the isolation of the two extant populations 56 . In turn, in comparison to A. cisternasii and A. muletensis, A. dickhilleni is present in colder and drier areas. Finally, A. muletensis occupies warm and dry areas, although the distribution of this species was much wider in the past 44 ; in addition, the current distribution may be restricted to highly isolated populations due to non-climatic factors such as intense predation pressure by the introduced water snake Natrix maura 43 . This scenario may bias the output of our correlative model, which relies on distributional and climatic factors. A  Table 3. Schoener's D and p values (similarity and specificity for niche conservatism (C) and niche divergence (D) hypotheses of intraspecific lineages for A. obstetricans. *Significant p values (boldface). www.nature.com/scientificreports/  www.nature.com/scientificreports/ process to determine the robustness of our approach could be to implement mechanistic physiologically based models and to examine the congruence of both approaches 55 . When examining the degree of environmental niche evolution at the intraspecific level, we found contrasting patterns across and within species. The four clades of A. cisternasii tended to exhibit slight niche differentiation differences between the four genetically distinct clades, although correlative models showed different predicted distributions of the four lineages, with the western lineage being the most different, as it is associated with more humid conditions, than those the three other lineages that form a complex with reduced niche differentiation, in agreement with the phylogenetic tree proposed by Gonçalves et al. 13 . However, a contrasting pattern was found in A. dickhilleni, whose southern lineage exhibited significant non-equivalence with respect to the other three lineages, reflecting a more separated evolutionary history for this clade. Interestingly, this southern lineage differentiation in climatic requirements, characterized by a reduced diurnal range and mild winters (temperature of the coldest month) whereas the other three lineages exhibited differentiation in their climatic niches, fits well with the observed pattern of genetic divergence, by which the southern clade formed a sister group to a complex containing the other three lineages 12 . Finally, clades of the widespread Alytes obstetricans exhibited the highest diversity in the pattern of climatic niche evolution with high niche conservatism. The two southernmost lineages between the north-western and south-eastern lineages (in the similarity test) suggest a diminished importance of niche divergence for the intraspecific lineages described for this species. In addition, we did not find support for the existence of phylogenetic signals in the age-range correlation tests (both at interspecific and intraspecific levels). This allowed us to consider the flexibility in the niche evolution hypothesis as opposed to niche conservatism reported for other groups 36. Climatic conditions are important factors influencing both the inter-and intraspecific evolution of Alytes and consequently its ecological niche segregation. The evolutionary history of this genus seems to be the result of a combination of vicariant factors influenced by both landscape and geographic factors 17 . Our results for Alytes reinforce the idea that intraspecific variability can be one of the major drivers of biodiversity 46 . Our results also match the conclusions of Maia-Carvalho et al. 9 about the ongoing processes of differentiation in A. obstetricans but provide a more general, wider view of the generation of diversity. Thus, vicariant and geographic barriers explain the current patterns (inter-and intraspecific) of diversification; environmental and geographical factors can act synergistically to drive differentiation at multiple scales. The intraspecific differences may be explained by the most commonly accepted evolutionary alternatives: (A) niche conservatism that may be the consequences of natural selection 47 . We observed a lower intraspecific niche divergence in A. cisternasii than in A. dickhilleni and, at the same time, a higher phylogenetic influence (p value lower). This may possibly be induced by its specialized thermophilic physiology that may constrain uplift dispersion to mountain ranges 48 (Rodríguez-Rodríguez et al. unpublished data). B) Other selective sources affect relatively lower dispersal ability due to ecomorphological constraints. A. cisternasii is the shortest-limbed species in group 25 , and indirectly, this may be associated with a lower dispersal ability than that of the other Alytes species 49,50 .

Alytes obstetricans Schoener's D Equivalency p value (C/D) Similarity p value(C/D)
Regarding the phylogenetic analysis results, we conclude that no phylogenetic signal was detected at either interspecific or intraspecific levels. This fact is congruent with the conclusions the remaining tests that supported the rejection of retained niche similarity (conservatism), suggesting a nom parallelism of phylogenetic inertia and niche evolution.
Our results support a model of hierarchical niche differentiation in midwife toads. This model helps to understand the evolution of this primitive genus of amphibians, but most importantly, this approach has widespread application in conservation biology. First, it demonstrates the need for a modelling strategy based on targets below the species level. Second, it shows that the identification of a relatively reduced series of bioclimatic variables can enable the identification of the most sensitive taxa or lineages 52 . This emphasizes the importance of evolutionary distinctiveness 53 and the need to connect species range projections with the concept of evolutionarily significant units (ESUs 54 ) to prioritize conservation efforts.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author upon reasonable request.