We collected DNA materials for 475 individuals from 42 populations of Chinese subalpine and temperate firs, representing 10 species (Figure S1, Table S2). Four individuals of one population of A. fansipanensis QP Xiang, LK Fu, and Nan Li from Mt Fan Si Pan in northern Vietnam are included in the present study as a typical example of a highly endangered fir with only limited individuals (Xiang, 1997). Within each population, we sampled needles from eleven–twelve individuals growing at least 100 m apart except for six populations (pop. 2, 9, 14, 22, 36, and 43; Figure S1), where we could only sample four–eight individuals. We sampled two to 13 populations each of the following seven species: A. chensiensis Tiegh., A. delavayi Franch., A. fabri (Mast.) Craib, A. forrestii Coltm.-Rog., A. georgei Hand.-Mazz., A. recurvata Masters, and A. squamata Masters. Among them, 181 individuals from 17 populations (Nos. 1–13 and 33–36) of A. chensiensis and A. recurvata have been characterized for mtDNA and cpDNA sequence variation in a previous study (Shao and Xiang, 2015). To address the conservation implications of four other highly endangered fir species, we sampled only one population of each of them due to their limited availability and restrict geographic distributions: A. fanjingshanensis WL Huang, YL Tu, and SZ Fang (IUCN: endangered), A. yuanbaoshanensis YJ Lu and LK Fu (IUCN: endangered), A. ziyuanensis LK Fu and SL Mo (IUCN: endangered) and A. fansipanensis (IUCN: critically endangered; IUCN, 2016). Our sampling also included the following outgroup taxa: A. bracteata endemic to North America as outgroup in the mitotype analyses, and A. alba, A. nordmanniana (from Euro-Mediterranean), and A. mariesii from Japan as outgroup taxa in the chlorotype analyses. Our outgroup selection was based on a previous phylogenetic study (Xiang et al., 2015); only one individual was sampled for each outgroup. We present all sampling efforts in Table S2 and indicate samples that resulted from this study vs. a prior one comprising 181 individuals from 17 populations of A. chensiensis and A. recurvata (Nos. 1–13 and 33–36). Population sampling in the wild was strictly observed according to international rules (Nagoya, 2010).

We extracted total genomic DNA using a modified cetyltrimethyl ammonium bromide (CTAB) protocol (Rogers and Bendich, 1988) or DNA secure Plant Kit (Tiangen). We used the DNA extractions to isolate two mtDNA fragments, nad1 intron 2 and nad5 intron 4, and one cpDNA intergenic spacer, trnS/G, following procedures in Shao and Xiang (2015). We performed sequencing reactions using a DYEnamic ET Terminator Kit (Amersham Pharmacia Biotech), and detected labeled sequencing products on a MegaBACE 1000 (Amersham Biosciences, Buckinghamshire, UK). We deposited all newly generated sequences in GenBank under accession numbers KP635976–KP636041.

We estimated molecular diversity indices, including the average pairwise differences per base pair between sequences (π) (Nei and Li, 1979), and the haplotype diversity (h) (Nei, 1987) for each taxon and population using DNASP version 4.10 (Rozas et al., 2003). Using h and π, we detected patterns of spatial genetic diversity using the GDivPAL function in Spads 1.0 (Dellicour and Mardulyn, 2014). We also calculated haplotype networks among mitotypes and chlorotypes, individually, in NETWORK version 4.2.0.1, with gaps treated as single events (Bandelt et al., 1999).

We used the program PERMUT (Pons and Petit, 1996) to calculate the average gene diversity within populations (H_S), total gene diversity (H_T) and two measures of population differentiation (G_ST and N_ST). When N_ST estimates were significantly higher than the G_ST-values, we inferred a phylogeographic structure. We performed an analysis of molecular variance (AMOVA) using the program ARLEQUIN 3.5 (Excoffier and Lischer, 2010) to determine the amount of genetic variation within and among the taxa and populations. In order to test the correlation between genetic and geographic distances, the mantel test was performed (Mantel, 1967; Nei and Li, 1979). We identified population groups that were geographically homogenous and maximally differentiated from each other using spatial analysis of molecular variance (SAMOVA; Dupanloup et al., 2002). We detected distinct populations (K) via a simulated annealing approach to maximize among population variance (F_CT), and we chose the K between 2 and 10 where F_CT began to plateau.

