To collect the occurrence records of F. cirrhosa, F. unibracteata, and F. przewalskii, respectively, we conducted an extensive literature search using online databases and referred to the Chinese Virtual Herbarium database (CVH)¹ and the plant photo bank of China (PPBC),² global biodiversity information facility (GBIF),³ and field investigation as well. Then, records lacking latitude and longitude were georeferenced using Google maps module that is included in the Chinese Satellite Map.⁴ Using the above resources, the distributional localities of F. cirrhosa, F. unibracteata, and F. przewalskii were added to a database. After deleting the duplicated records and filtering spatially by Buffer analysis, only one point was left within each grid cell (buffer radius is 20 km). Furthermore, to avoid sampling bias, spatial autocorrelation was measured by calculating Moran‘s I using GeoDa v1.14 software (Wu et al., 2021). Finally, the remaining records were used for constructing the models. With the help of ArcGIS 10.2, the sighting point map was obtained (Figure 2).

Initially, 22 environmental variables (Supplementary Table 1) that might influence the distribution of F. unibracteata, F. cirrhosa, and F. przewalskii were selected to model the current species distribution patterns. These variables included 19 bioclimatic variables that described annual and seasonal temperature and rainfall trends (WorldClim version 2.0, Fick and Hijmans, 2017) and three topographic variables (slope, altitude, and aspect extracted from elevation data through the spatial analysis module in ArcMap 10.2.). The 19 WorldClim variables have been broadly used for generating ENMs, including successful use with plants in the Himalayan–Hengduan Mountains (Ma and Sun, 2018; Zheng et al., 2021). These three topographic variables have been indicated to play a role in the distribution and density of Fritillaria species (Lakey and Dorji, 2016).

For future climate scenarios, we used BCC-CSM1.1 climate change modeling data underlying the shared socioeconomic pathways SSP126, SSP245, SSP370, and SSP585 scenarios released by the CMIP6 and IPCC Assessment Report 6 (AR6). BCC-CSM 1.1 is among the most-used models currently available for simulating the global climate response to increasing greenhouse gas concentration. Our four selected future climate data sets, including 2021–2040, 2041–2060, 2061–2080, and 2081–2100, were downloaded from the World Climate Database (Abou-Shaara et al., 2021) with a spatial resolution of 2.5 arc-min. All the environmental variables were in raster format and were prepared in ArcGIS 10.2 to align in the geographic space using a WGS84 datum system with default parameters.

As the high collinearity of 22 environmental variables may lead to the model overfitting, the Spearman correlation analysis was used to test multicollinearity among these variables. Before the correlation analysis, the environmental variables with a contribution rate <0.3 were removed. According to the present distribution points of Fritillaria, the associated environmental variables were extracted and with those values, a correlation matrix in SPSS 20 software was performed. Then, following the methods of Jiang (2017) and Wang et al. (2018), we identified environmental variables that had a coefficient of correlation >0.8 with other environmental variables and eliminated the one with a smaller contribution rate (Supplementary Table 2). Thus, the remaining environmental variables were used in the further prediction of models (Feng et al., 2019).

We used the MaxEnt (MaxEnt version 10.2) and GARP models. For each species, various training/testing sets (80/20, 75/25, and 70/30) were used for validating the models. The algorithm runs 500 iterations of these processes. We selected test samples by the bootstrap method. The Jackknife method was checked in the setting of environmental parameters and the response curve of environmental variables was created. The remaining settings were left at the default setting; finally, the test was repeated 10 times. The outputs were transformed into raster format using the ArcMap tool in ArcGIS software for further analysis.

