Climatic‐niche evolution with key morphological innovations across clades within Scutiger boulengeri (Anura: Megophryidae)

Abstract The studies of climatic‐niche shifts over evolutionary time accompanied by key morphological innovations have attracted the interest of many researchers recently. We applied ecological niche models (ENMs), ordination method (environment principal component analyses; PCA‐env), combined phylogenetic comparative methods (PCMs), and phylogenetic generalized least squares (PGLS) regression methods to analyze the realized niche dynamics and correspondingly key morphological innovations across clades within Scutiger boulengeri throughout their distributions in Qinghai–Tibet Plateau (QTP) margins of China. Our results show there are six clades in S. boulengeri and obvious niche divergences caused by niche expansion in three clades. Moreover, in our system, niche expansion is more popular than niche unfilling into novel environmental conditions. Annual mean temperature, annual precipitation, and precipitation of driest month may contribute to such a shift. In addition, we identified several key climatic factors and morphological traits that tend to be associated with niche expansion in S. boulengeri clades correspondingly. We found phenotypic plasticity [i.e., length of lower arm and hand (LAHL), hind‐limb length (HLL), and foot length (FL)] and evolutionary changes [i.e., snout–vent length (SVL)] may together contribute to niche expansion toward adapting novel niche, which provides us a potential pattern of how a colonizing toad might seed a novel habitat to begin the process of speciation and finally adaptive radiation. For these reasons, persistent phylogeographic divisions and accompanying divergences in niche occupancy and morphological adaption suggest that for future studies, distinct genetic structure and morphological changes corresponding to each genetic clade should be included in modeling niche evolution dynamics, but not just constructed at the species level.


