Concealed truth: Modeling reveals unique Quaternary distribution dynamics and refugia of four related endemic keystone Abies taxa on the Tibetan Plateau

Abstract Understanding the factors driving the Quaternary distribution of Abies in the Tibetan Plateau (TP) is crucial for biodiversity conservation and for predicting future anthropogenic impacts on ecosystems. Here, we collected Quaternary paleo‐, palynological, and phylogeographical records from across the TP and applied ecological niche models (ENMs) to obtain a profound understanding of the different adaptation strategies and distributional changes in Abies trees in this unique area. We identified environmental variables affecting the different historical biogeographies of four related endemic Abies taxa and rebuilt their distribution patterns over different time periods, starting from the late Pleistocene. In addition, modeling and phylogeographic results were used to predict suitable refugia for Abies forrestii, A. forrestii var. georgei, A. fargesii var. faxoniana, and A. recurvata. We supplemented the ENMs by investigating pollen records and diversity patterns of cpDNA for them. The overall reconstructed distributions of these Abies taxa were dramatically different when the late Pleistocene was compared with the present. All Abies taxa gradually receded from the south toward the north in the last glacial maximum (LGM). The outcomes showed two well‐differentiated distributions: A. fargesii var. faxoniana and A. recurvata occurred throughout the Longmen refuge, a temporary refuge for the LGM, while the other two Abies taxa were distributed throughout the Heqing refuge. Both the seasonality of precipitation and the mean temperature of the driest quarter played decisive roles in driving the distribution of A. fargesii var. faxoniana and A. recurvata, respectively; the annual temperature range was also a key variable that explained the distribution patterns of the other two Abies taxa. Different adaptation strategies of trees may thus explain the differing patterns of distribution over time at the TP revealed here for endemic Abies taxa.