The effectiveness of MaxEnt and GARP models was evaluated by three standard performance criteria, including area under the curve (AUC) of the receiver operating characteristic (ROC) curve, true skill statistics (TSS), and Kappa values. AUC of the ROC curve is the threshold-independent measure of model accuracy, which juxtaposes correct and incorrect predictions over a range of thresholds (Zhu et al., 2012). The AUC value was calculated as the area enclosed by the ROC curve and generally ranges between 0 and 1. A larger AUC value indicates better model precision. The AUC values could be grouped as failing (AUC ≤ 0.6), poor (0.6–0.7), moderate (0.7–0.8), good (0.8–0.9), and excellent (0.9–1). TSS scores range from -1 to 1, where +1 indicates a perfect ability to distinguish suitable from unsuitable habitat, while values of zero or less indicate a performance no better than random (Assefa et al., 2021). Kappa (Schatz et al., 2017; Fernández and Morales, 2019) measures how closely model predictions match the truth, controlling for random accuracy by estimating the expected accuracy. Kappa values range from -1 to 1, where 1 indicates 100% agreement, and -1 indicates 100% disagreement.

We imported prediction results into ArcGIS and then drew global maps of ecologically appropriate zones for three Fritillaria species, respectively. During the ArcGIS mapping, artificial grading was used to classify different grades based on their ecological similarity. The final potential species distribution map had a range of values from 0 to 1 which were regrouped into four classes of suitable habitats, including highly suitable habitat (0.75–1), moderately suitable habitat (0.5–0.75), low suitable habitat (0.25–0.5), and not suitable habitat (0–0.25) (Ren et al., 2020). These three levels of suitable habitats (25–100%) indicated that species have a probability of survival in these areas (Swets, 1998). The suitable ranges of the environmental variables were defined as those ranges where species inhabited survival areas with low, moderate, and high suitability (≥0.25) (Zhang et al., 2018). The range of appropriate eco-factor values was derived from the response curves of the MaxEnt model results.

The niche breadth in geographical and environmental space was calculated for each species as the mean Levins’ B1 (inverse concentration) and B2 (uncertainty) values (Levins, 1968) from habitat suitability maps of each species with the software ENMTools v1.3 (Warren et al., 2008). Levins’ B1 and B2 values range from 0 to 1, with values closer to 0 representing narrow niche breadth, and values closer to 1 representing wide niche breadth. The ENMTools package was further used to evaluate niche overlap based on the values of two indexes, Schoener’s D (Schoener, 1968) and Hellinger’s I (Warren et al., 2008). Schoener’s D and Hellinger’s I values range from 0 to 1, with values closer to 0 representing a small degree of niche overlap, while values closer to 1 representing a high degree of niche overlap.

To strengthen our speculation, the overlapping degree between highly and moderately suitable habitats of paired species is defined as follows, which is the proportion of the records in the overlapping region divided by their union.

The A_overlap and B_overlap represent the number of records of species A and species B in the overlapping regions, respectively. While A_total and B_total are the numbers of all existing records that are related to the respective individual species.

Internal transcribed spacer 1 (ITS1) and internal transcribed spacer 2 (ITS2) are the internal regions between 18S and 5.8S, 5.8S, and 28S, respectively. Both have been widely used as molecular markers to explore phylogenetic relationships (Cheng et al., 2016; Khan et al., 2019). While, as a promising method, phylogenetic analysis based on the Chloroplast genome (CP) was also applied in our research. We then searched the GenBank database for DNA sequences of ITS1, ITS2, and CP data, some of which were obtained by our laboratory previously (Zheng et al., 2019). Our strategy included (1) marker sequences available for all the species; (2) sequences that could be linked for all species; (3) sequences contributing to extensive gaps in the conserved regions of the alignment were discarded; and (4) sequences that could be assigned to a set of geographic coordinates.

The ITS1+ITS2, ITS1, ITS2, and CP trees were constructed using a maximum likelihood using the MEGA 7.0 (Deng et al., 2018). To obtain statistical support on the resulting clades, a bootstrap analysis with 1,000 replicates was performed. The evolutionary distance was defined as the number of nucleotide substitutions occurring between various sequences based on the Tamura–Nei model (Tamura and Nei, 1993). In the maximum likelihood tree, the pairwise distance was exported to calculate the evolutionary distance (Kumar et al., 2018), which was used for the further correlation analysis.

