Four firs, presenting comparable high-low elevations and distinct continents, were used in this study. Abies lasiocarpa and Abies balsamea are widely distributed in North America (Figure 2; Hosie, 1969). These two species were located in a monophyletic clade according to mtDNA and cpDNA phylogeny (Aguirre-Planter et al., 2012; Xiang et al., 2015, 2018). Abies lasiocarpa tends to grow in the montane region of New Mexico and Arizona, more extensively in Utah and Colorado (Alexander et al., 1990). It is typically found at 1100–2300 m elevation in coastal subalpine conifer forests, and grows to the alpine treeline in most of its range (Little, 1971). Abies balsamea always occurred in cold regions (around 4°C) of northeastern America and eastern Canada, and the northeastern United States, at 0–1700 m elevation in boreal and north temperate forests (Little, 1971; Frank, 1990). It often forms clumps, where branches from an older central stem have layered, developed apical dominance. Abies chensiensis, as an endangered species, is distributed in Qinling-Daba and Hengduan mountains region of central China (Cheng and Fu, 1978; International Union for Conservation of Nature [IUCN], 2021). It grows in cold, moist forests at 2100–3500 m elevation, in regions with annual precipitation of 1000–2000 mm. Abies nephrolepis is a species occurring in Northeast Asia. Its elevation ranges from 0 to 1200 m, charactered by short, cool summers, and long cold winters with much snow (Cheng and Fu, 1978). Abies lasiocarpa and Abies chensiensis are widely distributed in high-elevation regions, as Abies balsamea and Abies nephrolepis occur in much lower elevations (Figure 2).

Here, the species sampling and genomic data were derived from three published studies (Jiang et al., 2011; Wang et al., 2011; Cinget et al., 2015b). Altogether, we sampled 585 individuals from 44 natural populations, including seven A. lasiocarpa populations, 12 A. balsamea populations, ten A. chensiensis populations, and 15 A. nephrolepis populations (Supplementary Table 1). Our sampling well presented the current distribution of these species. Field sampling was strictly conducted following international rules. Tsuga canadensis and Keteleeria evelyniana were chosen as outgroups.

In total, one mtDNA fragment (nad5-4) and one cpDNA fragment (trnL-trnF) were derived from three published studies (Jiang et al., 2011; Wang et al., 2011; Cinget et al., 2015b). Even all previous population-based studies in the genus Abies were widely investigated, only these two fragments and four firs could be recognised. To be mentioned, even though only two makers were used, they were almost the most informative sites available in the phylogeographical analyses of the genus Abies so far (Cinget et al., 2015a,b; Shao and Xiang, 2015; Shao et al., 2017, 2020). In these studies, more than then mtDNA (e.g., nad1-2, nad5-1, nad5-4, nad7-1, and cox1) and cpDNA (e.g., trnS-trnF, trnL-trnF, and trnS-trnG) had been surveyed, nad5-4 and trnL-trnF were the most informatic markers. Thus the markers used could be sufficient to test the relevant phylogeographical history. The relevant Genbank numbers were recorded in Supplementary Tables 2, 3.

We used DNASP version 4.10 to calculate the genetic diversity values, including h (genotype diversity) and π (differences between genetic sequences) (Nei, 1987; Rozas et al., 2003). We run the DNASP for each species and population. We used the molecular diversity values to explore the spatial distribution patterns by using GDivPAL (Dellicour and Mardulyn, 2014). Such analyses in Spads were not performed using mtDNA, because of its extremely low diversity in each species. By using the software NETWORK, we further inferred the network relationships within mitotypes and chlorotypes (Bandelt et al., 1999). Moreover, each gap was treated as one mutation.

To survey the population differentiation, we used the PERMUT to calculate N_ST and G_ST, the genetic diversity within populations (H_S), and the total gene diversity (H_T) (Pons and Petit, 1996). A phylogeographic structure could be indicated if G_ST was significantly lower than the N_ST. Using the program ARLEQUIN 3.5, we run the AMOVA (molecular variance analysis) to survey the genetical variation (Excoffier and Lischer, 2010). We performed the Mantel test to estimate the correlation between genetic distance and geographical distances (Mantel, 1967; Nei and Li, 1979; Dupanloup et al., 2002).