| INTRODUC TI ON
The Tibetan Plateau (TP) is one of the world's major hotspots of plant species diversity, harboring 12,000 species of vascular plant, including 3,673 endemic species (Wen, Zhang, Nie, Zhong, & Sun, 2014;Yu et al., 2019;Zhang et al., 2016). From a geographical standpoint, the TP is unique (Favre et al., 2015;Liu et al., 2016): Over time, great fluctuations in temperature Shen et al., 2016) and precipitation Shen, Piao, Cong, Zhang, & Jassens, 2015) have occurred, in addition to the effects of elevation on generating different climatic conditions (Favre et al., 2015). These factors, along with the past geological history and uplifting of the TP, together explain its high level of species richness and endemism (Favre et al., 2015;McCain & Grytnes, 2010;Xing & Ree, 2017). The unique features that characterize this region have allowed Abies taxa to persist and evolve over an extended period of time (Guan & Zhou, 1983;Xiang, Cao, & Zhou, 2006). The speciation events of Abies on the TP, and related migrations of the tree genus have become a hot research topic. Recent studies indicated the most important speciation events had occurred during the late Eocene, due to regional drying that led to a geographic isolation of drought-sensitive taxa (Wang, Xu, et al., 2017;Xiang et al., , 2006 Farjon, and A. recurvata Mast.) constitute an unresolved biogeographic and paleo-botanical conundrum. These four taxa vary in their morphological traits and have distinct environmental requirements (China Forest Editorial Board, 1997;Guan & Zhou, 1983). Presently, A. fargesii var. faxoniana and A. recurvata occur exclusively in the northeastern area of the Hengduan Mountains ( Figure 1), which are located in the eastern part of the TP (Guan & Zhou, 1983). That region is characterized by a relatively warm climate with a wet summer (He, He, & Wu, 2015). By contrast, both A. forrestii var. georgei and A. forrestii are distributed only in TP's southeastern part (Figure 1), specifically in the southern portion of the Hengduan Mountains range (Fang, Wang, & Tang, 2011), featuring a cold and extended sunny climate with a dry winter (China Forest Editorial Board, 1997;He et al., 2015). The current distribution areas of these Abies taxa presumably arose from species adaptations, as the species have shifted their ranges over time (Liu, Fang, & Piao, 2002;WWF, 2017).
These four related taxa of Abies in the TP face serious threats due to habitat destruction and climate change (Xiang & Rushforth, 2013a. Rising temperature and extreme precipitation have been highlighted as major factors driving the population dynamics of Abies taxa (He et al., 2015;Wang, Xu, et al., 2017). Further, as explained by Ordonez and Svenning (2016), historic climate data have proved instrumental for describing past events and explaining the current distribution of species. Paleoecology can provide further insight for understanding the long-term influence of environmental changes on Abies taxa and their ecosystems (Tinner & Valsecchi, 2013). Also, accurate responses of plant distributions to historical climate change could serve as good predictors of ongoing climate change (Waltari et al., 2007;Willis & McElwain, 2002). Both Asian Monsoon and uplifting of the TP during the Eocene which created a unique climatic conditions have influenced the distribution and diversification of Abies in this area (Cao, Herzschuh, Ni, Zhao, & Böhmer, 2015;Xiang, Cao, & Zhou, 2007). Yet, studies of the driving factors or dynamics of spatial distribution patterns of Abies taxa and their response to habitat changes at the TP remain surprisingly limited (Xiang et al., 2006). Knowing which factors drove historical changes in Abies distributions might help enhance mitigation strategies for its species and contribute to the conservation of the Tibetan forests in the face of projected climate change scenarios. Ecological niche models (ENMs) can offer perspectives for answering inquiries concerning species' historic and future potential distributions (García-Callejas & Araújo, 2016;Shrestha & Bawa, 2014). Not surprisingly then, ENMs are widely used for species distribution modeling (Dan & Seifert, 2011;Elith & Leathwick, 2009) and for assessing the degree of ecological segregation among different co-occurring taxa (Yi, Cheng, Wieprecht, & Tang, 2014). By applying ENMs to paleo-climate data, we can elucidate the potential past distributions of species (Du, Hou, Wang, Mao, & Hampe, 2017).
Based on pollen fossil records supported by other biological evidence, the climatic oscillations during the Quaternary period significantly influenced biome shifts at the TP (Zhang, Fengquan, & Jianmin, 2000). Evidence from its fossil records revealed the TP served as a vital refugium for numerous plant species that survived the Quaternary glaciations, including those of Abies (Wang & Liu, 1994). The southeast region of the TP provided refugia for the northerly biome during the Quaternary glaciation while it gradually retreated from north to south (Du et al., 2017). Identification of Quaternary refugia is based on various forms of historical biogeographic evidence, particularly those derived from paleo-ecological respectively; the annual temperature range was also a key variable that explained the distribution patterns of the other two Abies taxa. Different adaptation strategies of trees may thus explain the differing patterns of distribution over time at the TP revealed here for endemic Abies taxa.

K E Y W O R D S
Abies forest, ecological niche models, fossils, phylogeography, Quaternary refugia studies. Such evidentiary studies helped to identify glacial refugia in the TP over the most critical periods of the Pleistocene for exemplary taxa including Ginkgo biloba, Pedicularis longiflora, Primula secundiflora, evergreen oaks (Quercus aquifolioides), and Primula sikkimensis (Du et al., 2017;Gong, Chen, Dobes, Fu, & Koch, 2008;Wang, Gong, Hu, & Hao, 2008;Yang, Li, Ding, & Wang, 2008). Knowledge of Quaternary refuge distributions of species long ago (Cao et al., 2015), as well as remaining refuge distributions under current climate change conditions (Tinner & Valsecchi, 2013), can inform the future protection of Abies trees against climate change. Several studies have addressed the biogeographic history of Abies species in the TP (e.g., Peng et al., 2015;Shao, Zhang, Phan, & Xiang, 2017;Xiang et al., 2015), but a thorough understanding of it still eludes us. This may be due, in part, to the lack of using integrative approaches that incorporate both palynological and phylogeographical data. The inherent biases and complexities of the palynological approach-its limitation for providing Abies taxonomic identifications to the species level, dearth of taxonomic precision, inability to account for underrepresented taxa in fossil records, and underestimation of migration distances and shifts in some Abies species distribution-have collectively hindered its ability to infer the location and timing of refugia. This has added to the difficulty in defining the spatial and temporal range of distribution of various species in the past.
This study applied ENMs to integrated paleo-climatic data, fossil pollen records, and previous phylogeographic data of Abies taxa for modeling and describing the current and Quaternary distributions of four endemic Abies taxa. Through the use of this integrative approach, we sought to provide a better understanding of how climate fluctuations have impacted the distribution dynamics of Abies taxa in Hengduan Mountains biodiversity hotspot.
Ecological niche models, Abies paleo-records, and phylogeography each presumes certain stability to some extent in ecological niche dimensions. Nonetheless, the integration of these employed approaches may offer better insight while also improving the accuracy of ENM applications intended to predict the potential Quaternary refugia of the study area. We hypothesized that the different adaptation strategies driven by environmental factors resulted in differing distribution patterns for the four endemic Abies in this area in the late Pleistocene. A better understanding of divergent adaptation strategies of Abies species may thus provide an effective instrument to robustly identify those vulnerable areas for which proactive conservation measures are needed in the high-altitude region of the TP.

| Study area and species
The TP is surrounded by the Himalayan Mountains to the south, the Kunlun and the Qilian Mountains to the north, the Pakistan Karakorum Mountain range stretches along its west side, while the mountains of Hengduan border the east side of the plateau (Zhang, Li, & Zheng, 2002; Figure 1). Subalpine coniferous forest composed of Abies and Picea is the dominant forest type mainly covering the southeastern region of the TP (Li, 1999), particularly in the Hengduan Mountains ranges. Abies species are widespread in this area, dominating the subalpine forests at elevations of 2,500-4,100 m a.s.l. (Li, 1993;Mackinnon et al., 1996).
Abies fargesii var. faxoniana is the representative tree taxon of dark coniferous forests distributed in the northern part of the Hengduan Mountains area, east of the TP (Guan & Zhou, 1983), bordered to the east by the Longmen Mountains, to the west by Jinchuan County, and to the north by the Balang Mountains and Tao River Basin, with an elevation range of 2,500-3,000 m a.s.l. (Xiang et al., 2006). Abies recurvata is an endemic species in the TP; it grows along the basin of Bailongjiang River south of Gansu Province and in the north and northwest of Sichuan Province, across a broad elevation of 2,300-3,600 m a.s.l. Both A. forrestii and A. forrestii var. georgei are the main taxa of the forest south of Hengduan Mountains, at the southeastern side of the TP (China Forest Editorial Board, 1997;Fang et al., 2011). Abies forrestii grows mainly in southwest F I G U R E 1 Current distribution localities of the four endemic Abies across the Tibetan Plateau Sichuan Province and the northern part of Yunnan Province, eastern Tibetan, as well as in Kachin State in Myanmar and Bhutan (Fang et al., 2011). Abies forrestii var. georgei is distributed mainly in southwestern Sichuan and northwestern Yunnan and occurs from the eastern section of the Brahmaputra River in Tibet to downstream of the Yalong River located in Sichuan Province (China Forest Editorial Board, 1997). Abies forrestii var. georgei is found at 3,500-4,500 m a.s.l., while A. fargesii var. faxoniana grows at a lower elevation, from 2,500 to 3,000 m a.s.l. (China Forest Editorial Board, 1997). Recently, A. recurvata and A. fargesii var. faxoniana have been accorded conservation status, being judged vulnerable taxa (Xiang & Rushforth, 2013a, while the status of A. forrestii var. georgei is currently of least concern (Zhang, Katsuki, & Rushforth, 2013b) and that of A. forrestii (syn. A. forrestii var. forrestii) is considered as near-threatened (Zhang, Katsuki, & Rushforth, 2013a). Studies have suggested the taxonomic boundaries are well-established at the species and intraspecific levels of the Abies genus (e.g., Aguirre-Planter et al., 2012;Cinget, Lafontaine, Gérardi, & Bousquet, 2015;Shao et al., 2017;Wang, Abbott, Ingvarsson, & Liu, 2014); however, their intrageneric classification is an active area of research (Eckenwalder, 2009;Farjon, 2010). These four taxa selected for study vary in their morphological traits and have dissimilar environmental requirements, according to information from other studies (e.g., China Forest Editorial Board, 1997;Guan & Zhou, 1983;Shao et al., 2017).

| Taxa occurrence data
Occurrence records of the four Abies taxa were obtained from a field survey conducted in September 2015 and June 2016, and from collected data archived in numerous databases, namely. the GBIF (Global Biodiversity Information Facility, http://www.gbif.org), the CVH (Chinese Virtual Herbarium, http://www.cvh.org.cn/), the NSII (National Specimen Information Infrastructure, http://www.nsii.org. cn/), and other related sources (see all references in the "Additional Reference List" in Appendix S1). Redundant points were removed to ensure a minimum of 1,000 m was maintained as the distance separating data records, to minimize the influence of spatial autocorrelation on the modeling results. In total, 212, 69, 184, and 220 valid effective occurrence records for these four taxa were collected ( Figure 1).

| Environmental data and selection of variables
To represent the environmental conditions in different climatic periods, 19 bioclimatic predictors obtained from the WorldClim database (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005) were used, in addition to two indirect gradients predictors (slope and aspect). used to obtain climatic data for the LGM (21 kyr BP) and the mid-Holocene (6 kyr BP). The LIG data (120-140 kyr BP), with a spatial resolution of 30 arc seconds, were obtained following Otto-Bliesner, Marshall, Overpeck, Miller, and Hu (2006) from the WorldClim 1.4 (http://www.world clim.org; Hijmans et al., 2005). The last two variables corresponded to indirect abiotic gradients: The slope and aspect derived from the digital terrain model (DTM) were extracted from the Shuttle Radar Topography Mission (SRTM) (http://srtm.csi. cgiar.org/) at 30 arc-second resolution. Principal component analysis (PCA) was performed on standardized and centered data of the 19 current bioclimatic predictors and two indirect gradients (Raes et al., 2014). The first three components (PCs) accounted for 98.48% of the variance in the dataset employed in the analysis ( Figure S1).
Multicollinearity can cause the mode to overfit the data used.
Therefore, the Pearson correlations were tested for all pairwise combinations of the 19 climate datasets (Raes et al., 2014). The SDMtoolbox (Brown, 2014;Brown, Bennett, & French, 2017) was used within the framework of ArcGIS 10.0 (ESRI, Redlands, CA, USA) to conduct a correlation analysis among the 19 climate datasets' variables' grids that represented the study area. Only uncorrelated variables with a spatial correlation value below 0.75 were retained (Table S1; Raes et al., 2014). To calibrate the ENMs for each species, nine selected predictive variables were used: mean diurnal range (MDR), temperature annual range (TAR), mean temperatures of driest (TDQ) and wettest (TWQ) quarter, precipitation of wettest (Pmax), driest month (Pmin), precipitation of seasonality (PS), and warmest (PHQ) and coldest (PCQ) quarter, along with another two predictors: slope and aspect.

| Modeling algorithm: Maxent
The maximum entropy (Maxent) modeling technique (Phillips, Anderson, & Schapire, 2006) was applied to the above data, it being one of the most commonly used and accurate ways to predict species' distributions and habitat suitability (Dan & Seifert, 2011;Elith & Leathwick, 2009). The models for each taxon were first calibrated for the current distributions relative to the current climate; the most accurate models were then applied to the historical data (LIG, LGM, mid-Holocene) to get predictions for those past periods. The occurrence records datasets were randomly split (75%-25%), where the 75% portion for each taxon was used to calibrate the algorithm, and the remaining 25% used for evaluating the produced ENMs (Bueno et al., 2017). The version of Maxent used was obtained from http:// biodi versi tyinf ormat ics.amnh.org/open_sourc e/maxen t/; for each taxon, n = 200 replicate runs were undertaken in Maxent 3.4.0.

| Model calibration and evaluation
The ENM models developed to determine the potential distribution of Abies taxa had a logistic output format, which quantifies habitat suitability as a continuous probability value ranging from 0 to 1 (Phillips & Dudik, 2008). The potential distribution curves were estimated for each predictor, and the percentage corresponding to the relative contribution of each one to the habitat suitability models were assessed using Maxent's built-in jackknife test. To evaluate the models and assess their performance, the receiver operating characteristic (ROC) curve was plotted, for which the area under the curve (AUC) was estimated at the possible thresholds (Raes et al., 2014). The ROC curve is generated by plotting the sensitivity versus 1-specificity for the entire set of runs. The AUC provides a widely used metric of model accuracy (Fielding & Bell, 1997). Models attaining an AUC value > 0.7 are deemed of acceptable performance (Swets, 1988).
The response curves for the predictors used in building the models were then generated. Through them, we could gain insight into the quantitative relationships among investigated environmental predictors and the logistic probability of species existence. More specifically, such response curves mirror the favorable ecological niche of the species (Dan & Seifert, 2011).
To minimize the biases to the objective methods that integrated sensitivity with specificity, the maximum training sensitivity plus specificity (Tmax) and the minimum training presence (Tmin) logistic thresholds (Vedelsørensen, Tovaranonte, Bøcher, Balslev, & Barfod, 2013) were each retrieved from Maxent. This was done in addition to using the normal rank grading method (Raes et al., 2014).

| Multivariate analysis of bioenvironmental data
Redundancy analysis (RDA) was conducted for a relative assessment of the different environmental requirements of A. forrestii, A. fargesii var. faxoniana, A. forrestii var. georgei, and A. recurvata. The RDA was performed on the set of environmental variables, based on the their relative contribution percentage to the species habitat suitability models produced by Maxent (mean diurnal range, temperature annual range, mean temperature of wettest/driest quarter, precipitation of seasonality, precipitation of warmest/coldest quarter, and slope), ranked according to their quantitative importance by a forward selection process. Those variables were used to analyze the effect of climate and topography on the pattern of four Abies taxa distributions. The R software environment and its "vegan" package (R Development Team, 2008) were used to carry out this statistical analysis and produce its graphics.

| Paleoclimate, fossil, and phylogeographic records
Climate variables representing the LGM and the mid-Holocene paleoclimates were obtained from CCSM4 and MIROC-ESM general circulation models (GCMs) for inclusion in the analysis. Furthermore, data from the last interglacial period (Otto-Bliesner et al., 2006) were obtained from sporopollen, carbon isotope measurements, calcium carbonate, and ice cores for the area under study (Table   S2). According to our review, the fossil and pollen sequences dated from or around the 6 kyr BP, 21 kyr BP, and 120-140 kyr BP periods in the TP (Table 1). Indications of when the assessed pollen sequence chronology overlapped with the LIG, LGM, or mid-Holocene periods are given in Table 1. Published cpDNA data were used to derive diversity patterns for use in the phylogeographical validation of the model hypotheses (Peng et al., 2012;Zhan, 2011). Since phylogeographical investigations of Abies taxa are limited at the TP, the evaluation of published cpDNA trnL-F, trnS-G, and ndhkC data was possible for only two Abies taxa, namely A. fargesii var. faxoniana (Zhan, 2011) and A. recurvata (Peng et al., 2012) (Table S3; Appendix S1). Additional details on the genetic diversity index and expansion calculations are given as Appendix S1.

| Climate scenarios
The summarized data for 14 deposits of the LGM and the mid-Holocene obtained from the final interpolated models (CCSM4 and MIROC-ESM), and for one point of the LIG climate data for TP, are presented in Tables S2 and S4. Also provided in Table S2 are the discrepancies between GCMs' values through the CCSM4 and MIROC-ESM simulations and the paleoclimates, which were obtained from calcium carbonate, sporopollen, and ice cores for the different climate scenarios applied. Annual mean temperature (Tann) and annual precipitation (Pann) had consistently higher accuracy in the CCSM4 than MIROC-ESM simulation, but with no significant differences for the mid-Holocene in the CCSM4 and MIROC-ESM simulations. In general, the temperature simulations by the CCSM4 and MIROC-ESM models for the LGM were better than those for the mid-Holocene, while precipitation simulations were better than temperature simulations for all periods.   (Tables S4 and 3) strongly affected the aforementioned Abies distributions (Figures 3 and S3), allowing their extension from glacial refugia with an inclination toward shifting to higher latitudes. Yet, this expansion was mostly constrained over the LIG, which indicated that the warm period could have featured unfavorable climatic conditions for these two Abies species (Figures 3 and S3).

| Potential distributions of Abies taxa in the LIG, LGM, and mid-Holocene
The stability surfaces steadily predicted an area extending along the central corridor of the A. forrestii and A. forrestii var. georgei distributions, and stretching from the eastern edge of the Salween River, northwards to the western edge of the Yangtze River (Figure 3), hereon referred to as the Heqing refuge (Figures 1 and 3

| Fossil and phylogeographic records
A classification of the fossils and pollen is given in Table 1.
Concurrence and spatial constituencies were found between the potential distribution and pollen records of Abies species across the TP ( Figures S4 and S5 LGM-CCSM4  Figure S6; Table S3). More details regarding the haplotype, haplotype diversity, and nucleotide diversity for each of the sampled locations are included in Table S3. Each of the aforementioned two taxa has a relatively high total genetic diversity (H T ) and lower within-population diversity (H S ), respectively.
According to the PERMUY analysis, the N ST exceeded the G ST for Abbreviations: MDR, mean diurnal range; PCQ, precipitation of coldest quarter; PHQ, precipitation of warmest quarter; Pmax, precipitation of wettest month; PS, precipitation of seasonality; TAR, temperature annual range; TCQ, mean temperature of coldest quarter; TDQ, mean temperature of driest quarter; TWQ, mean temperature of wettest quarter.

F I G U R E 4
Ordinate plot from redundancy analysis (RDA) on the four Abies change trends (red arrows) and their relationships with environmental variables (blue arrows). MDR, mean diurnal range; Pann, annual precipitation; PHQ/WCQ, precipitation of warmest/ coldest quarter; PS, precipitation seasonality; TAR, temp annual range; TD, mean temp of driest quarter; Tmin, min temp of coldest month; TWQ/TDQ, mean temp of wettest/driest quarter each species, and a significant phylogeographical structure was detected for the A. fargesii var. faxoniana and A. recurvata populations (Table S3). Furthermore, Tajima's D and Fu's Fs tests of neutrality all produced significant negative values (Table S3). These outcomes indicate strong agreement between the observed and expected distributions 9. The elapsed times since the species extension were projected to be 1,542, and 8,086 kyr BP for A. fargesii var. faxoniana and A. recurvata, respectively.

| Environmental predictors of the studied Abies taxa
Estimates of the relative contribution (percentage) of each environmental predictor used to build Maxent models are presented in Table 3.

| Model accuracy and prediction uncertainty
The calibrated potential distribution models for A. forrestii, A. fargesii var. faxoniana, A. forrestii var. georgei, and A. recurvata attained high AUC values ( Figure S2), thus indicating all models had excellent performance (Raes et al., 2014). The modeling technique we used relied on a powerful method (Maxent) that deals strictly with species' occurrence records (Elith & Leathwick, 2009;Phillips et al., 2006). The outcomes from these ENMs might be among the most accurate attainable for the used dataset (occurrence records and environmental data) (Raes et al., 2014 (Li, Wang, & Li, 2013). Therefore, the models presented here could have overestimated the Abies distribution during the LGM and the mid-Holocene (Tinner & Valsecchi, 2013). Yet overfitting should not be construed as low performance, because the models used climatic data alone to predict the potential distributions of the Abies taxa, consequently overlooking other contributing factors (e.g., soil type, land cover, and interspecific competition; Wang, Jia, Wang, Zhua, & McDowell, 2017) that might have been involved in shaping the pattern of distributions; this was clearly suggested by the PCA and RDA results ( Figures S1 and 4). The spatial resolution of the applied models may have posed an additional constraint, mostly in those regions with intricate topography (Adhikari, Barik, & Upadhaya, 2012).

| Abies distribution and refugia based on ENMs, paleorecords, and phylogeography
By combining ecological niche characteristics derived from the environmental characteristics of identified records of Tibetan endemic Abies taxa with their paleoecological investigations, we obtained a more refined depiction of the distribution, discontinuities, and segregation among these tree taxa. This will help provide a greater understanding of the ecology and climatic niche of Abies taxa and could also contribute to improving relevant climate change mitigation strategies. At the LIG, the combined effects of increasing fluctuations of temperature and precipitation throughout the TP (see Table S4) resulted in a contraction and disintegration of the Abies populations (Figures 3 and S3, Table 2). However, these events did not reach the level of severity that would lead to total extinction of the Abies taxa. Notably, topography was a critical factor, and this could have driven the shifting of tree populations in the basins rather than in the plains or mountains, albeit under a stable climate (Willis & McElwain, 2002), when seeking to persist and avoid extirpation. Thus, the Heqing Basin region, in addition to providing refugia for other species (Du et al., 2017), had sufficient topographical heterogeneity to provide many suitable microhabitats for the persistence of A. forrestii and A. forrestii var. georgei. Geographical overlapping among the Abies taxa was widespread throughout the glacial and interglacial times (Du et al., 2017).
During the LGM, the collective effects of reduced annual precipitation and lower summer and winter temperatures across the TP (Table S4) resulted in a shortened growing season and lowered level of atmospheric CO 2 , down to 200 ppm (Braconnot et al., 2007). In addition, and most remarkably, it apparently caused the contraction and fragmentation of A. fargesii var. faxoniana and A. recurvata populations (Figure 3; Table 3). Conversely, the other two Abies taxa expanded their populations in this area during the same period, because smooth fluctuations of temperature and precipitation were more important for them than were changes in other environmental variables ( northward but also westward. The plateau uplifting and climatic fluctuations during the Quaternary overwhelmingly influenced the endemic biome, which shifted from south to north and was associated with glaciation in the late Pleistocene period (Zhang et al., 2000). In parallel, topography was a crucial factor in setting and identifying the range within which populations could expand unimpeded along the continuous mountain ranges (Xiang et al., 2006). In this respect, the Longmen and Hengduan Mountains ranges provided sufficient topographical variability for Abies to spread into other regions during the LGM (Figure 3).
During the mid-Holocene, an extension of the Abies taxa populations is well reverberated by their present outcomes, in that those endemic Abies taxa displayed larger and continuous potential distribution areas relative to those at the LIG and LGM (see Figure 3).
The increased precipitation and warming (  (Du et al., 2017;Gao et al., 2007;Yu et al., 2014). By also considering the indicators of genetic diversity and nucleotide diversity (Zhan, 2011), the populations of A. fargesii var.
faxoniana that attained a high h value and a low Pi value had probably experienced an expansion in their distribution following a long-standing period of suffering from a low effective population size. Accordingly, all of these outcomes converge in supporting the hypothesis that both A. fargesii var. faxoniana and A. recurvata populations at the Longmen areas have expanded ( Figure S6). The Longmen Mountains, however, were not projected to have retained an extensive area offering stable habitat for A. fargesii var. faxoniana and A. recurvate, but rather it served as a temporary refuge for both tree species only at the LGM.

| Ecological variables affecting the distributions of the four related Abies taxa
The four related Abies taxa differed significantly in their climatic niche dimensions (Figure 4). This result is consistent with climate variables driving the distribution strategies of closely related plant species; that is, mean temperature of the coldest quarter significantly impacted the distribution of vegetation type (Wang, Xu, et al., 2017) and the minimum temperature of the coldest month was the main limiting factor for the growth of coniferous temperate forests in southwest China (Li, Peng, & Higa, 2016). In contrast, it was reported that Abies species were cold-tolerant, hydrophilic, sensitive to high temperature, and intolerant to aridity (Wang, Jia, et al., 2017). Such inconsistent results may arise from differences in the adaptation strategies of endemic Abies species to the environmental variables (Table 3; Figure S7). Also, divergent and convergent strategies might lead to an over-or underestimation of the suitable distribution area for combinations of several different co-occurring species, as noted by Willis and McElwain (2002) and Du et al. (2017). Our RDA suggested that A. fargesii var. faxoniana and A. recurvata were better capable of adapting to more magnitude of fluctuation in temperature and precipitation, while the other two taxa are adapted to the higher precipitation of the Eastern Tibetan Plateau (Figure 4). Still, the effect of temperature fluctuations on Abies' distributions was greater than that of either increases or decreases in temperature alone (Table 3).
This conclusion is consistent with the current distribution of these four endemic Abies taxa: In the south, A. forrestii and A. forrestii var. georgei occur in wetter areas, while in the north, A. fargesii var. faxoniana and A. recurvata occupy an area of the TP that undergoes greater temperature magnitude of fluctuation (He et al., 2015;Wang, Jia, et al., 2017). Since most studies tend to focus on the impacts of climatic changes, they simulate increases in average temperature, thereby ignoring the effects of changes in magnitude of fluctuation of temperature (Gazol et al., 2015;Koo et al., 2017;Tinner & Valsecchi, 2013;Xiong et al., 2016). Our results suggest the importance of temperature fluctuations for affecting the growth and distribution of trees, and this variable should be considered when assessing the effects of future changes in climate on plant species. Slope is typically an important ecological factor that influences the distribution of local plant species (Adhikari et al., 2012). Consistent with other regional studies and investigations (Adhikari et al., 2012;Alba-Sánchez & López-Merino, 2010), we found that slope was an important contributor to the distribution of Abies taxa studied here across the TP.

| CON CLUS IONS
The strong agreement of paleorecords with our model predictions demonstrated the merit of using the Maxent modeling approach for projecting the distributions of endemic Abies taxa in the TP. the slope also contributed substantially to the above Abies distribution models. Varied adaptation strategies might well be the reason for differences found in the past and current potential distribution patterns of these four related Abies taxa at the TP.
Incorporating ecological niche characteristics from Abies species combined with their paleoecological investigations provides a useful approach to better understand species migrations under climate changes. It could help guide rational management strategies of forests whose keystone species exist in high-altitude regions.
Future work should focus on whether other endemic tree species (such as Picea, Metasequoia), used the same refuge as, for example, the Heqing refuge from this study. 2019QZKK0303). We thank Professor Liu Guohua for helpful comments and discussions regarding the history of the Tibetan Abies species. We also thank Charlesworth Author Services (http://www. cwaut hors.com) for polishing language in this manuscript. We also appreciate the insightful observations of Ms. Chen Jun, Mr Xiong

ACK N OWLED G M ENTS
Kaiming, Prof. Liang Pinghan, and anonymous reviewers of earlier versions of this paper.

CO N FLI C T O F I NTE R E S T
The contributing authors declare no conflict of interests regarding the publication of this article.