A linear trend line has been drawn through the evolutionary distances and the occurrence records in the overlapping area. The R² value indicates how well data fit the line. To test the statistical significance of the correlation results, a Pearson Product Moment Correlation test was performed using SigmaStat v3.5 software. A correlation coefficient (ρ) positive and a p-value lower than 0.05 means that the two variables tend to increase in a concerted manner (Santamaría et al., 2014).

In the initial step, 112, 51, and 73 native occurrence records of F. cirrhosa, F. unibracteata, and F. przewalskii, respectively, were obtained. After filtering data, the records of F. cirrhosa, F. unibraceate, and F. przewalskii were decreased to 77, 45, and 65, respectively (Supplementary Table 3). The corresponding Moran I of three Fritillaria species were 0.0621465, -0.01864, and 0.0016039, respectively. The procedure greatly reduced sampling bias and spatial autocorrelation, resulting in evenly distributed occurrence points across space, and thus the remaining records could be used for further model construction.

With the given training and testing sets (80/20, 75/25, and 70/30), the MaxEnt model performed better than the GARP model by evaluation metrics (AUC, KAPPA, and TSS). When training/testing set was 75/25, F. unibracteata (AUC = 0.972 (SD = 0.006), KAPPA = 0.913, and TSS = 0.890) and F. przewalskii (AUC = 0.981 (SD = 0.006), KAPPA = 0.871, and TSS = 0.913) had better performance, while 70/30 set was more suitable for F. cirrhosa [AUC = 0.969 (SD = 0.006), KAPPA = 0.731, and TSS = 0.730] (Supplementary Table 4). Therefore, MaxEnt was used as the final model for three Fritillaria species, in which F. unibracteata and F. przewalskii employed 75/25 set, while F. cirrhosa employed 70/30 set. The internal jackknife test in the MaxEnt software package verified the importance of various environmental variables on the three Fritillaria species, respectively (Figure 3). When used independently, Bio12 (23.7% of contribution), Bio11 (17.7% of contribution), Elev (16.5% of contribution), and Bio4 (15.5% of contribution) provided the greatest contributions to the distribution model for F. cirrhosa relative to other variables. The cumulative contributions of these four factors reached values as high as 73.4%, indicating that these variables contained more useful information than the other variables. Bio12 (28.2% of contribution), Elev (27.9% of contribution), Bio4 (22.9% of contribution), and Bio2 (10.6% of contribution) made the highest contributions to the distribution model of F. unibracteata. The cumulative contributions of the variables reached values as high as 89.6%. Furthermore, the environmental variables that provided the greatest contributions to F. przewalskii were Elev (44.3% of contribution), Bio12 (21.5% of contribution), and Bio4 (10.3% of contribution). Therefore, the important environmental factors that influenced the suitable habitat were identified as Bio12, Bio11, Elev, and Bio4 for F. cirrhosa, Bio12, Elev, Bio4, and Bio2 for F. unibracteata, and Elev, Bio12, and Bio4 for F. przewalskii, respectively.

The response curve of environmental variables takes the specific value of environmental variables as Abscissa and the existence probability of species as ordinate, which reflected the suitable range of environmental variables in the study area (Figure 4). The suitable ranges of the environmental variables were defined as those ranges where three Fritillaria species inhabited survival areas with low, moderate, and high suitability (≥0.25) (Zhang et al., 2018).