| INTRODUC TI ON
A major goal of ecology is in the inspection of niche evolution dynamics to explain rapid lineage diversification and mechanisms of morphological evolution across clades, especially in complex mountainous regions. As reviewed by Smith et al. (2018), most niche models have been constructed at the species level to model a species responding to the environment as a single undifferentiated entity. For example, several studies have found evidence for climatic-niche conservatism among species (Crisp et al., 2009, Kozak & Wiens, 2006, Kozak and Wiens, 2010, Liu et al., 2020a, while others have shown evidence for niche divergence (Evans et al., 2009;Graham et al., 2004;Hu et al., 2015;Knouft et al., 2006). Either way, these practices ignore whether occurrence data represent a single evolutionary entity or a collection of evolutionary lineages that can vary in age, evolutionary independence, and genetic distinctiveness (Pearman et al., 2010). Moreover, the potential effects and implicit meanings of intraspecific niche evolution dynamics across clades within species level are seldom known (Tingley et al., 2016).
In fact, niches of species or clades clearly do evolve, and niche shifts in range limits as a result of such evolution (Peterson & Holt, 2003). Both ecological (available empty niches) and evolutionary changes (genetic drift or through selection) can potentially allow a species or clade to shift into a novel niche, and an observed shift can equally result from a change of the realized niche and the fundamental niche (Broennimann et al., 2007). Realized niche shifts between native and non-native populations can be accurately evaluated by niche expansion (i.e., species colonizing novel environmental conditions in their non-native range Petitpierre et al., 2012;Tingley et al., 2016) and niche unfilling (i.e., partial filling of the native niche in the native range, Petitpierre et al., 2012). The fundamental niche depicts the ecophysiological requirements of species, which can be viewed as the envelope of environmental (abiotic) conditions allowing populations to sustain themselves in an n-dimensional environmental space (Guisan et al., 2014;Soberón, 2007). In addition, over some temporal and spatial scales, intraspecific niche evolution and ecological innovation have taken place, such as in Mexican birds (Peterson & Holt, 2003). A growing number of cases indicate the evolutionary shifts occurred in range limits with rapidly changing environments (Davis & Shaw, 2001;Evans et al., 2009;Peterson & Holt, 2003;Thomas et al., 2001). Moreover, researchers have documented morphological evolution is strongly influenced by ecological niche shifts in passerine birds (Alström et al., 2015), chestnut-capped brushfinches (Moreno-Contreras et al., 2020), Eurasian perches (such as Perca fluviatilis, Bartels et al., 2012), and bivalved scallops (Sherratt et al., 2017).
Recently, combining molecular information and niche evolution models to analyze niche shifts has provided new insights into the roles of abiotic climate and geographical conditions in shaping range limits. There are three main reasons to incorporate evolutionary processes into niche modeling to assess niche and morphological evolution dynamics across clades within Scutiger boulengeri under climate change. First, the existence of cryptic species and frequent local adaptation suggest that cryptic niche architecture exists within the species-level taxa that are the focus of studies of clade evolution process, ecological niche, and biotic responses to climate change (Pearman et al., 2010). Second, persistent morphological divergence is caused by genetic drift or through selection under local adaptation to environmental heterogeneity. Third, niche models, pooled from the entire range of the species, assume that species respond to the environment as an undifferentiated entity along their entire distribution but underestimating differences between distinct niches caused by range limits and local adaption (Peterson et al., 2011). In fact, spatial heterogeneity in environments coupled with reduced gene flow (resulted from intraspecific competition or dispersal limitations) can encourage local adaptation, leading to divergence in niches among closely related lineages (Smith et al., 2018). However, numerous researchers constructed niche models just by pooling but ignoring intraspecific lineages for widely distributed species owning phylogeographical structures may lose sight of considerable variation in morphological, physiological, and life-history traits under niche evolution dynamics across clades within species (Barria et al., 2020).
The Qinghai-Tibetan Plateau (QTP)-the largest continental highland on Earth-is a major barrier to airflow in the atmosphere, which triggers the onset of the Indian summer monsoon (Molnar et al., 1993). Tibet continuously grew northward over millions of years in response to the thickening of Earth's crust associated with the collision of the Indian and Asian continental plates (Harrison et al., 1992), which is a long-standing topographic feature that arose from the collision between India and Asia (Rowley & Currie, 2006).
The orogeny of high mountain ranges separating deep valleys might have created geographical barriers reducing gene flow between isolated populations and promoted allopatric divergence (Favre et al., 2015). Meanwhile, novel environmental spaces released from biotic and abiotic constraints (Callaway & Maron, 2006;Hierro et al., 2005) would have provided key opportunities for occupation of novel niches especially in the early stages of clade divergence.
Scutiger boulengeri, an endemic Tibetan toad occurring in mountain streams from the South-Tibetan (Hofmann et al., 2017), has a wide range of distributions along the eastern and southern slopes of the QTP at elevations between 2,400 and 5,270 meters above sea level Subba et al., 2015). Several geographically structured haplotypes have also been identified using mitochondrial DNA and grouped into 3 major clades due to incomplete sampling (Li et al., 2009). Intraspecific clade diversity of S. boulengeri implies each clade resulted from unique patterns of limited migration, isolation, and local adaptation (Potter et al., 2013) across drainages along the margins of QTP. Such a promising case presents us an attractive system to study the link between intraspecific niche evolution and phenotypic evolution dynamics. The capacity for rapid phenotypic evolution may directly facilitate species diversification by increasing the ability of a radiating clade to exploit ecological opportunities (Parent & Crespi, 2009). Moreover, it also provides an ideal model to study phenotypic plasticity, which as the main means to cope with changing ambient conditions on a shorter time-scale in ectotherms (Angilletta, 2009;Seebacher & Franklin, 2011).
Given unique intraspecific lineages, the potential for distinctive climate niches and the need to model niche evolution dynamics of S. boulengeri at the intraspecific clade level have been recognized.
Herein, based on a set of climatic and morphological data, we applied multiple robust models and methods for incorporating evolutionary processes into niche modeling to assess niche and morphological dynamics across clades within S. boulengeri. Specifically, we focus to address four key issues: (a) Is there niche divergence caused by niche shifts across clades? (b) Is such a divergence caused by niche unfilling or niche expansion? (c) Which climate variables contribute most to such niche evolution dynamics? (d) Is there related trait evolution accompanied by a shifted niche when controlling for phylogenetic relatedness? We hypothesize that genetically isolated S. boulengeri clades would exhibit clearly segregated niche patterns and corresponding morphological variations in this system.

