Species Delimitation of Asteropyrum (Ranunculaceae) Based on Morphological, Molecular, and Ecological Variation

Objectively evaluating different lines of evidence within a formalized framework is the most efficient and theoretically grounded approach for defining robust species hypotheses. Asteropyrum Drumm. et Hutch. is a small genus of perennial herb containing two species, A. cavaleriei and A. peltatum. The distinction of these two species mainly lies in the shape and size of leaf blades. However, these characters have been considered labile and could not differentiate the two species reliably. In this study, we investigated the variation of the leaf blades of 28 populations across the whole range of Asteropyrum using the landmark-based geometric morphometrics (GMM), sought genetic gaps within this genus using DNA barcoding, phylogenetic reconstruction and population genetic methods, and compared the predicted ecological niches of the two species. The results showed that the leaf form (shape and size) was overlapped between the two species; barcode gap was not detected within the genus Asteropyrum; and little ecological and geographical differentiation was found between the two taxa. Two genetic clusters detected by population genetic analysis did not match the two morphospecies. The results suggest that there are no distinct boundaries between the two species of Asteropyrum in terms of morphology, genetics and ecology and this present classification should be abandoned. We anticipate that range-wide population genomic studies would properly delineate the species boundaries and help to understand the evolution and speciation within Asteropyrum.


INTRODUCTION
Delimiting species boundaries correctly is crucial to the discovery of life's diversity because it determines whether or not different individual organisms are members of the same entity (Dayrat, 2005). Traditionally, species are mostly established based on morphology which Cain (1954) called "morphospecies." However, under the unified species concept that defines species as segments of population-level evolutionary lineages (De Queiroz, 2007), distinctive morphology, reciprocal monophyly for haplotypes, reproductive isolation, and divergent ecology can arise at different times and in different orders during the process of speciation (De Queiroz, 1998;Leaché et al., 2009). Thus morphospecies can only be seen as hypotheses that should be tested via different approaches and by using multiple lines of evidence (De Queiroz, 2007). Furthermore, some taxonomists describe species typologically that ignores intraspecific variation and morphological plasticity, this practice could lead to oversplitting by misinterpreting individual variants as members of new specific entities (Dayrat, 2005). The oversplit morphospecies may result in an overestimation of biodiversity and then a waste of conservation resources (Hong, 2016). Therefore, morphospecies with great intraspecific variation and morphological plasticity are particularly in need of reevaluation with respect to their morphology, genetics, and ecology.
Leaf morphology is very important to plant taxonomy and systematics and the variation of leaf form has mostly been studied using morphometrics (Jensen, 2003). Traditional morphometrics involves the application of multivariate statistical analyses to collections of distances, angles, or distance ratios (Adams et al., 2004), seldom considering the shape information (Rohlf and Marcus, 1993). Landmark-based geometric morphometrics (GMM), however, provides a viable alternative for analyzing complex shapes in multivariate space by retaining information about the relative spatial arrangements of the landmarks, allowing for the visualization of shape differences among groups (Zelditch et al., 2004). During the last decade, leaf form variability has fruitfully been investigated using GMM to accurately discriminate species and their hybrids (Viscosi and Cardini, 2011).
DNA sequence-based methods for species delimitation have become popular in recent years. DNA barcoding, an approach initially being developed to identify species based on sequences from a short, standardized DNA region (Hebert et al., 2003;Hebert and Gregory, 2005), are proving to be a robust tool to help unveil biodiversity in a fast, accurate way relying on a distance threshold or "barcoding gap" (e.g., Rossini et al., 2016). However, species delimitation based on a single barcode may be an inaccurate portrait of speciation history (Maddison, 1997). DNA barcoding approaches of species delimitation are often complemented with population genetic analyses, especially those based on multiple nuclear loci (e.g., Lu et al., 2021).
Ecological niche modeling (ENM) utilizes associations between environmental variables and known species' occurrence localities to define abiotic conditions within which populations can be maintained (Guisan and Thuiller, 2005). ENM has already been integrated into a broad variety of research disciplines including species delimitation (Sites and Marshall, 2003;Raxworthy et al., 2007). If a set of populations is geographically separated from closely related species with similar ecological niche by unsuitable habitats, this pattern would support the hypothesis of allopatric speciation with niche conservatism (Peterson et al., 1999). Alternatively, if a set of populations occurs under climatic conditions that are different from closely related species, then gene flow between these populations and other species may also be unlikely or restricted, and these populations may represent a distinct species (Wiens and Graham, 2005).  (Fu and Robinson, 2001). Originally, the two species were described as members of Isopyrum L., however, they have simple leaves, differing obviously from the bi-ternately compound leaves of Isopyrum species. Thus, Asteropyrum was later established by Drummond and Hutchinson (1920) to accommodate the species with simple leaves. The leaf blade of A. cavaleriei is deeply five-lobate, with a width of 4-14 cm, but that of A. peltatum is inconspicuously lobate and 1-3.7 cm wide. In addition, the scape of the former (12-20 cm) is much longer than that of the later (6-10 cm) (thereafter the two species are referred to as "morphospecies"). However, with regard to their gross morphology and karyotype, A. peltatum and A. cavaleriei are very similar ( Figure 1A and Supplementary Figure 1), some populations are intermediate in leaf shape and size (Yuan and Yang, 2006). Accordingly, Yuan and Yang (2006) treated the two morphospecies as subspecies of A. peltatum. Although Yuan and Yang (2006) recognized that the morphological differentiation between A. peltatum and A. cavaleriei is obscure, they had not conducted statistical analyses on the leaf morphology and had not evaluated their genetic and ecological differentiation. In our pilot investigations, we even found that the two morphospecies occur in the same populations occasionally. These field observations cast new doubt on the taxonomic status of the two morphospecies or even morpho-subspecies. Therefore, the boundary between the two Asteropyrum species should be tested by different lines of evidence, including GMM, DNA barcoding, population genetics, and ecological niche modeling.
In this study, we extended our investigation to 28 wild populations, then conducted a geometric morphometric analysis on leaf form and quantified genetic and ecological niche divergence of the two morphospecies. Our specific objectives are: (1) to investigate the morphological variation of Asteropyrum across its distribution range; (2) to test if the two morphospecies of Asteropyrum are independent by means of GMM, DNA barcoding, population genetics, and ecological niche modeling.

