We collected occurrence records of the Chrysanthemum zawadskii complex from our field surveys and relevant databases, such as the Chinese Virtual Herbarium and the Global Biodiversity Information Facility ( Figure 2 ). All the occurrences were manually checked to remove incorrect records and records without herbarium images. To conduct ENM analyses, duplicate records per raster cell (2.5 arc-minute) were filtered out with ENMtools v. 1.4.3 (Warren et al., 2010). In addition, the records of four widespread species (i.e., C. zawadskii, C. maximowiczii, C. naktongense, and C. chanetii) were further thinned using the ‘spThin’ v. 0.2.0 R package (Aiello-Lammens et al., 2015) in R v. 4.1.3 (R Core Team, 2022) to keep pairwise distances between the records greater than 10 km for avoiding spatial sampling biases. Overall, we retained 197 points of occurrence that covered the entire distribution range of this complex in China ( Table S1 , Figure 1 ).

We extracted 19 current climatic variables from WorldClim v. 1.4 at 2.5 arc-minute resolution (Hijmans et al., 2005) for ENM analyses and reduced the geographical extent of climatic layers to our study area. To avoid potential multicollinearity effects, we calculated Pearson’s correlation coefficients between the climatic variables with ENMTools. We then selected four relatively independent variables (Pearson’s coefficient < 0.8), including annual mean temperature, temperature annual range, annual precipitation, and precipitation seasonality, for further analyses ( Table 1 ).

To provide potential distributions of species for subsequent niche comparisons, we generated ENMs for each species and lineage of the complex. According to the phylogeny of the C. zawadskii complex (Lu et al., 2022), each of two narrowly distributed species, C. oreastrum and C. hypargyrum, contains two distinct geographical lineages. Therefore, we conducted ENM analyses for both species and intraspecific lineages in this study ( Figure 2 ). We here chose the well-tuned boosted regression trees (BRT) model, which is regarded as one of the top modeling algorithms for ENM analyses (Valavi et al., 2022). We first generated 10,000 pseudo-absences for this complex based on the “random” strategy in our study area with the ‘biomod2’ v. 4.1.2 R package (Thuiller et al., 2022). Second, we tuned BRT models for each species and lineage of the complex. In detail, three hyperparameters of the BRT model (i.e., “interaction.depth”, “shrinkage”, and “n.trees”) were tuned following the guide of Elith et al. (2008) in the ‘SDMtune’ v. 1.1.6 R package (Vignali et al., 2020). We set five alternative values for “interaction.depth” (1, 3, 5, 7, 9) and three values for “shrinkage” (0.001, 0.005, 0.01), and values of “n.trees” were set at a range from 5,000 to 10,000 in increments of 500. We selected the best combination of those hyperparameters according to the area under the receiver operator characteristic curve (AUC) metric. After model tuning, we finally conducted BRT models for each species and lineage to predict their current potential distributions in ‘biomod2’. More specifically, we set other model parameters as default values (“distribution”: bernoulli; “bag.fraction”: 0.5) in ‘biomod2’ and ran the well-tuned BRT model for each species or lineage with ten replicates of cross-validation. Each replicate of models was trained with 80% of the occurrences and tested based on the remaining data. We assessed the performance of each model according to the true skill statistic (TSS) and AUC metrics. Given that TSS values > 0.7 indicate high model performances, we chose the models whose TSS values are greater than 0.7 to generate the ensemble model for each species or lineage to make subsequent projections. The ensemble models of each species or lineage were generated by using the probability-weighted mean method, and then those ensemble models were projected into the current climatic layers to predict the potential geographic distribution of species in this complex.

To compare the niches among different species or lineages of this complex, we quantified their niche overlap and breadth from two different perspectives: in environmental space (E-space) and geographic space (G-space), respectively ( Figure 2 ). Niche comparison in G-space is based on geographic ranges of species which are often inferred by ENMs in many studies, while niche comparison in E-space often uses the environmental information obtained directly from the occurrence locations to quantify the niches of species (Brown and Carnaval, 2019).