In F. cirrhosa, the suitable ranges of Bio12, Elev, Bio4, and Bio11 (see definitions in Table S1) were 468.68–1263.40 mm, 949.40–4631.11 m, 403.91–865.23, and -9.28 to 12.33°C, respectively. When Bio12 exceeded 468.68 mm, the probability for F. cirrhosa increased gradually, and reached the maximum at 707.59 mm with a probability of existence as high as 0.67. When Elev exceeded 949.40 m, the survival probability of F. cirrhosa first increased and then decreased, and reached 3716.54 m with the highest probability of existence as 0.78. When the Bio4 reached 532.05, F. cirrhosa obtained the highest probability of existence as 0.69. Finally, it achieved the top probability of existence as 0.65 once Bio11 reached -0.17°C (Figure 4A). In F. unibracteata, the optimum values of Bio12, Elev, and Bio4 ranged from 421.41 to 1338.27 mm, 1453.41 to 4653.12 m, and 513.42 to 888.86, respectively (Figure 4B). Furthermore, in F. przewalskii, the optimum ranges of Bio12, Elev, and Bio4 were from 335.98 to 1073.34 mm, from 1576.31 to 4342.94 m, from 537.14 to 967.75, respectively (Figure 4C). The existence probability of three Fritillaria species increased at first and then decreased with the increase of environmental factors. Table 1 shows the suitable range of environmental variables corresponding to the three Fritillaria species.

The prediction of current habitats for F. cirrhosa showed that the highly suitable habitat occurred in Sichuan, Tibet, Yunnan, Gansu, Shanxi, Henan, Qinghai, Guizhou, and Hubei provinces. The evaluated moderately suitable habitat distributed in Shanxi, Ningxia, Chongqing, and Hebei. The current areas of suitable habitats for F. cirrhosa were up to 0.68 million km², accounting for 7.10% of the total area of China (Figure 5A1).

In F. unibracteata, the evaluated highly suitable habitat was located in Sichuan, Qinghai, Gansu, and Tibet. The evaluated areas of moderately suitable habitat included Sichuan, Tibet, Qinghai, Gansu, Shaanxi, Henan, Hubei, and Yunnan. The current areas of suitable habitat for F. unibracteata encompassed 0.48 million km² (Figure 5B1). As for F. przewalskii, the highly suitable habitat occurred in Gansu, Qinghai, Sichuan, and Tibet and the moderately suitable habitat distributed around the highly suitable growth area. The total area suitable habitat for F. przewalskii was 0.35 million km², accounting for 3.64% of the total area of China (Figure 5C1).

The centroid of the current habitat of F. cirrhosa was located at the site of 100°699′E and 30°275′N in south-west Sichuan province. The centroid of the suitable area shifted to 101°130′E, 30°399′N, 101°639′E, 31°704′N, 100°200′E, 30°656′N, 100°264′E, 30°3′N under SSP126 from 2021 to 2100. Under SSP245, the centroid of the future suitable area moved to 100°235′E, 30°584′N, 100°35′E, 31°507′N, 101°315′E, 31°331′N, 102°292′E, 31°302′N from 2021 to 2100. Under SSP370, the centroid of the future suitable area changed to 100°577′E, 30°486′N, 99°121′E, 30°558′N, 100°917′E, 30°611′N, 99°827′E, 31°089′N. Moreover, on the condition of SSP585, the centroid of the suitable area shifted to 98°741′E, 30°406′N, 99°782′E, 30°505′N, 101°460′E, 31°334′N, 100°540′E, 30°859′N. Overall, the core distribution of F. cirrhosa shifted toward the north-east under four scenarios (Figure 5A2). The elevation of the centroid under four scenarios (4,143–4,601 m) was higher than that of the current centroid (3,907 m).

In F. unibracteata, the centroid of the current habitat was located at the position of 101°124′E, 33°092′N in north-west Sichuan province. The centroid of the suitable area shifted to 101°132′E, 33°468′N, 100°991′E, 32°789′N, 101°087′E, 32°745′N, 100°645′E, 33°041′N from 2021 to 2100 under SSP126. Under SSP245, the centroid of the future suitable area moved to 101°257′E, 33°018′N, 100°918′E, 33°022′N, 101°114′E, 32°407′N, 100°653′E, 32°673′N from 2021 to 2100. Under SSP370, the centroid of the future suitable area changed to 101°081′E, 33°238′N, 101°236′E, 33°189′N, 101°038′E, 32°691′N, 101°238′E, 31°838′N from 2021 to 2100. Finally, under the condition of SSP585, the centroid of the suitable area shifted to 101°150′E, 33°078′N, 101°078′E, 32°070′N, 101°192′E, 32°732′N, 101°877′E, 31°742′N from 2021 to 2100. Overall, the core distribution of F. unibracteata shifted toward the south under four scenarios (Figure 5B2).