Field Investigation and Sample Collection
We searched for distribution information of Asteropyrum from Chinese Virtual Herbarium (CVH) and websites of nature reserves or national parks. Most recorded locations have been investigated during the last 3 years. In total, fresh leaves of 28 populations were immediately dried with silica gel (Supplementary Table 1). Sampled individuals were spaced by >10 m to avoid sampling the ramets of the same genet. Geographical information (latitude, longitude, and altitude) was recorded with a smart phone's GPS.
Four to six specimens in each population, and three mature leaves per specimen were sampled for morphological analyses. These mature leaves, each is the largest one from a randomly selected individual, were selected for avoiding developmental plasticity. Because the distinction between the two morphospecies is obscure, species identity was predefined with blind test by three experienced taxonomists according to the majority rule. All voucher specimens were stored at Jiangxi Agriculture University (JXAU).

Morphological Analysis
A total of 437 dry leaves from 27 populations were collected for morphological studies. The PT population was not measured owing to lack of samples and specimen. Leaves were pressed, dried and scanned with the adaxial uppermost using a scanner (Brother DCP-B7530DN) with a resolution of 300 dpi. Scanned images were used to record 15 landmarks (Figure 2 and Table 1 for detailed landmark definitions) following Viscosi and Cardini (2011). Cartesian coordinate data of 15 landmarks were acquired using software from the TPS Series (Rohlf, 2010) and then imported into the free software MorphoJ v. 1.04a (Klingenberg, 2011). For the landmark data, a generalized Procrustes analysis (GPA) was performed to extract shape and size components of form variation (Viscosi and Cardini, 2011), calculating the Procrustes distance (Rohlf and Slice, 1990) and Centroid distance (Loy et al., 1998), respectively. Procrustes distance and centroid distance are used to test the shape variation (Holling, 1992;Mitteroecker and Gunz, 2009) and size variation (Loy et al., 1998) between groups, respectively. For Procrustes distance, Kolmogorov-Smirnov normal test (K-S test), MANOVA test (multivariate analysis of variance) at intraspecies and interspecies levels were carried out (Hair et al., 1998;Limsopatham et al., 2018;Hošková et al., 2021). For centroid distance, K-S normal test, ANOVA test (univariate analysis of variance) were carried at the same levels as well. SPSS 22.0 software was used for these statistical analyses.
Using the variance-covariance matrix of the GPA shape coordinates, we carried out a principal component analysis (PCA), a canonical variate analysis (CVA) and a discriminant analysis (DA) for leaf shape. All of the analyses were accomplished in software MorphoJ v. 1.04a (Klingenberg, 2011).