To detect historical distribution range expansions, the mismatch distribution analyses were performed in Arlequin 3.5. If recent expansions were indicated, further analyses could estimate potential periods. The estimated formula was s = 2 ut (t means expansion age), u = μkg (μ means the substitution rates of each nucleotide site every year, k means the average length of nucleotide sequences, and g means the potential generation years in Abies) (Rogers and Harpending, 1992). Here, we designed the value of μ as 2.61–4.02 × 10^–10 ss^–1 y^–1, the value of k as 1111 bp, according to numerous published studies in Pinaceae and Abies (e.g., Gernandt et al., 2008; Shao et al., 2020). The g value was designed as 25 years (option 1) or 50 years (option 2) for the occurrence of two different effective generation times of Abies (Xiang, 2001; Brown et al., 2004; Fan, 2006). Furthermore, such mismatch distribution analyses were performed in three species except for A. nephrolepis, because the chlorotype diversity within all the populations of A. nephrolepis was relatively low, which was fixed by C8 (Figure 3B).

The occurrence data for the whole of geographic distribution range were calculated from field observations, the Global Biodiversity Information Facility,¹ and the Chinese Virtual Herbarium.² Although many herbarium records did not specify latitude and longitude, we were able to geo-reference detailed locations, using only locations specified in the label information. In total, after removing duplicate records within each pixel (30 arc s, ca. 1 km), we obtained 267 presence records (see Supplementary Table 4 for the list of all occurrences). These records involved 53 A. lasiocarpa records, 107 A. balsamea records, 64 A. chensiensis records, and 43 A. nephrolepis records.

The species distribution modelling (SDM) of these four firs was modelled using the MAXENT (version 3.2.1; Phillips et al., 2006). The environmental data were provided by the WorldClim database³ (Fick and Hijmans, 2017), including 19 BIOCLIM variables and three different periods in 30 arc s resolution. These different periods are current (for the period 1,950∼2,000) at a spatial resolution of 30 arc s, the last glacial maximum (LGM; ∼21,000 years BP) with a spatial resolution of 2.5 arc min resolution simulated under the CCSM and MIROC models, and the last interglacial period (LIG; ∼120,000–140,000 years BP) with spatial resolution in 30 arc s resolution. Following Zurell et al. (2020) and Nualart et al. (2021), we removed those highly correlated BIOCLIM variables (pairwise Pearson correlations, r greater than 0.8), seven variables were retained for subsequent analysis (bio2: mean diurnal range of temperature, bio4: temperature seasonality, bio5: the max temperature of warmest month, bio9: mean temperature of driest quarter, bio11: mean temperature of coldest quarter, bio16: precipitation of warmest quarter, bio17: precipitation of driest quarter). The pairwise Pearson correlations among 19 BIOCLIM variables were illustrated in Supplementary Table 5. Finally, the SDMs of these four firs under different periods were conducted in 100 replicate runs, maximum iterations (1500). To ensure the consistency of the model prediction, 70% of the localities were used to train the model and 30% were used to test it. The area under the curve (AUC) value was estimated to evaluate the model performance (Phillips et al., 2006). All predictions were visualised in ArcGIS v. 10.2 (ESRI, Redlands, CA, United States). The suitable areas (in km²) (l logistic value > 0.25) for all models at each time slice were also calculated in ArcGIS.

The mtDNA nad5-4 was 246 bp in length. Four distinct mitotypes resulted from three point mutations and two indels (M1-4) (Supplementary Table 2). These four mitotypes were fixed in single species. Based on the NETWORK results, M1 and M2 were probably the ancestral haplotypes (Figure 3A). M1 and M2 were fixed in the East Asia species (A. chensiensis and A. nephrolepis) (Figure 3A). M3 and M4, which derived from M1 and M2, were exclusive to A. balsamea (Pop. 1–7) and A. lasiocarpa, respectively (Pop. 8–19).

The values of π and h were recorded in Supplementary Table 1. In AMOVA analyses, 95.17% mtDNA variation was scattered among species, and the H_T was much higher than H_S (Tables 1, 2). A significant phylogeographic structure of isolated geographic distributions were indicated (G_ST < N_ST; P < 0.05), and the geographic distances were remarkably related to mtDNA genetic diversity (r = 0.201, p = 0.037) (Table 2). In the SAMOVA analyses, there was no significant geographic division.