For niche comparisons in G-space, we first quantified the niche breadth of each species or lineage using Levins’ B metric, and a higher B value represents a broader ecological niche (Levins, 1968). Then, niche overlap between different species was estimated by Schoener’s D metric (Schoener, 1968) in which high D values indicate that different species may occupy similar niches. It is important to note that the niche overlap and niche breadth in G-space could be overestimated if the areas with very low habitat suitability were included in the calculation (Rödder and Engler, 2011). Thus, we need a threshold to determine the presence or absence of species per raster cell. The “maximizing the sum of sensitivity and specificity” (maxSSS) statistic is regarded as a well-performing threshold selection method for transforming suitability continuous predictions to binary outputs (Liu C. et al., 2013). As a result, we calculated the maxSSS value for each species or lineage with the ‘dismo’ v. 1.3-5 R package (Hijmans et al., 2022), and the habitat suitability value per raster cell below the maxSSS value was set to 0. Finally, we estimated Levins’ B and Schoener’s D based on the binary suitability maps of species and lineages with the ‘ENMtools’ v. 1.0.6 R package (Warren et al., 2021).

We then compared the niches among different species or lineages in E-space. As the Hutchinsonian niche concept describes the niche of a species as an n-dimensional hypervolume, we here used the functional richness (FRic) metric to quantify the niche breadth of species in E-space. FRic is defined as the amount of functional space occupied by the species, and it can be calculated by using probabilistic hypervolumes in trait space (Carmona et al., 2019). Given that some species of the complex, such as C. chanetii, C. naktongense, and C. oreastrum often occur at different elevations on the same mountain, we consider that data from WorldClim may not reflect the niche differences of those species. Therefore, we used ClimateAP v. 2.20 which integrates the local elevation information of species to provide climatic variables with increased spatial accuracy (Wang et al., 2017). We obtained 16 scale-free climatic variables ( Table 1 ) using ClimateAP based on the latitude, longitude, and altitude of species occurrences. We then performed principal component analysis (PCA) for those 16 climatic variables to summarize these variables on two principal component axes with the ‘FactoMineR’ v. 2.6 R package (Lê et al., 2008). We calculated the FRic value for each species and lineage by the ‘TPD’ v. 1.1.0 R package (Carmona et al., 2019) according to the values of PC1 and PC2. To estimate niche overlaps between the species pairs in E-space, we selected the PCA-env approach based on climatic data of presence records and background points (Broennimann et al., 2012). We generated a minimum convex polygon (MCP) for each species or lineage and then defined the geographical backgrounds as MCPs with a one-degree width buffer area. We randomly selected 1,000 points for each species in the background areas and obtained the climatic data of these background points in ClimateAP. We subsequently estimated Schoener’s D and tested the niche similarity for different species pairs in the ‘ecospat’ v. 3.4 R package (Di Cola et al., 2017). In detail, the niche similarity of each species pair was assessed with 1000 randomization tests. Those tests generated a null distribution based on Schoener’s D values between a species and background points of another species in the pair. If the Schoener’s D values of the null distribution are significantly lower than that of the target species pair, these two species may share similar niches than expected by chance.

To reconstruct a dated phylogeny at the species level, we downloaded a set of nuclear loci that included 100 exons with the highest number of parsimony-informative sites (Lu et al., 2022). Because there is no fossil data for Chrysanthemum species, we reconstructed the species tree using the uncorrelated relaxed clock model in *BEAST v. 1.8.4 (Suchard et al., 2018) with the clock rate of 8.37 × 10⁻⁴ – 2.61 × 10⁻³ substitutions/site/million years. This clock rate is based on the standard rate of nucleotide substitution for the Asteraceae (Kay et al., 2006) and three years of the generation time for species of this complex. We selected substitution models based on the Bayesian Information Criterion with ModelFinder (Kalyaanamoorthy et al., 2017), and the optimal model was TN + F + G4. We used the Yule process as the species tree prior and selected the Piecewise linear and constant root as the population size model. We ran two independent analyses with 1 × 10⁷ Markov chain Monte Carlo generations based on a random starting tree, and the trees were sampled every 1,000 generations. The results of these two runs were merged using LogCombiner v. 1.8.4 (Suchard et al., 2018), and then the convergence of these results was evaluated with Tracer v. 1.7 (Rambaut et al., 2018). Finally, a maximum clade credibility tree was summarized by TreeAnnotator v. 1.8.4 (Suchard et al., 2018).