Molecular Procedure and Species Delimitation Based on DNA Sequences
Total genomic DNA was extracted from silica-dried leaves using a modified cetyltrimethylammonium bromide (CTAB) method (Doyle and Doyle, 1987). We amplified and sequenced DNA barcoding markers (intergenic spacer psbA-trnH and internal transcribed spacer 2, ITS2) across all sampled Asteropyrum populations. To complement the barcoding analyses, five lowcopy nuclear loci were also amplified to evaluate the genetic divergence at nuclear genome, with A. cavaleriei and A. peltatum represented by nine and four populations, respectively. The primers used for the polymerase chain reaction (PCR) followed Sang et al. (1997) for psbA-trnH and White et al. (1990) for ITS2, and the primers of five nuclear loci were designed from transcriptome sequences following our previous study (Kou et al., 2020). Amplification reactions were carried out in a volume of 20 µl, containing 10 µl 2 × Taq PCR MasterMix (Tiangen, Shanghai, China), 1 µl each forward and reverse primer (0.2 µM), 1 µl template DNA (ca. 50-100 ng) and 7 µl ddH 2 O. Amplification was carried out in a Bioer XP cycler (Bioer, Hangzhou, China) programmed for an initial step 5 min at 94 • C, followed by 36 cycles of 94 • C for 50 s, 50-53 • C for 50 s and 72 • C for 1 min or 1 min and 40 s, with a final extension for 10 min at 72 • C. Sanger sequencing reactions were conducted with the corresponding forward and reverse primers commercially by Sangon Biotech Co., Ltd. (Shanghai, China).
The best substitution model was respectively, selected using jModeltest 2 (Darriba et al., 2012) for psbA-trnH and ITS2 haplotype matrices and phylogenetic inferences were carried out using PhyML v3.0 (Guindon et al., 2010) with 1,000 bootstrap pseudo replicates for estimation of the branch support. Bayesian phylogenetic analysis using MrBayes 3.2 (Ronquist et al., 2012) were also performed (mcmcp ngen = 1000000, burnin = 250). The resulting maximum likelihood (ML) and Bayesian trees were visualized with FigTree v1.4.4 (Rambaut, 2018). To define species on the basis of psbA-trnH and ITS2 haplotypes, the Automatic Barcode Gap Discovery (ABGD) species delimitation tool was employed (Puillandre et al., 2012). Briefly, ABGD is an automatic procedure that sorts the sequences into hypothetical species based on the barcode gap. It uses a range of prior intraspecific divergence to infer from the data a model-based one-sided confidence limit for intraspecific divergence. The method then detects the barcode gap as the first significant gap beyond this limit and uses it to partition the data. Inference of the limit and gap detection are then recursively applied to previously obtained groups to get finer partitions until there is no further partitioning (Puillandre et al., 2012). We used MEGA v5.05 (Tamura et al., 2011) to compute pairwise sequence comparison matrices using the p-distance, Kimura-2-parameter (K2P) corrected distance (D) (Kimura, 1980); the matrices were then analyzed in the ABGD online species delineation tool 1 .
Pairwise differences between different populations and haplotypes were calculated by software Arlequin v3.0 (Excoffier et al., 2005). These were used for mantel test with GenALEX v6.5 (Peakall and Smouse, 2012) and to construct a minimanspanning network by Network 10.3 2 , respectively. Structure 2.3.4 (Hubisz et al., 2009) was used to assess population structure with the admixture model and the assumption of correlated allele frequencies using the dataset of five low-copy nuclear loci. The number of clusters, K, corresponding to the number of populations, was explored using 20 independent runs per K. Burn-in was set to 20000 and Markov chain Monte Carlo (MCMC) run length to 200000. We ran Structure with K varying from 1 to 10, with 10 runs for each K value. The most likely number of clusters was estimated using LnP (D) (Pritchard et al., 2000) and K statistics (Evanno et al., 2005). The population clusters were visualized using the program Distruct 1.1 (Rosenberg, 2004).