In total, fourteen chlorotypes were resulted by nine substitutions and three indels (C1–14) (Supplementary Table 3). Among them, four chlorotypes had a frequency above 5% (C1, C8, C12, and C14), while twelve chlorotypes (C1–10 and C13–14) were fixed in single taxa. Moreover, C11 and C12 are shared two North American species (A. lasiocarpa and A. balsamea) (Figure 3B). Based on the NETWORK analysis, C1 and C8 might be the ancestral haplotypes (Figure 3B). In A. lasiocarpa and A. chensiensis, only one or two step mutations lain in the “star-like” network (Supplementary Figure 1).

In AMOVA analyses, 94.43% of the mtDNA variation was located among species, and the H_T was much higher than H_S (Tables 1, 2). A significant phylogeographic structure isolated geographic distributions were also indicated (G_ST < N_ST; P < 0.05). In our SAMOVA analyses, there was no significant geographic division, but the results of Mantel analyses indicated remarkable correlations between geographic distances and genetical distribution (r = 0.226, P = 0.010) (Table 2).

Compared to the mtDNA, higher genetic diversity appeared in the cpDNA genealogy (Figure 3 and Supplementary Table 1). Based on the Spads analyses on the species level, the mtDNA diversity did not show any latitudinal trend (Supplementary Table 1 and Figure 4). Based on π and h diversity, only one diversity centre was recognised in two low-elevation species (A. balsamea and A. nephrolepis), as three or two diversity centres were revealed in high-elevation species (A. lasiocarpa and A. chensiensis) (Figure 4).

The analyses of mismatch distribution were performed on species level using cpDNA sequences. In total, three species were included in the comments (A. lasiocarpa, A. balsamea, and A. chensiensis), excluded A. nephrolepis because of the low cpDNA variation (Variation of cpDNA = 0; Table 1).

Our mismatch distribution results did not support the historical expansion hypothesis of A. balsamea (P_SSD = 0.020) (Supplementary Figure 2). Unlike A. balsamea, the expansion hypothesis was supported in two high-elevation firs, including A. lasiocarpa (P_SSD = 0.147) and A. chensiensis (P_SSD = 0.734) (Supplementary Figure 2). Further, our results estimated the relevant expansion periods of A. lasiocarpa in 43.8–28.4 kya (option 1: g = 25 years) or 21.9–14.2 kya (option 2: g = 50 years) (Figure 5 and Table 3). Moreover, the historical expansion of A. chensiensis probably occurred in 53.1–34.5 kya (option 1) or 26.6–17.2 kya (option 2).

The SDMs were built for four firs, separately, based on seven environmental variables under the present climate, the LGM (CCSM), the LGM (MIROC) and the LIG, respectively. All SDMs had high AUC values (greater than 0.9), indicating that the model predictions were much better than the chance (AUC = 0.5). The area of potential habitats, after accounting for uncertainty, was calculated using the logistic value > 0.25 (Table 4).

In North America, models for A. lasiocarpa showed a heavily decrease of suitable habitats while global warming, from the LGM (CCSM) (6,568,140 km²) to the present (4,710,507 km²) to the LIG (1,158,792 km²) (Figure 6A and Table 4). The potential habitats for the LGM (MIROC) (4,872,627 km²) were almost the same as the habitats under the current climate. These models predicted the lowest distribution areas under the warmest climate during the LIG. Models for A. balsamea showed obviously range expansion when global warming, from the LGM (CCSM) (172,001 km²) and LGM (MIROC) (1,789,170 km²) compared to the LIG (6,672,244) and present (5,995,647 km²) (Figure 6B and Table 4). In addition, the suitable habitat of A. balsamea during the LGM was distributed more southerly than the LIG and present models.

Compared to A. lasiocarpa and A. balsamea in North America, the suitable habitat shifting within A. chensiensis and A. nephrolepis in East Asia were unremarkable. For A. chensiensis, the models showed a slightly decreasing of suitable habitat while global warming, from the LGM (CCSM) (1,868,430 km²) to the present (1,823,649 km²) to the LIG (1,670,956 km²) (Figure 6C and Table 4). However, the LGM (MIROC) model (1,216,388 km²) was relatively smaller than others. And its potential habitats shifted along latitude during different periods. For A. nephrolepis, the models showed a slightly expansion of suitable habitat while global warming, from the LGM (CCSM) (2,324,014 km²) to the present (2,363,103 km²) to the LIG (2,537,519 km²) (Figure 6D and Table 4). But the LGM (MIROC) model (1,513,332 km²) was also relatively smaller than others. Moreover, the LGM (MIROC) models were quite different from the LGM (CCSM) models in three firs except for A. balsamea. Unlike the LGM (MIROC) models, the LGM (CCSM) models were consistent with the expansion or contraction trends from the LIG to the present models.