| Phylogenetic analysis
Based on previous studies (Hofmann et al., 2017;Li et al., 2009) and our own field works in recent years, we compiled 2 published mtDNA cytochrome b (cytb) genomes (GenBank IDs FJ463132 and EU180928) plus 5 newly obtained cytb genomes (GenBank IDs MW600725-MW600729). To construct a phylogeny for S. boulengeri clades, we used MEGA-X (Kumar et al., 2018) to align six selected mtDNA cytb genomes with one genome from outgroup Oreolalax omeimontis. Phylogenetic trees were constructed separately by using maximum likelihood (ML) and Bayesian inference (BI) methods, both of which were implemented in PhyloSuite v1.2.1 (Zhang et al., 2020).
The best-fit BIC substitution model (TPM2 + F + G4) was selected in ModelFinder (Kalyaanamoorthy et al., 2017). Divergence time for the reconstructed trees was estimated with the RelTime ML method using MEGA-X (Kumar et al., 2018). Based on previous researches, we choose the most recent common ancestor (MRCA) of Scutiger and Oreolalax (53 Ma) as the calibration point (Hofmann et al., 2017).

| Occurrence and environmental data
We obtained occurrence records for Scutiger boulengeri from our own field works and the published literatures (Hofmann et al., 2017;Li et al., 2009). Localities cover Himalayas, QTP, Hengduan Mountains, Min Mountains, and adjacent mountains ( Figure 1).
We compiled 19 BIOCLIM variables and elevation for each period from the WorldClim database with a resolution of 30 s (~1 km) for each environment layer (Fick & Hijmans, 2017). We included the dissimilarity of the enhanced vegetation index variable, drawn from the Global Habitat Heterogeneity project (Tuanmu & Jetz, 2015).
Because strong colinearity between environmental variables could inflate model accuracy Veloz, 2009), we examined pairwise correlations among the 21 variables within each clade distribution. We reduced autocorrelation of input environmental data by removing highly correlated variables with the threshold of Pearson's correlation tests |r| > .8 (Dormann et al., 2013). Finally, we chose 8 variables with lower correlation for subsequent analyses.