To reconstruct the ancestral niche of the complex, we first tested the phylogenetic signal for the 16 climatic variables obtained from ClimateAP based on the dated species tree in the ‘phylosignal’ v. 1.3 R package (Keck et al., 2016). Three metrics, including Moran’s I, Blomberg’s K, and Pagel’s λ, were calculated, and the significance of these metrics was assessed by 1,000 randomization tests. Moran’s I = 0 indicates that a trait evolved under Brownian motion (BM); I < 0 indicates no phylogenetic signal, and I > 0 suggests that the related species have more similar traits. Similarly, Bloomberg’s K of 1 indicates that trait distribution follows BM; K < 1 indicates no phylogenetic signal, and K > 1 implies that the traits of related species are more similar than those predicted under BM. Pagel’s λ usually ranges between 0 and 1; λ= 1 also indicates that the trait evolved according to the phylogeny, and λ = 0 indicates no phylogenetic signal. However, λ can > 1 when a trait of a related species is more similar than expected under BM (Münkemüller et al., 2012). We then tested the optimal model of niche evolution for those 16 climatic variables using the ‘geiger’ v. 2.0.10 R package (Harmon et al., 2008). Four models including BM, Ornstein–Uhlenbeck (OU), Early Burst (EB), and White Noise (WN) were selected, and we used the Akaike’s Information Corrected Criterion (AICc) value to evaluate model performance. The BM model was the optimal model for all variables ( Table S2 ). Finally, we selected the BM model to estimate the ancestral niches under a maximum likelihood framework by the ‘phytools’ v. 1.2-0 R package (Revell, 2012).

After reconstructing the ancestral niche of species in the complex, we further used phyloclimatic modeling to estimate the ancestral distribution ranges of those species for exploring their historical shifts in ranges ( Figure 2 ). Firstly, we obtained the paleo-climatic data of the same four climatic variables that were used in previous ENM analyses at seven past time periods (i.e., 3.79 Mya, 3.36 Mya, 2.74 Mya, 2.24 Mya, 1.68 Mya, 1.15 Mya, and 0.80 Mya) from the Oscillayers datasets for phyloclimatic modeling (Gamisch, 2019). We focused on those seven time periods because they correspond to the divergence times of species in the complex. All of those variables were at 2.5 arc-minute, and their geographical extents were also reduced to our research area ( Table 1 ). Next, we constructed Bioclim models for each species and lineage of the complex using their current occurrences and climate data in the ‘machuruku’ v. 1.8.3 R package (Guillory and Brown, 2021). The niches of species were described by three parameters including the mean, standard deviation, and skewness of climatic variables, and those parameters were used as inputs for subsequent analyses. We then estimated ancestral climatic parameters of species based on the dated phylogeny and those climatic parameters. Finally, we converted those ancestral climatic parameters into ancestral niche models of species with the ‘machuruku’. All those niche models were projected into corresponding paleo-climate layers of the above seven time periods to reconstruct the ancestral geographic ranges of species in the complex, and those ancestral ranges were converted to binary suitability maps with the ‘machuruku’.

The AUC and TSS values of all the ENMs were greater than 0.90 ( Table S3 ), indicating that all the models performed well. The ENM results for most of the species and lineages in the current period were consistent with their actual distribution in China, except for C. hypargyrum Qinling lineage and C. oreastrum Taihang-Lyuliang lineage ( Figures 1 , S1 ). In addition, the maxSSS values of those species varied from 0.305 (C. zawadskii) to 0.742 (C. mongolicum). Given that BRT models were not effective for those two lineages, we did not estimate their niche breadth in G-space. The Levins’ B values in G-space for other taxa of this complex showed that narrowly distributed species possessed lower Levins’ B values than widespread species ( Table S3 ). The estimations of niche breadth for all species and lineages in E-space also obtained similar results. Species with wider ranges of distribution, such as C. zawadskii, C. naktongense, and C. chanetii, had higher FRic values, and the FRic values of some narrowly distributed species or lineages were even less than one-tenth of the former ( Figure 3 ).