Ecological Niche Modeling
According to the geographical coordinate of the presence sites, the corresponding 19 bioclimatic variables were got through DIVA-GIS software Version 7.5.0 3 . After the standardized processing of the bioclimatic matrix data, a principal component analysis was carried out in SPSS 22.0 to test whether the morphospecies are discernible based on the climatic variable. The eigenvalue >1 was selected as the principal component extraction method. The "factor score" was saved as a variable and converted into "principal component score" for scatter plot analysis. Meanwhile, the correlation analyses were also performed among 19 bioclimatic variables and applied a correlation threshold of 0.8 to select variables for ecological niche modeling.
In order to examine niche divergence between the two morphospecies of Asteropyrum, we used Maxent v3.3.3k (Phillips et al., 2006) to generate predictive distribution models based on known occurrence records and their corresponding environmental variables (compiled from the WORLDCLIM database with a resolution of 2.5 arc-min) using Maxent (Elith et al., 2006). Niche modeling was constructed for present-day, the Last Glacial Maximum (LGM) and the Last Interglaciation (LIG) periods using 150 presence records that include the sampling sites of this studies and herbarium specimen records from Chinese Virtual Herbarium (CVH). The models were run 10 times using the parameters (convergence threshold of 10 −5 , maximum iterations of 5000 and regularization multiplier of 10) and the following user-selected features: application of a random seed, duplicate presence records removal and logistic probabilities used for the output. Eighty percent of the presence records being used for training and 20% for testing the model. The models were then projected onto a broader geographic area encompassing East Asia. The calculated AUC (the Area Under the Receiver Operating Characteristic Curve) values of the Receiver Operating Characteristics (ROC) curve were used to assess the accuracy of the simulation (Elith et al., 2006). A score above 0.9 indicates that the simulation model has high accuracy and reliable results (Fielding and Bell, 1997). Maps were generated using ArcGIS Version 9.3 4 . In addition, we estimated the degree of niche overlap between the two morphospecies with ENMTools Ver. 1.4.4 (Warren et al., 2010) based on Schoener's D (Schoener, 1968).