Previous studies concerning arctic-alpine species in Northern Hemisphere indicate two distinct expansion hypotheses in response to global climate change. These two proposed hypotheses have never been simultaneously documented in comparable high-low elevations or distinct continents (Moritz and Agudo, 2013; Theodoridis et al., 2017; Lin et al., 2021). Using four closely related firs, this study evaluated their significance on distinct elevations and continents for the first time.

Our results certainly supported the glacial expansion of high-elevation firs and interglacial expansion of low-elevation firs, either in North America or East Asia. Firstly, recent expansions were revealed in two high-elevation firs (A. lasiocarpa, P_SSD = 0.147; A. chensiensis, P_SSD = 0.734), but not in A. balsamea (P_SSD = 0.020) (Supplementary Figure 2). Similar analyses were not performed in A. nephrolepis because of its extremely low cpDNA variation (Table 1). Such recent expansions were confirmed by the Networks of chlorotypes. The chlorotypes of A. lasiocarpa and A. chensiensis showed a “star-like” pattern, which could be interpreted as characteristic of recent population expansions (Supplementary Figure 1; Hwang et al., 2003; Shao et al., 2020). Furthermore, the recent expansions’ age of A. lasiocarpa and A. chensiensis were estimated to be around 43.8–28.4 kya and 53.1–34.5 kya (option 1: g = 25 years), or 21.9–14.2 and 26.6–17.2 kya (option 2: g = 50 years), respectively (Table 3). It should be noted that the above time durations were highly consistent with the last glaciation of East Asia and North America (Figure 5). In East Asia and North America, the last glaciation was believed to occur around 73–10.4 kya, accompanied by maximum glacial around 21–18 kya (Hewitt, 1996, 2000, 2004; Yi et al., 2004). To be mentioned, potential genetic variation within continuous range species (A. balsamea and A. nephrolepis) might be difficult to capture for the generally low-variable in plastid markers. However, even though only two makers were used, they were almost the most informative sites available in the phylogeographical analyses of the genus Abies so far (Cinget et al., 2015a,b; Shao and Xiang, 2015; Shao et al., 2017, 2020). Similar analyses had all been successfully conducted in previous related studies (Cun and Wang, 2010; Shao and Xiang, 2015; Shao et al., 2017, 2020).

Our SDM models also depicted historical glacial expansions in high-elevation firs and interglacial expansions in low-elevation firs. The LGM (MIROC) models were quite different from the LGM (CCSM) models in three firs except for A. balsamea (Figure 6 and Table 4). Till now, the incongruence between the LGM (MIROC) and LGM (CCSM) had been widely surveyed (Guevara et al., 2019). The MIROC estimated noticeably more different climatic conditions than CCSM and suggested a remarkable different prediction of potential distribution (Rico et al., 2021). Compared to the LGM (MIROC), the LGM (CCSM) model was much more effective and widely used in related Pinaceae species, such as the East Asian white pines (Liu et al., 2022). In this study, the LGM (CCSM) models were highly consistent with the expansion or contraction trends from the LIG to the present models. For A. lasiocarpa and A. chensiensis, the distribution area estimations showed that the distribution range was more extensive in glacial period (the LGM CCSM) than in interglacial periods (the LIG and present) (Figures 6A,C). For A. balsamea and A. nephrolepis, the distribution range was reduced in glacial period (the LGM CCSM) than in interglacial periods (the LIG and present) (Figures 6B,D). These results suggested that A. lasiocarpa and A. chensiensis underwent expansions from the LIG till the LGM during global cooling and extensive reconstructions from the LGM to the present. Such patterns were also determined in the future geographic distribution modelling of four endemic fir species in Japan (Tanaka et al., 2012). In the context of global warming, three Japanese higher-elevation distributed firs (Abies mariesii, Abies veitchii, and Abies homolepis) were predicted to lose 45% suitable habitat but greater than 73% for lower-elevation distributed firs (Abies firma). Moreover, unlike North America, the habitat shifting in East Asia were unremarkable (Figure 6 and Table 4). For A. chensiensis and A. nephrolepis, the models showed slightly expansion or contraction of suitable habitats during different time periods. Such distinct spatial patterns could be explained by the difference in Quaternary climate between East Asia and North American (Shi et al., 1987; Shi, 2002). During Quaternary climatic fluctuations, North America has been repeatedly covered by large continuous ice sheets during glaciations, and only a few montane glaciers or valley glaciers occurred in East Asia regions. Our models indicated that the Asian range of the genus Abies was characterised by much more stable conditions, which might have facilitated the survival of many woody species that represent genera with highly disjunct distribution (e.g., Aesculus, Pterocarya, Fagus, Castanea) survived in East Asia (Zhang et al., 2020; Song et al., 2021).