Analyses of niche overlap in the E-space suggested that the niche differentiation greatly varied among the taxa of this complex in which the Schoener’s D values ranged from 0 to 0.462 ( Table 2 ). The niche overlaps between the three widespread species (C. zawadskii, C. naktongense, and C. chanetii) were relatively high (D = 0.279 to 0.462). Moreover, the niche similarity tests showed that C. chanetii and C. zawadskii possessed significantly more similar environmental conditions than expected by chance (p < 0.05, Table S4 ). Schoener’s D values between intraspecific lineages of two narrowly distributed species, i.e., C. hypargyrum and C. oreastrum, were relatively low (D = 0 to 0.089). Meanwhile, the results of the niche similarity tests suggested that Schoener’s D values between intraspecific lineages of those two species were not significantly higher than that of the null distribution ( Table S4 ). C. maximowiczii, a species that occurs in grasslands, showed relatively low niche overlap with other species. Similarly, another endemic species, C. mongolicum, also had little or no overlap with other species. Interestingly, the niche overlap in G-space showed a similar pattern to that in the E-space ( Table S5 ).

The mean divergence times of the main nodes ( Figures 1 , S2 ) estimated under a relaxed clock model were similar to those estimated by Lu et al. (2022). In particular, the C. hypargyrum Qinling lineage diverged first about 4.16 million years ago (Mya), and other lineages or species have then emerged since the mid-Pliocene (ca. 3.78 Mya, Figure S2 ). Based on Moran’s I, Blomberg’s K, and Pagel’s λ, we detected statistically significant phylogenetic signals for three of the 16 climatic variables, including the annual heat: moisture index, precipitation as snow, and Hargreaves climatic moisture deficit ( Table 3 ). We identified the BM model as the evolution model that best fits these climatic variables according to the AICc statistics ( Table S2 ). The results of ancestral state estimation under the BM model suggested that there is niche divergence between the four early-diverged alpine lineages and the core C. zawadskii complex ( Figure 4 ). In comparison to the alpine lineages, the species of core C. zawadskii complex have evolved toward climatic conditions with less precipitation and higher temperature ( Figures 4 , S3 ). The kernel density plots in E-space also suggested that there was differentiation in the climatic conditions occupied by different species ( Figure S4 ).

Our results of phyloclimatic modeling showed that this complex originated in the Qinling mountains and its surrounding areas during the early Pliocene ( Figures 5 , 6 ). These results also showed that this complex diverged with its expansion to northern China along mountain ranges ( Figures 6 , S5 ). Additionally, the geographical ranges of C. hypargyrum and C. oreastrum showed little change and seem to have contracted over evolutionary time ( Figure 5 ). For example, the suitable habitats of these alpine lineages decreased from 2.24 Mya to the present, except for the C. oreastrum Taihang-Lyuliang lineage. In contrast, there were several range expansions for the core C. zawadskii complex during its diversification ( Figure 6 ). From 2.74 to 2.24 Mya, we detected a large shift in the geographical range of the core C. zawadskii complex. Another two range expansions occurred at 1.68 Mya and 1.15 Mya for the core C. zawadskii complex. We then checked the response curves and statistical values of four climatic variables for the core C. zawadskii complex at four different past periods ( Figure 6 , Table 4 ). The 95% confidence interval of those climatic variables indicated that the niche breadth of the core C. zawadskii complex did not increase obviously from 2.74 to 2.24 Mya ( Figure 6 , Table 4 ). Notably, the temperature-related climatic variables had high confidence interval overlap, and their mean values showed little difference between the above two past periods. Conversely, the mean values of the precipitation variables changed more clearly. In addition, as the core C. zawadskii complex diversified, the mean values and 95% confidence interval of those climatic variables greatly varied among species of the core complex.