Geometric Morphometric Analysis
Kolmogorov-Smirnov normal test was carried out on Procrustes distance and centroid distance, respectively. The results showed that the distance data used in this study were normally distributed (P > 0.05). Multivariate analysis of variance (MANOVA) test on Procrustes distance showed very significant multivariate differences at both inter-and intraspecific levels (P < 0.001, Table 2). Differences between individuals (63.4% of explained sum of squares) were greater than between species (10.1%) and populations (26.5%, Table 2). For centroid distance, univariate analysis of variance (ANOVA) test showed that inter-and intraspecific differences were both significant as well, and difference between species (43.0%) were greater than between populations (24.6%) and individuals (32.4%, Table 2). These results are in contrast with the prediction that interspecific divergence would be larger than intraspecific divergence if there were morphological gaps between species.
Principal component analysis (PCA) on the variancecovariance matrix of the Procrustes shape coordinates showed that the first three principal components explained 68.1% of the total variance. The scatter plots of PC1 vs. PC2 showed that the two morphospecies were partially overlapped ( Figure 3A). Similarly, the first three canonical variables explained 69.8% of the total variance and the two morphospecies cannot be differentiated in CVA because their scatter plots were also partially overlapped ( Figure 3B). As in the PCA and CVA, the results of the DA indicated that the range of leaf shapes ( Figure 3C) and leaf size (data not shown) partially overlapped between the two morphospecies. The cross-validation classification showed that the accuracy of leaf shape in predicting species is better than a 86.99% random chance.
We arranged the picture of one representative leaf of each population together, the pictures visually display that the leaf of Asteropyrum varies continuously in shape and size ( Figure 1B) and there is no morphological gap that are fit for delimiting species. We also generated a distribution map of leaf form (Figure 1C), the map shows that the leaf form of Asteropyrum has no obvious geographical pattern.

Molecular Species Delimitation
We obtained sequences of the cpDNA intergenic spacer psbA-trnH from 313 individuals of 28 populations. The aligned sequence matrix was 219 bp in length, containing 15 substitutions. The substitutions defined 10 different haplotypes. ITS2 sequences were obtained from the same individuals and the aligned sequences were 260 bp in length. Forty-nine substitutions determined 32 haplotypes. All haplotype sequences are deposited in GenBank. Both haplotype diversity (H d = 0.849, 0.909 for cpDNA and ITS, respectively) and nucleotide diversity (π = 0.011, 0.030 for cpDNA and ITS, respectively) are relatively high ( Table 3). The neutrality tests including Tajima's D and Fu and Li's test showed that the genetic variabilities were not due to natural selection (Table 3).
Because psbA-trnH and ITS2 have a relatively rapid evolution rate in Ranunculaceae, it is difficult to align the psbA-trnH and ITS2 sequences between Asteropyrum and outgroups, especially when Cimicifuga (psbA-trnH phylogeny) and Megaleranthis (ITS2 phylogeny) were included. However, the topology of the haplotypes of Asteropyrum did not change irrespective of which outgroup(s) were used, thus all the outgroups were kept in further analyses without any exclusion of sequence regions. Maximum likelihood (ML) and Bayesian trees of psbA-trnH and ITS haplotypes indicated that the two Asteropyrum  morphospecies were not reciprocally monophyletic, albeit most clades did not receive strong bootstrap support ( Figure 1D). Likewise, the network analyses of psbA-trnH chlorotypes and ITS2 ribotypes also did not support the monophyly of the two morphospecies (Figure 4). Distributions of chlorotypes and ribotypes did not display obvious geographical or species-specific pattern (Figure 4).
The result of mantel test based on ITS2 sequences showed that the correlation coefficient was very low (r = 0.303, P < 0.01) between genetic distances and geographical distances. The mantel test of cpDNA data was not performed due to lack of informative sites.
Five low-copy nuclear loci were sequenced for 134 individuals from 13 populations. The total aligned length was 1,536 bp, with loci ranging from 231 to 389 bp. A total of 99 segregating sites, including two indels and 68 singleton sites after excluding the sites of significant linkage disequilibrium were used to population genetic structure analysis. The Bayesian clustering algorithm indicated that the most likely number of clusters was 2, however, the two genetic clusters did not correspond to the two morphospecies ( Figure 1C). Meanwhile, ABGD analysis showed that there was no barcoding gap within Asteropyrum as well and the identified gaps occurred between Asteropyrum and outgroups (Figures 5A,C). The recursive partition of 12 for both markers were obviously unrealistic (Figures 5B,D).