Additionally, in F. przewalskii, the centroid of the current habitat was located at the position of 102°236′E, 33°99′N in north-west Sichuan province. The centroid of the suitable area shifted to 101°769′E, 33°627′N, 102°511′E, 34°256′N, 101°459′E, 33°536′N, 101°708′E, 33°719′N from 2021 to 2100 under SSP126. Under SSP245, the centroid of the future suitable area moved to 102°328′E, 34°331′N, 101°911′E, 33°697′N, 100°996′E, 33°39′N, 101°299′E, 33°834′N from 2021 to 2100. Under SSP370, the centroid of the future suitable area changed to 102°059′E, 34°102′N, 101°766′E, 33°711′N, 101°355′E, 33°904′N, 101°377′E, 33°246′N from 2021 to 2100. Finally, under the condition of SSP585, the centroid of the suitable area shifted to 101°559′E, 33°647′N, 101°85′E, 34°141′N, 101°036′E, 33°852′N, 101°039′E, 33°564′N from 2021 to 2100. Generally, the core distribution of F. przewalskii shifted slightly toward the south-west under four scenarios (Figure 5C2). The elevation of the centroid under four scenarios (3,503–4,284 m) also increased compared with that of the current centroid (3,425 m).

As shown in Supplementary Table 5, in F. cirrhosa (Figure 6A), the areas of highly and moderately suitable habitats, as well as suitable habitats (including high, moderate, and low suitable habitats) increased in a wave-like change under the SSP126 scenario. The maximum area of highly and moderately suitable habitats appeared in 2041–2060 with 0.48 and 2.32%, respectively. Under the SSP245 scenario, the highly and moderately suitable habitats, as well as the suitable habitat continuously increased from 2020 to 2100. Both the maximum area of highly and moderately suitable habitats appeared in 2081–2100 with 0.97 and 4.12%, respectively. Moreover, under the SSP370 scenario, the highly and moderately suitable habitats showed a similar increasing tendency as SSP126. The highly and moderately suitable habitats reached the maximum at 2061–2080 period with 1.39 and 4.01%, respectively. Finally, under the SSP585 scenario, the highly and suitable habitats as well as suitable habitat showed similar tendencies as those under the SSP126 scenario from 2020 to 2100. The highly and moderately suitable habitats reached the maximum at 2081–2100 with 1.53% and at 2061–2080 with 3.95%, respectively. Overall, SSP585 was the optimum scenario for F. cirrhosa, since the highly and totally suitable habitats reached the maximum areas during 2021–2100 at 146,758 and 133,884.2 km², respectively.

In F. unibracteata (Figure 6B), the areas of highly and moderately suitable habitats as well as suitable habitat increased first and then declined slightly (but higher than that of the current corresponding habitat) under the SSP126 scenario. The maximum area of the highly and moderately suitable habitats appeared in 2061–2080 with 0.36 and 1.66%, respectively. Under the SSP245 scenario, the highly and moderately as well as totally suitable habitat increased in a wave-like change. The maximum area of highly and moderately suitable habitats appeared in 2041–2060 with 0.46% and in 2081–2100 with 1.61%, respectively. Furthermore, under the SSP370 scenario, the highly and moderately, as well as totally suitable habitat showed a similar tendency as those under the SSP245 scenario. The highly and moderately suitable habitats reached the maximum at 2081–2100 with 0.54% and at 2061–2080 with 1.96%, respectively. Finally, under the SSP585 scenario, the highly and moderately as well as the totally suitable habitat increased continuously. Both highly and moderately suitable habitats reached the maximum at 2081–2100 with 0.36 and 2.08%, respectively. Overall, SSP370 was the optimum scenario for F. unibracteata, since the highly and totally suitable habitats reached the maximum areas during 2021–2100 at 517.39 and 831.501 km², respectively.