We performed mismatch distribution analyses in the program ARLEQUIN 3.5 to detect historical population expansion events. When we detected an expansion event, we estimated its age using the formula s = 2 ut (Rogers and Harpending, 1992) where t is the age of the expansion, u = μkg (μ is the substitution rate per nucleotide site per year, k is the average sequence length, and g is the generation time in years). In this study, we estimated expansions using a range of neutral substitution rates for cpDNA (μ = 2.61 × 10⁻¹⁰–4.02 × 10⁻¹⁰ ss⁻¹ year⁻¹) (Gernandt et al., 2008), an average haplotype sequence length (1,111 bp) in Pinaceae (Brown et al., 2004). In Pinaceae, the minimum estimate of generation time was 25 years (Brown et al., 2004). However, it should be mentioned that sexual maturity and effective generation time (when seedlings are actually recruited) of fir species in wild populations was also likely to be around 50 years (Cheng and Fu, 1978; Fan, 2006). Here, we calculated the historical expansion time under the two different minimum estimate of generation time (option 1: g = 25 years; option 2: g = 50 years).

We performed ecological niche modeling (ENM) of Chinese subalpine and temperate firs in MAXENT (version 3.2.1; Phillips et al., 2006). The ENMS were not performed on species level for the extremely limited reliable occurrence records except for A. chensinesis (Fan, 2006; Shao and Xiang, 2015). For input into MAXENT, we obtained species occurrence data from the Chinese Virtual Herbarium (CVH, www.cvh.org.cn) and personal observations, including 112 records. These records involved 36 A. chensiensis records, 15 A. delavayi records, 11 A. fabri records, 8 A. forrestii records, 9 A. georgei records, 20 A. recurvata records, seven A. squamata records, 1 A. fanjingshanensis record, 1 A. yuanbaoshanensis record, 3 A. ziyuanensis records, and 1 A. fansipanensis record.