Niche Differentiation Between Two Morphospecies of Asteropyrum
According to the correlation analysis, fourteen climatic variables (except for bio 4, bio 5, bio 6, bio 10 and bio 15) were retained for niche modeling. Areas under the "Receiver Operating Characteristic (ROC) Curve" (AUC) had values >0.95 for both species, indicating good predictive model performance. The projected distributions of the two morphospecies across different periods changed very little (Figure 6). The Schoener's D values of three climatic periods between these two morphospecies were more than 0.58 (Current 0.614, LGM 0.582, LIG 0.581), reflecting a high overlap of ecological niches. Accordingly, the distribution ranges of two morphospecies were overlapped at different periods (Figure 6), indicating they could have few opportunities to speciate allopatrically. Likewise, principal components analysis (PCA) of 19 bioclimatic variables (Supplementary Table 2) revealed that the two morphospecies of Asteropyrum had similar climatic requirements (Figure 3D), although multivariate plots of PC1 and PC2 showed that A. peltatum adapts to slightly colder environments than A. cavaleriei ( Figure 3D)

High Variation of Leaf Shape and Size in Asteropyrum
We observed high variation of leaf shape and size in Asteropyrum, A. cavaleriei possessing larger and more obviously lobed leaves than A. peltatum. Leaves are photosynthetic organs with many other functions such as transpiration, dissipating heat, and pest defense (Niklas, 1988;Nicotra et al., 2011). Leaf traits, including shape and size, are the results of functional trade-offs that have been resolved in various ways by different species depending on their ecological settings (Nicotra et al., 2011). Such a high number of functional trade-offs is responsible for the tremendous morphological diversity in plant leaves, especially for dicots (Niklas, 1988). Remarkably, different leaf morphology may occur between closely related species, as within-species variants, or even in the same plant (Tsukaya, 2006;Nicotra et al., 2011).
Although high leaf form variation in Asteropyrum, we did not find any morphological gap and geographical trend of leave shape (Figures 1B,C). It seems that deeply lobed leaves that traditionally belong to A. cavaleriei, generally occur hotter environments at low elevations (below 1,700 m a.b.l.) or southern localities (Figure 1C and Supplementary Table 1). Accordingly, PCA analysis of the 19 bioclimatic variables in Asteropyrum indicates that A. cavaleriei occurs in environments with higher annual mean temperature ( Figure 3D). It has been noted that by adding lobes to a leaf, the rate of heat transfer across a leaf is greater than that of an unlobed leaf of the same area (e.g., Gurevitch and Schuepp, 1990). So, deeply lobed leaves may be selected for hotter environment to reduce the temperature of leaf blades. On the other hand, deeply lobed leaves have a lower ratio of mesophyll tissue to large, highly conductive veins, they have reduced hydraulic resistance relative to less or unlobed leaves (Sack and Holbrook, 2006). Therefore, leaf lobing could represent an effective removal of hydraulic stress-prone tissue and lobed FIGURE 4 | The distribution maps and network analyses of chlorotypes (psbA-trnH, a) and ribotypes (ITS, b). The haplotypes circled with red rectangle are shared by the two morphospecies. leaves in Asteropyrum might be an adaptation to warmer conditions. Compared to leaf shape that is highly labile and responsive to a range of biotic and abiotic factors, leaf size varies in a more straightforward way: small leaves are associated with harsh conditions such as cold, hot, dry, high light, exposed, nutrient poor and saline environments (Dkhar and Pareek, 2014). In this study, we found that populations with smaller leaves (usually belonging to A. peltatum) generally inhabit high elevations or northern localities (e.g., SNJ, PL, Figure 1C), possibly indicative of an adaptation to colder and harsher environments.
In addition to an adaptation to different environments (i.e., ecotypic variation or genetic variation), phenotypic plasticity could also contribute to the variation of leaf shape and size (Royer et al., 2008). We notice that there are both lobed and slightly lobed leaves within the same population such as LCG. Although this variation might be attributed to genotypic difference among individuals, phenotypic plasticity is more likely because we also found such variation in the same individual. In plants, well-developed plasticity of many traits is usually interpreted as an adaptive response to environmental heterogeneity as a consequence of immobility and modular growth (Sultan, 2000). Populations of Asteropyrum, predominantly occurring in riparian habitats which are highly heterogeneous in terms of environmental factors such as soil moisture, could have experienced strong balancing selection to maintain polymorphisms over long periods of time within each population (Zhu et al., 2020).