In F. przewalskii (Figure 6C), the areas of highly and moderately as well as totally suitable habitat increased in a wave-like change under SSP126, SSP245, SSP370, and SSP585 scenario, respectively. Under the SSP126 scenario, both the maximum area of highly and moderately suitable habitats appeared in 2061–2080 with 0.65 and 1.66%, respectively. Under the SSP245 scenario, the maximum area of both highly and moderately suitable habitats appeared in 2080–2100 with 0.89 and 2.24%, respectively. Additionally, under the SSP370 scenario, the highly and moderately suitable habitats reached the maximum at 2061–2080 with 0.84% and 2081–2100 with 2.00%, respectively. Finally, under the SSP585 scenario, the highly and moderately suitable habitats reached the maximum at 2081–2100 with 1.10 and 2.43%, respectively. SSP585 was the optimum scenario for F. przewalskii since the highly, moderately, and suitable habitats reached the maximum areas during 2021–2100 at 106,043, 232,835, and 784,581 km², respectively.

We carried out a prediction on the ecological niche overlap of three Fritillaria species based on the distribution maps (Figure 2). As a result, F. unibracteata showed the broadest niche (indicated by the highest value of B1 and B2: 0.187 and 0.908, respectively), while the F. cirrhosa showed a narrower niche (B1 and B2: 0.155 and 0.886, respectively) and F. przewalskii showed the narrowest (B1 and B2: 0.112 and 0.874). The most similar niches were of F. unibracteata and F. przewalskii, which is indicated by the highest values of niche overlap (I: 0.927, D: 0.723). Fritillaria cirrhosa was more similar to F. unibracteata (I: 0.896, D: 0.659) than to F. przewalskii (I: 0.879, D: 0.619). Then, based on the statistics and calculations, the overlapping degree in highly and moderately suitable habitats of F. cirrhosa and F. unibracteata was 34.43%, those of F. cirrhosa and F. przewalskii was 26.06% and those of F. unibracteata and F. przewalskii was 50.00% (Figure 7B). The results coincided with the results of prediction by ENMTools that F. unibracteata and F. przewalskii shared the highest niche overlap.

After predicting the overlapping and specific suitable habitats of the three Fritillaria species, we found that the densely overlapping highly and moderately suitable habitats were mainly distributed in the Himalayan–Hengduan Mountains region, and mainly concentrated in Sichuan, Tibet, Qinghai, Gansu, and Shanxi provinces (Figure 7B). The specific highly suitable habitat of F. cirrhosa was distributed at southwest Sichuan, north Yunnan, east Tibet, and some were distributed at Shaanxi, Henan, and Hubei. The specific highly suitable habitat of F. unibracteata was mainly concentrated in north Sichuan, south Qinghai, and west Gansu. The specific highly suitable habitat of F. przewalskii was mainly concentrated in Gansu, Qinghai, and Sichuan province as shown in Figure 7A.

To study the correlation between the overlapping extent and the evolutionary distance between different Fritillaria species, ITS1 (Figure 8A), ITS2 (Figure 8B), ITS1+ITS2 (Figure 8C), and CP (Figure 8D) were obtained, and further the phylogenetic trees were constructed (Figure 8) and the corresponding evolutionary distance was calculated (Supplementary Table 6). As shown in Table 2, the correlation between evolutionary distances based on ITS2 and CP, and overlapping extent showed weak accuracy (no correlation, R² = 0.289) with no significant correlation (P >> 0.05), whereas that of ITS1 obtained medium accuracy (R² = 0.766) but high P-value. The evolutionary distances based on ITS1+ITS2 and overlapping extent showed a positive correlation with relatively higher accuracy (R² = 0.93) and lower P-value (P = 0.12).