| Model evaluation of ecological niche models (ENMs)
We used the ENMeval R package to facilitate increased strictness in the development of Maxent models (Muscarella et al., 2014 Akaike information criterion corrected (AICc) for small sample sizes reflects both model goodness-of-fit and complexity. The model with the lowest △AICc value (i.e., △AICc = AICc − AIC min = 0) is considered as the best model out of the current suite of models. We explored models with regularization multiplier (RM) values ranging from 0.5 to 4.0 (increments of 0.5) and with six different feature classes (FCs) combinations as suggested by Muscarella et al. (2014).

| Suitable habitat prediction
We used a maximum entropy modeling algorithm implemented in the program Maxent v.3.4.1 (Phillips et al., 2006;Phillips & Dudík, 2008) to predict suitable habitat. Maxent uses environmental variables from localities at which a species has been documented previously to build a predictive model of where else the clades may occur due to the presence of similar environmental conditions (Elith et al., 2011).
To assess model performance, we calculated the average value of the area under the receiver operating characteristic curve (AUC) for training and testing datasets (Swets, 1988), AUC takes on values ranging from 0.5 (no better discrimination than random) to 1 (perfect discrimination).

| Niche overlap and null hypothesis test
One key assumption for applying ENMs is that species' niche changes very slowly across space and time (Warren et al., 2008). These tests are based on two similarity metrics (Warren's I and Schoener's D); we calculated these metrics in ENMTools v1.4.3, using 100 replicates to generate a pseudoreplicated null distribution (Warren et al., 2010).
The null hypothesis of niche equivalency is rejected when empirical values are significantly less than the critical values for both the niche equivalency and similarity tests (Warren et al., 2010).
Climate niche overlap in E-space between lineages of S. boulengeri was estimated using the PCA-env approach proposed by Broennimann et al. (2012). An unbiased estimate of the Schoener's D metric can be calculated for our data and is ensured to be independent of the resolution of the grid; statistical confidence in niche overlaps was then tested through a bidirections niche similarity test .

| Niche expansion or unfilling
Schoener's D on species occupancy disentangles climate availability and the extent of niche divergence of clade pairs, but does not take into account the difference between partial filling and expansion . However, expansions measured can characterize true niche shifts, when native and non-native ranges overlapped in climatic space, Following Petitpierre et al. (2012), three categories were considered in: (a) stable environments where species occurs in both ranges, (b) unfilled environments where species occur only in the native range, and (c) expansion environments where the species occur only in the non-native range.

| Niche evolution
We used phytools R package (Revell, 2013) to visualize niche evolution throughout the phylogeny. We calculated the K value for a given trait and phylogeny; phytools package provides a randomization test to assess the significance of the observed K value (Revell, 2020).
Calculations were conducted using the geiger package (Harmon et al., 2007), and the best-fitting model was chosen using the Akaike information criterion corrected (AICc) for small sample sizes and Akaike weights (ω) (Wagenmakers & Farrell., 2004).  Fei et al. (2005). Scutiger boulengeri is characterized by one or two pairs of keratinized spine patches on the chests of males, a reduced columella, and the absence of a tympanum . In data analysis, we removed females for their insufficient quantity. Prior to all statistical analyses, the variables were log-transformed to better meet the requirements of normality and homogeneity (Rabosky & Adams, 2012).

| Phylogenetic comparative methods (PCMs) and trait correlative analyses
Biologists have long recognized that closely related species are generally more similar to one another than they are to more distantly related, which is often termed phylogenetic conservatism (Martins & Hansen, 1997). Phylogenetic signals can be considered as the degree to which similarity in trait values between species can be predicted upon their relatedness (Harvey & Rambaut, 2000). If there are phylogenetic signals in the data, then PCMs are necessary for robust statistical analyses of trait correlations. To address the relationships between morphology and climate variables caused by niche evolution, we used phylogenetic generalized least squares (PGLS, Grafen, 1989;Martins & Hansen, 1997). The workflow, incorporating all the methods and processes for modeling climatic-niche evolution dynamics and key morphological changes along phylogenetic clades in S. boulengeri, is presented in the Supplementary material Appendix Figure A1.

| Occurrence and environmental data
Our final dataset comprises 96 georeferenced occurrence records:  (Table 1). PC1 is mainly represented by the increased temperature and reduced elevation, while PC2 explains increased precipitation, PC3 explains mean diurnal range and precipitation seasonality, PC4 and PC5 both represent Min temperature of coldest month and temperature annual range. The top three principal components are presented in Figure 3a; for more details, please see Table 1.

| Prediction of the Maxent distribution
We chose the best model and found the corresponding RM and FC parameters for Maxent model to facilitate increased rigor in the development of Maxent models ( Table 2). The distributions based on Maxent across S. boulengeri clades are characterized by high AUC statistics, indicating that these ENMs successfully discriminate real occurrences from background locations. The jackknife tests on variable importance for S. boulengeri clades reveal that precipitation of driest month in E. A, E. B, and W. a clades, while annual precipitation in E. C and E. D and annual mean temperature in W. b produce the greatest decrease in gain when excluded from the model, suggesting these climate factors limit distributions of clades correspondingly, which are likely the most important reasons for the next step niche divergence.

| Hypothesis tests based on ENMs and PCAenv approaches
ENM-based niche equivalency tests reveal that four paired compari- Only four paired comparisons show niche overlap categorized as a moderate overlap (0.307-0.527). On the other hand, the background similarity tests indicate generally little niche similarity among the pairwise comparisons in six clades (Table 3 and Supplementary material Appendix Figure A5). The ordination null tests of niche similarity show that niches are less similar than random expectations in 30% (five paired comparison cases) in reciprocal directions, while 26% (four paired comparisons) are significantly not to reject the null hypothesis of niche conservatism in bidirections.

| Niche expansion and unfilling dynamics
As predicted, examining patterns of niche expansion and niche unfilling demonstrate a gradient of realized niche change across clades in their shifted ranges (Supplementary material Appendix Figure A6).
There is considerable evidence of expansion when comparing the realized niche of clade to its native niche, such as E. A clade shows 99.7% niche expansion (vs. native clade E. D) (Figure 3h). From our results (Table 3), there is nearly 50% (3/6; clade E. A, E. B and E. C) of the clades' non-native niche exists in climates that are less occupied in its native range (i.e., niche expansion). While 33% (2/6; clade E. D and W. a) of the clades' native niches remain stable, only one clade W. b can be viewed as typical niche unfilling (Figure 3i).

| Niche evolution
The results of phylogenetic signal tests based on Blomberg's K show values in PC1 and PC2 are less than 1, best-fitted to WN model, which suggests that PC1 and PC2 dimensional climate changes are independent of phylogenetic relationships (data with no covariance structure among clades). But the values of more than 1 in PC3, PC4, and PC5 indicate PC3, PC4, and PC5 dimensional changes have closely related phylogenetic relationships (Table 4). The comparisons of model fit based on the AICc values and AICc weights (ω) indicate that the WN model is preferred in PC1 and PC2, where climatic components change in PC1 and PC2 regardless of shared ancestry between clades. Although niche conservatism is maintained in some clades for PC1 and PC2 (Figure 2c and Table 4), our tests of PC3, PC4, and PC5 fitted to a BM model, suggesting that following a divergence event along these climate PCs, clades branches may be subject to these variable environmental conditions, such as promoting the evolution of different climatic tolerances, which may be accumulated independently from ancestral ones. relatedness. Following a divergence event along phylogeny, that is less labile than expected under a BM model of evolution, clades branches may be subject to phylogeny, such as promoting the evolution of different morphological traits under climatic tolerances, which may be accumulated independently from ancestral ones, while the others fitted a WN model, which seems to be unrelated with shared ancestry but preference of an adaptively phenotypic plasticity (Table 5), which may be closely related to habitat types adapted by distinct genetic clades.    (p > .05). For other best-fitted models, please see Table 6.

| D ISCUSS I ON
Combining molecular information and niche evolution models, our results show there are six clades contained in S. boulengeri by mtDNA genetic marker. Geographic structure identified from mtDNA suggests some clades resulted from unique patterns of migration, isolation, or local adaptation (Potter et al., 2013). There are three clades departed from their native niche and show niche divergence, and all the three clades shift their niches due to niche expansion but not by niche unfilling. Some climate variables may contribute to such a shift, such as annual mean temperature and annual precipitation and precipitation of driest month according to jackknife test of variable importance. We found evolutionary changes (i.e., SVL) and phenotypic plasticity (i.e., LAHL, HLL, and FL) may together contribute to niche expansion toward adapting novel niche.
Our results agree with previous progressive uplifts of Tibet (Mulch & Chamberlain, 2006 increasing niche change. Similarly, Tingley et al. (2014) found that a native-range ENMs under-predicted the extent of the species' Australian invasion. An effective strategy to improve model predictability is to develop species-specific models or models for functional groups (Guisan & Thuiller, 2005).
In addition, we found E. A, E. B, and E. C clades within S. boulengeri have an obvious expansion of climatic niche (Figure 3h and Supplementary material Appendix Figure A6). Intriguingly, they shift beyond the realized niche with the new conditions but still overlap the fundamental niche, which provides positive proof that the niche conservatism hypothesis (Kozak & Wiens, 2006;Wiens, 2004) and niche divergence hypothesis (Evans et al., 2009;Graham et al., 2004) are not contradictory. We found a gradient of realized niche change in the non-native ranges across clades within S. boulengeri: niche stasis in E. D (96.3%) and W. a (96.6%), niche unfilling in W. b (75.2%), and niche expansion (vs. E. D and W. a separately) in E.
A (mean = 99.3%), E. B (mean = 97.1%), and E. C (mean = 90%). Our results seem to be inconsistent with the conclusion of previous studies in Petitpierre et al. (2012), in which realized niche shifts between the native and non-native ranges were largely due to niche unfilling.
Our results are also different from the results of cane toad in Tingley et al. (2014): the shift in the realized niche of the cane toad Rhinella marina was solely due to niche expansion. In our results, niche expansion into novel environments is more popular than niche unfilling, suggesting that our niche divergence due to niche expansion in the shifted range and thus represents true niche changes. Why did E. D and W. a fail to fill its fundamental niche in its native range? One possibility is that the presence of closely related species (S. glandulatus and/or S. mammatus) might have prevented S. boulengeri from colonizing suitable environments south of its present range. Indeed, previous study found there were low rates of interspecific hybridization . We cannot exclude dispersal limitation in the native range as a possible contributing factor, such as Jinsha River and Yalong River (Li et al., 2009), which can also enforce stable parapatric range boundaries. Future studies will be able to test this hypothesis using laboratory or field experiments.
Numerous examples of rapid adaptation in non-native niches suggest that rapid evolution may be common during invasions in species level (Alström et al., 2015;Bartels et al., 2012;Sherratt et al., 2017). The degree to which species adapt to novel environments is important to a range of topics in ecology and evolution (Warren et al., 2008), but is of special concern for the study of intraspecific niche evolution (Tingley et al., 2016). In our study, niche divergence caused by niche expansion indeed accompanied key morphological innovations of preadaption in novel climates versus niche unfilling and stability. Our finding of significant phylogenetic signals in SVL (Table 5) and elevation (Table 4 and Figure 2c) indicates that these acquired data are not random and our results are robust.
Furthermore, our findings of significant phylogenetic signals in these traits are consistent with previous studies (Blomberg et al., 2003;Freckleton et al., 2002;Oufiero et al., 2011). We found that eleva-  of warmest month (AIC = 29.31; p = .037) are significantly negative predictors of SVL under phylogenetic models, which suggest S. boulengeri toads from warmer and more arid environments tend to be larger, which is in concert with true records in our field works.
Several factors may underlie the observed pattern of SVL variations in S. boulengeri clades. One possibility pertains to the expected relationship between fasting endurance and SVL (Mautz, 1982). The second possibility is ecological release in novel shifted areas may allow for larger SVL (Losos & Queiroz, 1997;Yoder et al., 2010).
The third is likely that maintenance of preferred body temperature influences the evolution of SVL (Oufiero et al., 2011 These findings and potential issues provide us important inspiration and guidance for our future research.

| CON CLUS ION
Combining and analyzing distinct genetic clades from different geographic areas with correlative niche models and morphological evolution models, we have shown considerable variations in the degree of realized niche expansion and unfilling across the clades within S. boulengeri toads in the Qinghai-Tibet Plateau region. In the case of S. boulengeri toads, niche divergence occurs accompanied by niche expansion rather than niche unfilling, that is, niche expansion is more prevalent than niche unfilling in E. A, E. B, and E. C clades, while niche unfilling presents just in W. b clade.
Meanwhile, niche divergence caused by niche expansion indeed accompanies key morphological innovations of preadaption in novel climates than niche unfilling and stability. Factors such as enlarged body size and enhanced locomotor performance have been shown to increase expansion success by helping toads to cope with novel conditions.
Recognizing true niche shifts accompanied by key morphological innovations do exist; further assessments should seek to understand molecular mechanisms of key morphological innovations and/or related life strategies that have allowed these particular clades to expand their niches dramatically. It would be particularly interesting to use the same framework to test in the future whether the same patterns are found in other organisms, especially for other widespread species or clades.

ACK N OWLED G EM ENTS
We sincerely appreciate constructive suggestion and corrections from the Editor and two reviewers. We thank Dr. Chunlong Liu from

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

O PE N R E S E A RCH BA D G E S
This article has earned an Open Data Badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at https://doi.org/10.5061/ dryad.2ngf1 vhn6.