We used 19 BIOCLIM variables as environmental data for ENM representing five different time periods. Past and current data for these variables were available from the WorldClim database (www.worldclim.org). The data from WorldClim included baseline climate (BioClim layers for the period 1,950–2,000) at a spatial resolution of 30 arc s), data for the last glacial maximum (LGM; ~21,000 years BP) with spatial resolution of 2.5 arc min resolution simulated by CCSM model and MIROC model, and the last interglacial (LIG; ~120,000–140,000 years BP) period with spatial resolution in 30 arc s resolution. We obtained predictions for future climatic conditions in 2,050 and 2,080, from the International Center for Tropical Agriculture (www.ccafs-climate.org) based on the HadCM3. The HadCM3 model shows two scenarios (A2, and B2), which were represented in 2.5 arc min resolution Many climatic variables were highly correlated (|Pearson's R| ≥ 0.8). Therefore, our final ENMs were constructed using seven variables (mean diurnal range of temperature, temperature seasonality, max temperature of warmest month, mean temperature of driest quarter, mean temperature of coldest quarter, precipitation of driest, and warmest quarter) with limited covariance.

To construct ENMs, we conducted ten replicate runs with these parameters: default convergence threshold, maximum iterations (1,500), and 30% of the sites for model training (Waltari et al., 2007; Jezkova et al., 2009). We evaluated model performance using the area under the receiver operating characteristic curve (AUC) calculated by MAXENT.

We initially surveyed more than 200 Quaternary fossil pollen records from both publications and database sources involving distinct species and time span, following Cao et al. (2013). During the last glaciation, Abies spp. occurred in 10 pollen records with reliable chronologies and sufficiently high sampling resolutions, covering the subtropical and temperate China (Liu and Ye, 1977; Zhou et al., 1978; Zhang et al., 1997; Sun and Li, 1999; Wu et al., 2002; Xu et al., 2002; Feng et al., 2007; Tang et al., 2007; Cook et al., 2011; see details in Table S1).

The length of the mtDNA fragments, nad1 intron 2 and nad5 intron 4, was conservative; 675 bp and 220 bp, respectively. The fragments contained seven point mutations and four indels resulting in 11 distinct mitotypes (M1–11; Table S3). Our results show that M1 has the widest distribution and is likely the ancestral haplotype according to the NETWORK analysis. M1 was the dominant haplotype in the central Hengduan Mountains region (Figure 1) and also occurred in Vietnam and southeast China which harbored four rare and critically endangered species, A. fanjingshanensis, A. fansipanensis, A. yuanbaoshanensis, and A. ziyuanensis (Pop. 40–43; Figure 1).

The estimates of molecular diversity π and h, are summarized in Table S2. Diversity of mitochondria was high in several distinct geographic regions including the Hengduan Mountains, the Qinling-Daba Mountains, and Southeast China according to the results from Spads (green dots in Figures 2A,B). Within each species, the distribution of mitochondria diversity also didn't show any latitudinal trend (Table S2; Figures S4–S6). Our AMOVA analyses revealed that 12.36% of the total mtDNA sequence variation occurred among taxa, whereas 78.37% of the variation occurred among populations within taxa (Table 1). The mtDNA sequence variation showed H_S much lower than H_T. Furthermore, mitotypes exhibited strongly isolated geographic distributions and a significant phylogeographic structure (G_ST < N_ST; P < 0.05; Table 2). Genetic diversity was significantly correlated with geographic distances (r = 0.201, p = 0.037) according to our Mantel test (Table 2). However, no significant geographic division was revealed by the SAMOVA analyses.

In the chloroplast DNA data, we detected six indels and 17 substitutions resulting in 25 chlorotypes (C1–25; Table S4). Five chlorotypes had frequencies above 5% (C1, C7, C11, C14, and C19). Fifteen chlorotypes (C2, C4–6, C8–10, C12, C18, and C20–25) were fixed in single species, while the remaining ten were shared among species (Figure 1B). The chlorotype C1 possessed the widest geographic distribution and may be ancestral haplotype according to the NETWORK analysis (Figure 1).

The cpDNA exhibited higher genetic diversity than the mtDNA according to the Spads analysis (Figure 2; cpDNA in green). The cpDNA possessed multiple centers of diversity according to the indices h and π (Figures 2C,D; yellow and red dots). Within each species, the distribution of cpDNA diversity also didn't show any latitudinal trend (Table S2; Figures S4–S6). Of the variation among chlorotypes, 38.14% was among taxa and and 39.44% was within populations (Table 1). Chlorotypes exhibited significant phylogeographic structure (G_ST < N_ST; P < 0.05), and H_S was significantly lower than H_T (Table 2). However, we did not detect a significant division among populations according to SAMOVA analysis, but the Mantel test revealed a significant correlation between genetic and geographic distances (r = 0.227, P = 0.010; Table 2).

We performed mismatch distribution analyses on Chinese subalpine and temperate firs taken together and each species individually using cpDNA sequence data. We analyzed seven species for mismatched distributions (A. chensiensis, A. delavayi, A. fabri, A. forrestii, A. georgei, A. recurvate, and A. squamata) and excluded the four rare and endangered species (A. fanjingshanensis, A. fansipanensis, A. yuanbaoshanensis, and A. ziyuanensis) because of the small sample size (No. pop = 1).

The mismatch distribution results rejected the demographic expansion hypothesis for the Chinese subalpine and temperate firs as a whole (P_SSD = 0.000; Figure S2). However, the expansion hypothesis was supported for all of the seven species analyzed separately except for A. delavayi (P_SSD = 0.078; Figure S2). We determined that the demographic expansion of the six species probably occurred in 64.9–18.2 kya (option 1: g = 25 years) or 32.5–9.1 kya (option 2: g = 50 years; Figure 3, Table 3).

The value of the AUC test (mean ± SD) for the ecological niche modeling, averaged across all 10 runs, was very high (0.96 ± 0.005). Models for the present (Figure 4A) and last interglacial (Figure 4B) showed contracted ranges compared to the LGM (Figure 4C), indicating losing of suitable habitat in south, central, and east China. The future predictions of suitable habitat were narrower than for the present day, and ranges were further reduced between 2,050 (Figure 4D) and 2,080 (Figure 4E). Except the Models for the LGM, the suitable habitats of the four highly endangered firs in southeast China were always absent (Figure 4). As for other seven species, the suitable habitats were slightly changed in all time periods bun not in the LIG. During the LIG when globe warming, suitable habitat of fir species distributed in northern part of the Hengduan Mountains region was contracted heavily, such as A. recurvata, A. fabri, A. forrestii, and A. squamata (Figure 4, Figure S1).

Palynological data indicated that Abies were much more common in lowland forests of central (no. 2,3,8) and eastern (no. 1,4,5,7,9) China during the last glaciation (Figure S3). This pattern depicted wider distribution during the last glaciation than today, Except for three pollen records from northeastern China (no. 1, 4, and 8), the pollen-based distribution pattern was highly consistent with the ENMs for the LGM by CCSM model (Figure 4C).

Previous studies concerning arctic-alpine species from Europe and North America such as Sedum lanceolatum (Dechaine and Martin, 2005), Asplenium fontanum (Bystriakova et al., 2014), and Abies cilicica (Awad et al., 2014), indicate glacial expansion of elevational shift scenario. This “elevation shift” scenario had not been well-documented in arctic-alpine species of East Asia, and places constraints on its application to biodiversity conservation. By combining molecular phylogeography and ENM, we explored how the Chinese subtropical and temperate firs (typical alpine species) in eastern Asia responded to Quaternary climate change.

Two distinct expansion models have been proposed by Hewitt (1996) to distinguish elevational shift from latitudinal shift. The pioneer expansion model of latitudinal shift is characterized by decrease of genetic diversity along latitudes during interglacial expansion (Hewitt, 2004; Yan et al., 2012; Andreia et al., 2016; Nicolas et al., 2017). The phalanx expansion model of elevational shift supposes glacial expansions without significant reduction of genetic diversity along latitudes, resulting in multiple centers of diversity but without decreased genetic diversity along latitude (Galbreath et al., 2009; Awad et al., 2014; Bystriakova et al., 2014; Inoue and Berg, 2016). Based on current phylogeographical patterns, the phalanx expansion model via elevational shift was potentially active in Chinese subtropical and temperate firs. Firstly, multiple centers of fir diversity were discovered in the Hengduan Mountains region, Qinling-Daba Mountains region and in Southeast China, indicating no latitudinal decreasing trends (Figure 2). A similar pattern was also found within each species (Table S2; Figures S4–S6). Secondly, relatively low genetic diversity of mtDNA and cpDNA were found both within populations and species across Chinese subalpine and temperate firs (Table S2). Such low genetic diversity or strong homogeneity could be the result of rapid recolonization following severe bottleneck effects, accompanied by climate-driven elevation shifts along a local elevation gradient (Hewitt, 2000; Jaramillo-Correa et al., 2004; Jiang et al., 2011). This result is consistent with similar processes reported in the Mediterranean region and North America (Ziegenhagen et al., 2005; Jaramillo-Correa et al., 2008; Wang et al., 2011). In addition, the population diversity structure of Abies in subtropical and temperate China is similar to that of congeners in North America and Europe especially by having differentiation among populations (mtDNA: 90.73%; cpDNA: 60.56%; Table 1) significantly greater than the averages of conifers (mtDNA: 76.4%; cpDNA: 16.5%; reviewed in Petit et al., 2005) (Liepelt et al., 2002; Fady-welterlen, 2005; Jaramillo-Correa et al., 2008; Fady and Conord, 2010; Wang et al., 2011). Greater differentiation of conifers among populations may result from recent habitat fragmentations or severe bottleneck effects due to global warming since the Holocene epoch (Hewitt, 2004; Jaramillo-Correa et al., 2004; Jiang et al., 2011). Thus, the phylogeographical pattern of Chinese subtropical and temperate firs could be considered potential evidence of elevational shift, in general, and of the phalanx model in particular.

More directly evidences were further revealed by the expansion time estimation and ecological niche modeling. We estimated that the ranges of most firs (i.e., taken individually, except A. delavayi) expanded during approximately 64.9–18.2 kya (option 1: g = 25 years) or 32.5–9.1 kya (option 2: g = 50 years; Table 3; Figure S2). These two timing of the expansion were highly coincident with the last glaciation in China, which was 73–10.4 kya with a maximum during 21–18 kya (Yi et al., 2004, 2005, 2007; Zhao et al., 2011). Specifically, range expansion of Abies in subtropical and temperate China might have commenced shortly after the beginning of the last glacial cycle and continued throughout the last glacial maximum (option 1) or since the middle of the last glaciation till the end (option 2; Figure 3). Both the two options indicated the glacial expansions via elevational shift.

Curiously, we did not detect evidence for recent, historical range expansion in all subtropical and temperate China and Vietnam firs considered together or in A. delavayi. For A. delavayi, distribution fragmentation could be easily inferred by its typical high-altitude distribution, extending from the southern part of Hengduan Mountains region to north Vietnam. For the whole Chinese subtropical and temperate firs, multiple asynchronous expansions, multiple independent refugia in situ and complex topography in this region might be responsible for the undetectable in mismatch distribution analyses (Cun and Wang, 2010; Yan et al., 2012; Qin et al., 2013). Such analyses seemed to underestimate the possibility of historical expansions in some situations. Nevertheless, our historical distributional projections showed that firs of subtropical and temperate China and Vietnam as a whole underwent an expansion from the LIG to the LGM; that is, during a period of global cooling. Distributional projections for the present and future show ongoing range contractions (Figure 4) coincident with continued global warming. Although such analysis was not conducted at the species level due to insufficient occurrence records except for A. chensiensis, we believed that this tendency should be well-represented in respective taxa because of their conservative ecological niche requirement (Cheng and Fu, 1978; Farjon, 2001; Fan, 2006). The accuracy of ENMs was further validated by the high similarity between the LGM-CCSM prediction (Figure 4C) and pollen-based distribution during the last glaciation (Figure S3). Two prior studies modeled future geographic distributions of six fir species of China, Japan, and Pakistan, and determined that those species would all lose their suitable habitat in the lowlands region in response to ongoing global warming (Tanaka et al., 2012; Ali et al., 2014).

Many plant species have strong preferences for a narrow range of climatic conditions. Therefore, past climatic conditions can be inferred from the palynological records of plants (Tinner and Lotter, 2006). Palynological data indicate that Abies were much more common in lowland forests of central and eastern China during the last glaciation (Figure S3) when the mean annual temperature was 5–10°C lower than today (Xu et al., 1980; Shi, 1996; Cao et al., 2013). Moreover, the abundance of Abies pollen in central China decreased with the onset of the interglacial period, which was warmer and less humid than the glacial one. The current mountain top distribution of several highly endangered fir species (e.g., A. beshanzuensis, A. fanjingshanensis, A. yuanbaoshanensis, and A. ziyuanensis) is living evidence of interglacial fragmentation (Cheng and Fu, 1978; Xiang, 2001; Jiang et al., 2015). Even though deforestation in China had been significant in past decades, these fir species were not a main target due both to their distribution along mountain tops or along the tree line away from human settlements, and to the low-quality wood properties (Fan, 2006). In summary, our results including genetic diversity, mismatch distribution analyses, ecological niche models, combining evidence from other lines involving fossil records and the current distribution patterns, led to a conclusion that the Chinese subtropical and temperate firs experienced elevational shift during the Quaternary climate change.

Our results show that Chinese subtropical and temperate firs have undergone major habitat fragmentation and genetic isolation and differentiation since the Holocene. The main goal of biological conservation is to maintain as many of the important genetic building blocks of the species as possible so that the evolutionary progress is not constrained (Fraser and Bernatchez, 2001; D'Amen et al., 2013). Thus, defining conservation units with evolutionary differences is necessary for these economically and ecologically valuable threatened firs (e.g., de Querioz, 1998; Fraser and Bernatchez, 2001).

In subtropical and temperate China and Vietnam, four species of Abies are listed as critically endangered due to their restricted distribution on isolated mountain tops and small population sizes: A. beshanzuensis, A. fanjingshanensis (Pop. 40), A. yuanbaoshanensis (Pop. 41) and A. ziyuanensis (Pop. 42) in warm, temperate China and A. fansipanensis (Pop. 43) in Vietnam (Fu, 1992; Xiang, 2001). Abies beshanzuensis was not included in this study because its population comprises only three individuals and is too ecologically fragile for sampling. The other four species all possessed the M1 ancestral mitotype but had high cpDNA diversity (Figure 1). This observation implied a possible shared ancestral gene pool which was fragmented during interglacials (Jaramillo-Correa et al., 2006; Loehr et al., 2006; Dubreuil et al., 2010). Specifically, the total number of cpDNA haplotypes was three, two, four, and one in A. fanjingshanensis (C1, C3, C7), A. yuanbaoshanensis (C1, C3), A. ziyuanensis (C1, C5, C7, C11), and A. fansipanensis (C19), respectively (Figure 1). These distinct chloroplast haplotypes in each species highlight the distinct evolutionary histories of the species (Liepelt et al., 2009; Jay et al., 2012), thus each may represent separate conservation units.

For the species sensitive to environmental change, conservation planning must take into consideration the climatic envelope of the species and its present geographic distribution (Root et al., 2003; Thomas et al., 2004; Forest et al., 2007). The climatic envelope of Abies has been shrinking in geographic space since the last ice age and will probably continue to do so under present global warming conditions (Figure 4). The reasons are various, such as pest and disease damage, human disturbance, and habitat destruction. Unfortunately, Abies may have limited opportunity to escape global warming through elevational shift, because most species already occur at or nearly at the maximum available altitude: For example, A. beshanzuenesis occurs on Baishanzu Mountain (1,857 m) at 1,750 m, A. fanjingshanensis is restricted to 2,050–2,390 m on Fanjing Mountain (2,494 m), A. yuanbaoshanensis is distributed on Yuanbao Mountain (2,081 m) at 1,900–2,000 m, and A. ziyuanensis is located on Shunhuang Mountain (1,882 m) at 1,700 m (Fan, 2006; Li et al., 2012). While some species of Abies can move upwards, they may not be able to do so at a fast enough rate (Huck et al., 2012; Moritz and Agudo, 2013). Based on the meta-analyses, Chen et al. (2011) proposed that the distributions of many terrestrial organisms currently shifted to higher elevations at a median rate of 11.0 m per decade or to higher latitudes at a median rate of 16.9 km per decade to accomodate global warming. Such a high rate is difficult to achieve for conifers, which have a two (or longer) year reproductive cycle (Liu, 1971; Farjon and Rushforth, 1989; Xiang, 2001).

In situ conservation is the current govermentnal strategy for these Abies species, but it is clearly not a permanent solution. In fact, there is not any suitable habitat for these highly endangered firs in actual or 2,050/2,080 models (Figure 4). Moreover, a massive suitable habitat contraction occurred in the northern part of the Hengduan Mountains region during the LIG (Figure 4B). Facing the forecasted global warming, the suitable habitats of A. recurvata, A. fabri, A. forrestii, and A. squamata were believed to be heavily contracted (Figure 4B, Figure S1). Thus, we propose that these species urgently require ex situ management in botanical gardens in areas of suitable climate with dedicated breeding programs. According to our ENMs (Figure 4), the suitable habitats for breeding firs were stable in the southern Hengduan Mountains region during each time periods. At present, more than 30 National Nature Reserves are located in the Hengduan Mountains, and also some relevant state-grade research institutes, such as the Kunming Institute of Botany, Chinese Academy of Sciences (Chen et al., 2009). These specific gardens and their existing programs would strengthen the development of ex situ conservation.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.