The geographical distribution of genetic diversity could be another direct independent evidence. To distinguish interglacial and glacial expansion hypotheses, Hewitt (1996) proposed two different expansions models. One is the model of pioneer expansions, charactered by interglacial expansion hypothesis, which predicted along-latitude decreasing of genetic diversity while historical expansion (Hewitt, 2004). Another is the model of phalanx expansions, charactered by glacial expansion hypothesis, indicated along-elevation historical expansions, resulting in multiple diversity centres and no latitude trends (Oliveira et al., 2021; Wyatt et al., 2021). Based on current phylogeographical patterns, multiple diversity centers were revealed in only A. lasiocarpa and A. chensiensis, indicating no latitude decreasing trends (Figure 4). Thus the phalanx expansion model was probably activated in these two high- elevation firs. It should be noted that such pattern could also result from the fact that the isolation between potential refugia in the lowlands was less than in the mountains, which may have resulted in homogenization after the expansion period (Liu et al., 2022).

Moreover, the palynological records could be another convincing proof of species’ preferred past climatic conditions. In North America, the southernmost paleo records of A. lasiocarpa were US Rockies and Oregon, and could be its putative glacial refugia (Barnosky, 1981). Based on published macrofossils, A. lasiocarpa was originated in the Pleistocene and much common in the lower elevation regions, including Vancouver Island and Washington state (Bartlein et al., 1998; Cinget et al., 2015b). Nowadays, those lower elevation regions are dominated by Abies amabilis and Abies grandis (Farjon, 2001). For Abies balsamea, the palynological records reflected post-glacial migration from the refugium in southern and western Great Lakes towards central Alberta (Cinget et al., 2015a). In East Asia, numerous published palynological researches inferred a similar pattern (Xu et al., 1980, 2002; Cook et al., 2011). During interglacial periods, the pollen abundance in the high-elevation regions was obviously decreased in Eastern Asia (Tang et al., 2007; Yong et al., 2020). Indeed, those high-elevation fir species in central China (e.g., A. chensiensis) distributed more commonly in lowland forests during the LGM as the temperature decreased about 6–10°C (e.g., Li et al., 2019). Zhao et al. (2014) reported several remarkable impacts of repeating climate fluctuations since the LGM on the distribution of Abies in the mountaintop regions of central China, with obviously decreasing continuously with climate change from humid and cold (glacial) to dry and warm (interglacial). In addition, almost all the critically endangered Abies worldwide are restricted in the mountain top (Farjon, 2001; Xiang, 2001). Such current distribution pattern could be living evidence of interglacial fragmentation in high-elevation firs. Even though the significant deforestation during past decades, these critically endangered Abies species had never been the main targets for their high-elevation distribution in mountain tops and far away from human settlements (Fan, 2006). Thus climate-driven distribution changes and low seed viability could be the most critical factors.

In summary, the above results of five independent evidence lines, involving geographic-driven genetic patterns, mismatch distribution analysis, species distribution modelling, palynological records, and current distribution pattern, were highly uniform and reliable. Our findings supported the glacial expansion hypotheses of high-elevation firs and interglacial expansion hypotheses latitudinal of low-elevation firs during the Quaternary climate change in North America and East Asia.

The genus Abies Mill. (Pinaceae) discontinuous distributed in North America, Eastern Asia, and Europe (Xiang et al., 2018). Fir species are ecologically, and economically important plants that are dominant in Taiga and mountain conifer forests which maintain water and soil, protect the ecological environment, and conserve the biodiversity in the region (Farjon and Rushforth, 1989). However, eleven firs had been listed as endangered or critically endangered species (Xiang, 2001; International Union for Conservation of Nature [IUCN], 2021). Thus it is urgent to develop efficient conservation planning.