Integrative Taxonomy of Asteropyrum Based on Multiple Lines of Evidence
Plant taxonomists have long recognized the importance of leaf features for identifying taxa. Leaf characters are emphasized because floral features either illustrate little variation or are available only during the relatively short flowering season for each species (Jensen et al., 2002). In fact, for some groups of plants, e.g., Quercus, Betula, as well as Asteropyrum in this case, leaf characters are considered "the most important". However, as discussed above, the leaf shape and size of Asteropyrum are labile and vary among different populations and even within populations. This casts a doubt on the species boundaries of the two species in Asteropyrum. Based on a geometric morphometric approach, this study found that the variation of leaf shape and reconstruction and haplotype network analyses based on plastid psbA-trnH and nuclear ITS2 sequences all indicate that that there is no barcoding gap within Asteropyrum and the haplotypes of both morphospecies are not reciprocally monophyletic (Figures 1, 4, 5). In addition, population genetic assignment based on five low-copy nuclear loci indicates that the two genetic clusters do not correspond to the two morphospecies, enhancing that the two morphospecies in Asteropyrum are not real genetic entities. This situation is similarly observed between Bactrocera invadens and Bactrocera dorsalis (Diptera: Tephritidae) and then the invasive fruit fly B. invadens was recommended to be synonymized with B. dorsalis (Schutze et al., 2015).
The Schoener's D values and principal components analysis (PCA) of 19 bioclimatic variables ( Figure 3D) reveals a high degree of niche overlap between the two Asteropyrum morphospecies. The results indicate that the two morphospecies have very similar niches (niche conservatism), implying that lineage separation through niche divergence mediated by disruptive selection (ecological speciation) is unlikely within Asteropyrum (Graham et al., 2004). Given the niche conservatism of Asteropyrum, speciation within this genus would have fulfilled under a scenario of allopatry caused by a geographic barrier that consists of suboptimal environmental conditions for the species in question (e.g., deserts, mountains, or oceans) (Wiens and Graham, 2005). However, the projected distribution ranges of the two morphospecies at different historical stages (current, LGM, and LIG) exhibit a high degree of overlap, without a geographic gap that can potentially act as a physical barrier (Figure 6). These results suggest that both ongoing allopatric speciation and a secondary contact after historical allopatry within Asteropyrum are also unlikely.
It is widely accepted among taxonomists that objectively evaluating several lines of evidence within a formalized framework is the most efficient and theoretically grounded approach for defining robust species hypotheses (De Queiroz, 2007). In this study, we used different approaches (GMM, DNA barcoding, phylogenetic reconstruction, and niche modeling) to evaluate the two-species hypothesis of Asteropyrum, the results show that there are no distinct boundaries between the two morphospecies of Asteropyrum in terms of leaf shape and size, genetic data and ecological niche. Although two genetic clusters were detected by population genetic analysis, the two clusters mismatch with the two morphospecies. Taken together, the two species Asteropyrum defined by leaf morphology do not reflect the divergence pattern within the genus and the present classification should be abandoned. In the future, range-wide population genomic study and in-depth phenotypic investigations would be constructive for delineating the species boundaries and understanding the evolution and speciation within Asteropyrum.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
ZZ and DF conceived and designed the project. LL, HL, and HW collected the materials. YK and JW were responsible for the lab work. SC and WZ measured the morphological traits and performed the analysis. ZZ and SC wrote the manuscript. All authors read and approved the manuscript.

FUNDING
This study was supported by grants from the National Natural Science Foundation of China (Grant Nos. 31560064 and 31770236) and Strategic Priority Research Program of Chinese Academy of Sciences (XDA20050203).