According to our findings, the territories of Qinling Mountains showed favorable environmental conditions to host the ancestor of C. zawadskii complex in the early Pliocene, making this region eligible as probable land where C. hypargyrum, C. oreastrum, and the core C. zawadskii complex might have evolved. The formation of four alpine geographic lineages could result from mountain uplifts, such as those that took place 4.1–3.4 Mya for the Minshan Mountains and 5.3–3.6 Mya and 2.6 Mya for the Taihang Mountains, and the resultant geographical barriers between alpine populations since the early Pliocene (Zhang et al., 2005; Gong, 2010; Shi et al., 2020). In addition, species of the complex have a relatively weak ability to disperse owing to their achenes that lack pappi or wings (Shi et al., 2011). Both geographical barriers and limited dispersal capacity restricted the gene flow between alpine populations or lineages, and then local selection and/or drift facilitated the formation of four alpine lineages.

Notably, although these alpine lineages have diverged, they still prefer habitats with relatively humid climate conditions and low mean warmest month temperature ( Figures 4 , S3 ), which suggests their phylogenetic niche conservatism in some niche axes. Unlike alpine lineages, species of the core C. zawadskii complex can tolerate hotter and drier environmental conditions. The ancestral niche estimations indicate that the variables related to precipitation exhibit more pronounced divergence within the complex than those of the temperature factors ( Figure 4 , Table 3 ). In fact, the relative importance of temperature and precipitation in driving the distribution and diversification of species is not identical among different plant taxa (Moles et al., 2014; Häkkinen et al., 2022). A recent study indicated that many plants tend to adapt to novel precipitation regimes, whether wetter or drier, rather than adapting to hotter or colder climates, suggesting the importance of precipitation in the niche evolution of species (Häkkinen et al., 2022). In this study, we can deduce that niche differentiation within the complex might be primarily driven by the moisture gradient. According to our results, the common progenitor of the core C. zawadskii complex split from the alpine lineages at approximately 2.73 Mya, and then its distribution range extended to the Chinese Loess Plateau (CLP) and surrounding areas ( Figure 6 ). Previous studies showed that there was a major change in the precipitation gradient on the CLP during the Pliocene-Quaternary boundary, and the climate became dry and cold (Wu et al., 2007; Nie et al., 2013; Peng et al., 2018). The progenitors of the core C. zawadskii complex could have quickly adapted to such a change and then expanded into drier habitats. Additionally, analyses of paleosol pollen assemblages from the CLP suggested that open habitats increased during this period (Hui et al., 2021). Species of this complex generally occur in open habitats, such as grassland, slope, and forest edges (Shi et al., 2011), and such changes could provide more ecological opportunities for range expansion and diversification of the core C. zawadskii complex.

Comparative analyses of niches showed that there exists a distinct niche divergence among species of the complex ( Table 2 , Figure 3 ). Interspecific variations in the phenotypic traits, habitats, and distribution ranges of the C. zawadskii complex implied that niche evolution occurred (Lu et al., 2022). Multiple niche analyses in this study supported that shifts in niche breadth occurred in the complex. For decades, there was a consensus that specialization is an evolutionary dead end (Moran, 1988; Schluter, 2000; Vamosi et al., 2014). Nonetheless, a growing number of studies have indicated that such a “consensus” might not be general, and that it is possible to transition from specialists to generalists (Armbruster and Baldwin, 1998; Dennis et al., 2011; Donoghue and Edwards, 2014; Vamosi et al., 2014; Day et al., 2016; Sexton et al., 2017; Rombaut et al., 2022). Within this complex, early divergent alpine lineages possess narrow niches, and they are more likely to be specialists than generalists. Another small-range species, C. mongolicum, also has a narrow niche, and it is sympatric with its sister species, C. chanetii with a large range (Lu et al., 2022). Taking together the findings described above and the climatic tolerances of ancestral nodes under the phylogenetic framework, we found that specialists can be converted to generalists and vice versa in this complex. Notably, niche comparison in G-space could be affected by the results of ENMs, so such comparisons may not be applicable to some rare species. By comparison, niche comparisons in E-space based on the environmental information of occurrence locations appear to be more stable and reliable.