The primary goal in conservation biology is to determine how species respond to climate change, which generated conservation units to define (Petit et al., 1998). At present, all the endangered fir species had an extremely high elevation habitat (>1000 m), which was quite different from the two low-elevation firs used in this study (A. balsamea: 0–1700 m, A. nephrolepis: 0–1200 m) (Farjon, 2001). Furthermore, six fir species were listed as critically endangered (Abies beshanzuensis, Abies yuanbaoshanensis, Abies fanjingshanensis, Abies ziyuanensis, Abies fansipanensis, and Abies numidica) (Fu, 1992; Xiang, 2001; International Union for Conservation of Nature [IUCN], 2021). Their distribution ranges were restricted in limited mountain tops and extremely low effective breeding size. The critically endangered status might be resulted by historical rapid recolonization via elevational shifting in situ, following severe bottleneck effects (Cinget et al., 2015a,b; Jiang et al., 2015; Shao et al., 2020). Such evolutionary progress should not be constrained. For high-elevation firs, with A. lasiocarpa and A. chensiensis as example, more than one diversity centers were depicted (Figure 4). These different diversity centers and genotypes drawn from different populations cannot be ignored. For conservation purposes, it is necessary to define conservation units on an evolutionary level for these endangered firs.

In the context of global warming, our results showed that the high-elevation firs in Northern Hemisphere had undergone glacial expansion via elevational shift, which led to massive habitat fragmentation or isolation. Under this situation, the suitable habitat of high-elevation Abies could have begun to lose continually since the LGM and quicken steps, as depicted in the SDMs of A. lasiocarpa and A. chensiensis (Figure 6). Suitable habitats of these endangered Abies species continued to deteriorate, and the surviving individual numbers were continuously decreased (Xiang, 2001; Peng et al., 2012; International Union for Conservation of Nature [IUCN], 2021). At present, these critically endangered firs were vulnerable to the verge of extinction, for no enough elevational space to up-migrate in respond to future climate warming (Habel et al., 2011). By analysing available studies of latitudinal and elevation range shifts for multiple taxonomic groups, most species’ range shifting rate, in responding to climate change, were located around eleven meters along elevation per decade or 16.9 kilometres along latitudes per decade (Chen et al., 2011). As stated above, the endangered Abies species would be forced to contract up the slope by elevation shifting at a relatively high speed during the foreseeable global warming. However, these fir species have almost reached the mountaintop while no enough space to continue up-migrate, for example, Abies beshanzuenesis located in Baishanzu Mountain (1857 m) at 1750 m, Abies fanjingshanensis restricted in Fanjing Mountain (2494 m) with altitudes from 2100 to 2350 m, Abies yuanbaoshanensis distributes in Yuanbao Mountain (2081 m) at 1900–2000 m, and Abies ziyuanensis occurred in Shunhuang Mountain (1882 m) at 1700 m (Fan, 2006). To protect these critically threatened Abies, their up-migrate along elevation must be persisted and demanded critical attention by the government managements and protection organizations.

Till now, the general government protections of these critically endangered firs were fixed by in-situ conservation, but it was found to be not an effective and practical program. Indeed, the continued up-migrate of the six highly endangered firs was unsustainable. Facing the ongoing climate warming, they were believed to be near extinction or heavily contracted (Xiang, 2001; Xiang et al., 2009). Thus, we propose that the most helpful and urgent strategy was ex-situ management in more suitable areas, accompanied by breeding programs in botanical gardens. Even though these critically endangered firs were not included in our SDMs, the tendency was still reliable and they could be well represented for their conservative ecological niche requirement (Shao et al., 2017). As for the five critically endangered firs in China and Vietnam, the southern Hengduan Mountains region was always the suitable and ideal habitat during the whole Quaternary according to the SDMs results (Figure 6). And more than 35 National Botanical Reserves were constructed in the Heng-duan Mountains region and could play a meaningful role in the ex-situ conservation (Chen et al., 2009; Shao et al., 2020). Similar specific gardens and existing programs in other countries should be activated and strengthened the ex-situ conservation of defining conservation units. Our results are timely and significant to biodiversity conservation for endangered firs at current and future global climate change.