Climate variability is often regarded as a critical factor in determining the niche breadth of a species (Quintero and Wiens, 2013). Many studies examined the relationship among species range, niche breadth, and habitat climate variability, and they found that species that originated from regions with variable climates could evolve broader niche breadths than their allied species in stable habitats (Stevens, 1989; Addo-Bediako et al., 2000; Quintero and Wiens, 2013; Pintor et al., 2015; Li et al., 2016; Mumladze et al., 2017). Reconstruction of ancestral niche indicated that species of the core C. zawadskii complex have evolved broader niche breadths since 2.24 Mya ( Figure 6 ). During the early Pleistocene (2.2–1.7 Mya), the East Asian summer monsoon (EASM) gradually intensified in the CLP (Hui et al., 2021), and the dramatic seasonal climatic fluctuations in this region enhanced the niche breadth of the core C. zawadskii complex. Nonetheless, the alpine lineages did not evolve increased niches or undergo range expansion ( Figure 5 ). We hypothesized that alpine lineages could survive in situ mountain refugia with stable climates and adapt to alpine habitats rather than expanding to lower altitudes because of the geographical barriers and their limited ability to disperse (Bai et al., 2010; Zhao et al., 2013; Guo et al., 2014).

Speciation is usually a lengthy process with niche evolution and range changes (Stankowski and Ravinet, 2021), and the current pattern of distribution and niche of species represents the consequences of this process. When we integrate snapshots of niches and ranges of species across different historical periods, the complete picture of an eco-evolutionary process that shapes species distribution can be reconstructed. The progenitor of this complex probably adapted to relatively cold climates and diverged from the other Chrysanthemum lineages with the uplift of Qinba mountains since the late Miocene (Liu J. et al., 2013; Meng, 2017; Shi et al., 2020; Lu et al., 2022). Subsequently, the complex spread west- and northwards along the mountain ranges to northern China, and then rode the rising mountains and occupied alpine habitats. Some lineages (C. hypargyrum and C. oreastrum) adapted to survive in situ, and others expanded into drier climates in northern China since the early Pleistocene. Changes in vegetation and the EASM evolution in the CLP during the early Pleistocene enhanced broader niches and larger range sizes of the core C. zawadskii complex.

The spatiotemporal dynamics of the niches and ranges of the C. zawadskii complex indicate that the coupled changes between ranges and niches drive the range expansion of the complex, which is more consistent with the range shift-niche breadth hypothesis (Lancaster, 2022). Climate fluctuation and the decrease of interspecific competition since the Pleistocene provided an opportunity for this herbaceous species complex to expand its range (Wu et al., 2007; Nie et al., 2013; Peng et al., 2018; Hui et al., 2021), and the niche breadth increased in pace with the expansion of its range into areas that have a highly variable climate, which further promoted the range expansion. Notably, even though the species are closely related, the patterns of range shifts are diverse across taxa of the complex, which suggests that historical biogeographic processes also affect the distribution of species.

In sum, by integrating the findings from phyloclimatic modeling, niche comparisons, and biogeographic analyses, we can better understand the processes and mechanisms behind the changes in geographic range and niche evolution. Nonetheless, the reconstruction of ancestral distribution ranges by phyloclimatic modeling is still in the early developmental phase, and some factors, such as the reliability of paleoclimatic data and robustness of phylogeny, could dramatically influence the results (Gamisch, 2019; Rivera et al., 2020; Guillory and Brown, 2021). Although we weighed and chose the appropriate climatic variables for the complex in the current period, we could not say that the paleoclimatic data of the same variables are the best fit for estimating the ancestral distribution of the target species. In addition, we here focused on the niche of species on a larger scale, but microhabitats could also influence the niche differentiation of species and their distribution (Lawson et al., 2014). Therefore, we need to consider more factors that affect niche evolution and range dynamics of species and incorporate more data in our future works.