-
PDF
- Split View
-
Views
-
Cite
Cite
Guo-Hong Wang, Hai-Wei Zhao, Meng An, He Li, Wei-Kang Zhang, Detecting the driving forces underlying the divergence of spruce forests in China: evidence from phytocoenology, morphology and phylogenetics, Journal of Plant Ecology, Volume 13, Issue 1, February 2020, Pages 59–69, https://doi.org/10.1093/jpe/rtz047
- Share Icon Share
Abstract
We aimed to elucidate the driving forces underlying the geographical distribution of spruce forests, as well as the morphological and phylogenetic divergence among spruce species in China.
One hundred and seventy two sites across the entire range of spruce forests in China (23°–53° N, 75°–134° E, 250–4300 m a.s.l.) were sampled for species composition, geographical coordinates, and topographic and climatic variables. Sixteen spruce taxa, which are naturally distributed in China, were respectively grouped into morphologically defined sections and phylogenetically distinct clades. Multivariate approaches, including two-way indicator species analysis, principal components analysis, detrended correspondence analysis, canonical correspondence analysis (CCA), and partial CCA, were used for data analysis.
The 172 samples grouped into 13 spruce forests, the geographical distributions of which were closely related to climate and geographical location. The variation in species composition explained by the geographical coordinates (32.01%) was significantly higher than that explained by the climatic (27.76%) and topographic variables (23.32%). Of the three morphologically defined sections, sect. Omorica occurred mainly in wetter habitats with a mean annual precipitation of ca. 229 mm and 426 mm higher than the habitats of sect. Casicta and sect. Picea (P < 0.01), respectively. Of the two phylogenetically distinct clades, Clade-II (an older clade) occurred in habitats with warm winters and cool summers whose mean temperature in the coldest month was ca. 8–10°C higher, yet accumulated temperature during the growing season (≥ 5°C) was ca. 297–438°C lower, than the habitats of Clade-III (a younger clade) (P < 0.01). Our data support the hypothesis that geographical location may be a greater determinant of variation in species composition. In addition, moisture conditions tend to be the key determinants that account for the divergence among the morphologically defined sections, while the phylogenetic divergence among spruce species is mainly affected by temperature conditions. While the clades or sections of the spruce species in question carry strong climatic signals, their divergences are subject to different selective pressures.
摘要
本研究在中国云杉林的全部自然分布范围内(23°–53°N,75°–134°E,海拔250–4300 m),采样调查了172个样方,收集其群落的物种组成、地理坐标、地形信息和样地的气候数据。这16个自然分布的云杉类群,被分为不同形态学组和系统发育进化枝。采用多元分析方法,包括双向指标物种分析、主成分分析、除趋势对应分析、典型对应分析(CCA)和部分典型对应分析,对中国云杉林地理分布、云杉属形态分化和系统发育分异与环境因子间的关系进行研究,以揭示中国云杉林的地理分布规律,阐明各群系建群种的形态分化和系统发育分异的驱动机制。结果表明,172个样方可归为13个云杉林群系,各群系的地理分布与经纬度和气候因子密切相关;经纬度对中国云杉林物种组成变化的解释量为32.01%,显著高于气候因子(27.76%)和地形因子(23.32%)。在三个形态学组里面,鱼鳞云杉组(sect. Omorica)的物种主要分布在较湿润的区域,年均降水量比丽江云杉组(sect. Casicta)多229 mm,比云杉组(sect. Picea)多426 mm。在中国云杉属的两个系统发育进化枝中,较古老的分枝(Clade-II)主要分布于温度年较差较小的区域,其最冷月均温比较年轻的分枝(Clade-III)高8–10°C,但前者生长季≥5°C有效积温比后者低约297–438°C。本研究结果支持地理位置可能是物种组成分异的主要决定因素这一假设。此外,湿度是云杉属各形态组间分异的决定因子,而云杉属系统发育分异则主要受温度的影响。可见,尽管云杉属的形态分化和系统发育分异均受气候的影响,但两者的选择压力并不相同。
INTRODUCTION
It has long been established that climate is a key environmental variable that influences plant distribution (Transeau 1909; Woodward and Williams 1987; Woodward et al. 2004) and community composition (Gauthier et al. 2000; Økland et al. 2003; Zhao et al. 2018). However, there are noteworthy exceptions to this. For example, Ammopiptanthus mongolicus, which is an evergreen shrub that grows naturally in deserts, is likely a Tethyan relic as opposed to a product of climate selection (Sun and Li 2003). Similarly, the ranges of some European flora and faunal taxa are more strongly determined by the location of three glacial refugia than by current climate (Araújo et al. 2005). In fact, historical factors, such as the vicissitudes of past climate that are associated with the advance and retreat of ice sheets during the late Tertiary and Quaternary, have played an important role in determining plant distribution (Hewitt 2000; Walker 1986), with species becoming extinct, surviving in refugia or dispersing to new locations (Hampe and Petit 2010; Hewitt 2000). Species that dispersed to new locations might have undergone severe selection by climate, while those that survived in refugia should be highly compatible with their current geographical locations. Determination of the respective roles of current climate and geographical location in the characterization of current vegetation distribution would provide an improved foundation for predicting plant responses to future climatic changes (Araújo et al. 2005; Liang et al. 2018; Svenning and Skov 2004).
Spruce (Picea) forests, that is, coniferous forests dominated by a single spruce species or numerous species that include spruce, are important components of boreal (Taiga) vegetation and subalpine coniferous forests. Spruce forests have a wide geographical range that covers the Northern Hemisphere and extends from the Eurasian continent to Northern America (Toby and Chytrý 2002), of which eastern Asia is a diversity center of spruce species (Farjón 2001). Although the origin of spruce species is debated, an increasing body of evidence has shown that spruce likely originated in the early Tertiary or late Cretaceous, while the ancestor of present spruce dates to the Oligocene (LePage 2001; Lockwood et al. 2013; Sigurgeirsson and Szmidt 1993; Ran et al. 2006). The current distribution of spruce forests in Asia was mainly affected by the uplift of the Tibetan Plateau in the late Tertiary and the last glacial maximum in the Quaternary (Tsukada 1966; Wu et al. 2007; Xu et al. 1980).
Almost 16 spruce taxa are naturally distributed in China (Fu et al. 1999), most of which are dominant species in cold temperate coniferous forests. Spruce forests in China grow mainly in subalpine areas but have a wide geographical range, with different species having different ranges in different climatic areas (Zhou and Yang 1980). In addition, the 16 spruce taxa can be grouped into three sections according to their morphological divergence (Supplementary Appendix S1; Zheng and Fu 1978). Alternatively, they can be classified into phylogenetically distinct clades (Supplementary Appendix S3; Lockwood et al. 2013; Ran et al. 2006). Spruce species within a section, however, are not necessarily more similar in terms of their phylogenetic relationship than those between sections (Bailey et al. 2015). Given this, we hypothesized that morphological and phylogenetic divergence among the spruce species in question might be subject to different selective pressures. In addition, while current climate and geographical locations are mutually correlated, the latter may be a better determinant of vegetation distribution, as some relics of the glacial refugia that are restricted to specific locations tend to have poor relationships with current climate (Hampe and Petit 2010; Hewitt 2000; Svenning and Skov 2004).
In this study, we used data from field surveys conducted in spruce forests across China to test this hypothesis by answering three questions: Do species compositions differ between spruce forests in different climate regions? What are the relative contributions of current climate and geographical locations in determining the distribution of spruce forests? Are morphological and phylogenetic divergences among the spruce species in question subject to different selective pressures?
MATERIALS AND METHODS
Study area and field survey
Sampling sites were selected within the entire range of spruce forests on the Chinese mainland and Taiwan (23°–53° N, 75°–134° E, 250–4300 m a.s.l.; Fig. 1; Zhang 2007). A total of 172 sites were available, of which 163 were sampled between 1983 and 2012, and information on the additional nine sites was obtained from published data (Tseng 1991; Yang et al. 2002). All sampled sites had visually homogeneous vegetation with minimal signs of recent disturbance. At each site, geographical coordinates, altitude, terrain, slope gradient and slope aspect (see Supplementary Appendix S2) were accurately recorded. A few sites whose geographical coordinates were not available in the original records were revisited and resampled using Google Earth System with the aid of the originally recorded data. The species composition was surveyed at each site. At most sites, we surveyed one 20 × 20 m2 quadrat for the tree layer, and four to five quadrats for the shrub (5 × 5 m2 each) and herb layers (1 × 1 m2 each). The quadrat size for the tree layer (i.e. 20 × 20 m2) should be sufficient for most of the spruce forests in China because the number of tree species remained unchanged when the quadrat size was increased from 20 × 20 m2 to 1 ha. However, a larger sample size for the tree layer (i.e. 30 × 30 m2) was used for sampling spruce forests in northeastern China and Taiwan because of the higher species richness in these regions. We created a checklist of vascular species present at each site. Spruce species were the dominant components in the sampled communities. Consequently, a matrix of sample by species (presence and absence) was available for further analysis.

Sampling sites of spruce forests in China. Sites that were sampled across the entire range of spruce forests in China are shown as closed black triangles. Elevation gradients are indicated using color fields.
Climatic variables
The original observations of monthly mean temperature, monthly precipitation and percentage of possible sunshine hours were derived from 1814 meteorological stations across China. The data from 740 stations were observed from 1971 to 2000, and the remaining data were observed from 1981 to 1990 (China Meteorological Administration, unpublished data). These climate data were interpolated at a 10′ resolution using a smoothing spline method (Hutchinson 1989). The mean annual temperature (MAT) and mean annual precipitation (MAP) were extracted for each site. Seven bioclimatic variables were calculated from these interpolated data using the global equilibrium biogeography-biogeochemistry BIOME3 model (Haxeltine and Prentice 1996): actual evapotranspiration (AET), potential evapotranspiration (PET), mean temperature of the coldest month (TCM), mean temperature of the warmest month (TWM), ≥ 0°C and ≥ 5°C accumulated temperature during the growing season (GDD0 and GDD5) and soil moisture (SM). The input soil texture data for BIOME3 were digitized from published data (Xiong and Li 1987). We also calculated the climatic moisture index (MI = MAP/PET) and plant water availability (α = AET/PET) for each site.
Data analysis
In this study, we focused on the inter-alliance variation in species composition. The precision of vegetation classification was therefore set on the level of alliance. For simplicity, the short name of an alliance is used throughout, for example, ‘Alli. P. morrisonicola’ represents the alliance of Picea morrisonicola (Taiwan spruce forest).
We performed a two-way indicator species analysis (TWINSPAN) to classify individual samples into subgroups using WinTWINS software (Hill and Šmilauer 2005). The divisions of TWINSPAN were terminated when the samples within a specific subgroup shared the same spruce species. Such a subgroup therefore corresponds to an alliance. The final subgroups were largely determined by the results of TWINSPAN but were slightly modified by expert judgment. Specifically, Alli. P. jezoensis var. microsperma was merged into Alli. P. jezoensis var. komarovii, and similarly, Alli. P. brachytyla var. complanata was merged into Alli. P. brachytyla, while Alli. P. likiangensis var. linzhiensis was merged into Alli. P. likiangensis. These merged alliances either overlapped geographically or were closely located and had highly similar dominant species in terms of morphological characters.
Factor analysis was conducted to extract the major climatic variables. In doing so, a principal components analysis (PCA) of the nine climatic variables (MAT, MAP, TCM, TWM, GDD0, GDD5, SM, MI and α) was performed using the SPSS statistical package (SPSS, Chicago, IL, USA).
Two ordination methods were used to detect the relationship between environmental variables and spruce forests: detrended correspondence analysis (DCA), which is an indirect gradient analysis that reveals the gradients in species composition (Hill and Jr 1980), and canonical correspondence analysis (CCA), which is a direct gradient analysis that reveals how the variation in species composition is influenced by environmental variables (ter Braak 1986). Significance tests in CCA were conducted using the Monte Carlo permutation test.
To reduce the arch effect of the CCA, three climatic variables, that is, GDD5, MAP and TCM, which represented the three dominant climatic gradients revealed by the PCA (Table 1), were retained in the constrained analysis, while those that were strongly correlated with the three variables were excluded. Given that the variables MAP and SM were highly correlated, we selected the more widely used MAP to represent the second climatic gradient. In addition, the geographical location (longitude, latitude and altitude) of the constrained variables was set as a group to ensure the clarity of the boundaries between different spruce alliances. It should be noted that the overall ‘sample-environment’ pattern revealed by the biplot of the CCA tended not to vary with geographical location (data not shown).
PCA of climate variables used in the factor analysis, showing the correlations of climate variables with the three components (1, 2 and 3)
Component matrix . | Component 1 . | Component 2 . | Component 3 . |
---|---|---|---|
MAP | 0.039 | 0.884 | 0.428 |
SM | −0.181 | 0.901 | 0.085 |
MAT | 0.469 | 0.126 | 0.865 |
TCM | −0.157 | 0.189 | 0.954 |
TWM | 0.962 | −0.122 | −0.189 |
GDD5 | 0.985 | −0.097 | 0.060 |
GDD0 | 0.954 | −0.038 | 0.233 |
MI | 0.296 | 0.630 | −0.534 |
Α | 0.781 | 0.424 | −0.057 |
Eigenvalues | 3.112 | 2.249 | 1.834 |
% of Variance | 38.900 | 28.107 | 22.919 |
Cumulative % | 38.900 | 67.008 | 89.927 |
Component matrix . | Component 1 . | Component 2 . | Component 3 . |
---|---|---|---|
MAP | 0.039 | 0.884 | 0.428 |
SM | −0.181 | 0.901 | 0.085 |
MAT | 0.469 | 0.126 | 0.865 |
TCM | −0.157 | 0.189 | 0.954 |
TWM | 0.962 | −0.122 | −0.189 |
GDD5 | 0.985 | −0.097 | 0.060 |
GDD0 | 0.954 | −0.038 | 0.233 |
MI | 0.296 | 0.630 | −0.534 |
Α | 0.781 | 0.424 | −0.057 |
Eigenvalues | 3.112 | 2.249 | 1.834 |
% of Variance | 38.900 | 28.107 | 22.919 |
Cumulative % | 38.900 | 67.008 | 89.927 |
PCA of climate variables used in the factor analysis, showing the correlations of climate variables with the three components (1, 2 and 3)
Component matrix . | Component 1 . | Component 2 . | Component 3 . |
---|---|---|---|
MAP | 0.039 | 0.884 | 0.428 |
SM | −0.181 | 0.901 | 0.085 |
MAT | 0.469 | 0.126 | 0.865 |
TCM | −0.157 | 0.189 | 0.954 |
TWM | 0.962 | −0.122 | −0.189 |
GDD5 | 0.985 | −0.097 | 0.060 |
GDD0 | 0.954 | −0.038 | 0.233 |
MI | 0.296 | 0.630 | −0.534 |
Α | 0.781 | 0.424 | −0.057 |
Eigenvalues | 3.112 | 2.249 | 1.834 |
% of Variance | 38.900 | 28.107 | 22.919 |
Cumulative % | 38.900 | 67.008 | 89.927 |
Component matrix . | Component 1 . | Component 2 . | Component 3 . |
---|---|---|---|
MAP | 0.039 | 0.884 | 0.428 |
SM | −0.181 | 0.901 | 0.085 |
MAT | 0.469 | 0.126 | 0.865 |
TCM | −0.157 | 0.189 | 0.954 |
TWM | 0.962 | −0.122 | −0.189 |
GDD5 | 0.985 | −0.097 | 0.060 |
GDD0 | 0.954 | −0.038 | 0.233 |
MI | 0.296 | 0.630 | −0.534 |
Α | 0.781 | 0.424 | −0.057 |
Eigenvalues | 3.112 | 2.249 | 1.834 |
% of Variance | 38.900 | 28.107 | 22.919 |
Cumulative % | 38.900 | 67.008 | 89.927 |
Partial CCAs (pCCA) were conducted to partition the relative importance of different environmental variables that resulted in variation in species composition. The ‘Var-part-3 groups-Conditional-effects-tested’ method in Canoco 5 was used for this purpose. There were three types of environmental variables in this dataset: the seven bioclimatic variables mentioned above, the geographic coordinates of each site and the topographic variables. Given that the number of variables in each group and therefore the degree of freedom of each group may substantially affect the explained proportion, we performed a set of pCCAs for subsets of the original dataset. In each of the subsets, the number of independent variables was the same. Because the second and third environmental group contained three variables, we randomly selected three climatic variables from the total, which produced 35 subsets of the original dataset. Thus, the environmental data used for the subsequent 35 pCCAs shared the same data for the geographic coordinates and topography but differed with respect to climatic variables. The results from individual pCCAs were summarized using Post Hoc multiple comparison tests, selecting Tamhane’s T2 method, which does not assume equal variances. The ordination analyses were performed using Canoco 5 (ter Braak and Smilauer 2012).
The dominant spruce species from the 13 alliances can be classified into morphologically defined sections or phylogenetically distinct clades. Accordingly, the samples from the 13 alliances were reclassified into two sets of groups in the biplot of the CCA to show the relationship between the morphological and phylogenetic divergence of the dominant species and climatic variables. One-way ANOVA was conducted to examine the differences in climatic variables among the clades and among the sections. The phylogenetic clades of spruce are largely equivalent in different literature but are named differently (see Supplementary Appendix S3). We used the system presented by Ran et al. (2006) for simplicity.
RESULTS
Vegetation classification
The TWINSPAN resulted in 13 alliances when the analysis was terminated at the 12th division level. Alli. P. morrisonicola was identified at the first division, and the remaining alliances were identified at the subsequent 10 division levels. Specifically, from division 2 to 12, the alliances were Alli. P. jezoensis var. komarovii and Alli. P. koraiensis (division 3), Alli. P. obovata and Alli. P. schrenkiana (division 5), Alli. P. crassifolia and Alli. P. meyeri (division 7), Alli. P. brachytyla and Alli. P. likiangensis (division 9), Alli. P. likiangensis var. rubescens (division 10), Alli. P. asperata (division 11), and Alli. P. wilsonii and Alli. P. purpurea (division 12) (Fig. 2).

A dendrogram showing the TWINSPAN classification for the samples based on the level of alliance. The sections (◆ Sect. Picea, ▽ Sect. Omorica, ◇ Sect. Casicta) (Zheng and Fu 1978) and three phylogenetic clades of Picea distributed in Asia (Ran et al. 2006), to which the dominant spruce species in question belong, are listed and marked.
Variation in climatic variables
The PCA of the climatic variables across the sampling sites revealed three dominant climatic gradients that accounted for 89.9% of the variance. The first component, which had an eigenvalue of 3.112 and accounted for 38.9% of the variance, was a gradient characterized by variation in the cumulative growing season temperature (GDD5 and GDD0) and the mean temperature in the warmest month (TWM). The second and third components, which accounted for 28.17% and 22.92% of the variance, were characterized by variation in moisture (MAP and SM) and the mean temperature in the coldest month (TCM), respectively (Table 1).
Variation in species composition
The 172 samples contained 1181 vascular plant species belonging to 406 genera and 102 families. The DCA explained 5.9% of the total variation (43.27) in species composition. Given the large number of species in these datasets, the explained percentage should be reasonable. The eigenvalues of the first four DCA axes were 0.953, 0.714, 0.620 and 0.554, and their gradient lengths were 15.11, 6.76, 6.84 and 5.02, respectively, suggesting strong variation (beta diversity) of species composition between samples. A unimodal response model such as the CCA was therefore an appropriate choice for the following constrained ordination analyses.
Relationship between species composition and environmental variables
The CCA of species versus the set of environmental variables showed rather diverse variations, which were indicated by the relatively high eigenvalues of the first four axes. The variation in species data explained by the first four axes was relatively low, but nearly half of the variations in the species—environment relationship were explained by the first two axes. The permutation test indicated that the species—environment relationships were significant for all axes (P = 0.002) (Table 2).
Axis . | 1 . | 2 . | 3 . | 4 . |
---|---|---|---|---|
Eigenvalues | 0.878 | 0.844 | 0.619 | 0.533 |
Species–environment correlation (pseudo-canonical correlation) | 0.994 | 0.983 | 0.920 | 0.921 |
Cumulative % variation of species data | 2.03 | 3.98 | 5.41 | 6.64 |
Cumulative % variation of species–environment relations | 23.83 | 46.74 | 63.54 | 78.01 |
Significance (permutation test on all axes): pseudo-F = 2.5, P = 0.002 |
Axis . | 1 . | 2 . | 3 . | 4 . |
---|---|---|---|---|
Eigenvalues | 0.878 | 0.844 | 0.619 | 0.533 |
Species–environment correlation (pseudo-canonical correlation) | 0.994 | 0.983 | 0.920 | 0.921 |
Cumulative % variation of species data | 2.03 | 3.98 | 5.41 | 6.64 |
Cumulative % variation of species–environment relations | 23.83 | 46.74 | 63.54 | 78.01 |
Significance (permutation test on all axes): pseudo-F = 2.5, P = 0.002 |
Axis . | 1 . | 2 . | 3 . | 4 . |
---|---|---|---|---|
Eigenvalues | 0.878 | 0.844 | 0.619 | 0.533 |
Species–environment correlation (pseudo-canonical correlation) | 0.994 | 0.983 | 0.920 | 0.921 |
Cumulative % variation of species data | 2.03 | 3.98 | 5.41 | 6.64 |
Cumulative % variation of species–environment relations | 23.83 | 46.74 | 63.54 | 78.01 |
Significance (permutation test on all axes): pseudo-F = 2.5, P = 0.002 |
Axis . | 1 . | 2 . | 3 . | 4 . |
---|---|---|---|---|
Eigenvalues | 0.878 | 0.844 | 0.619 | 0.533 |
Species–environment correlation (pseudo-canonical correlation) | 0.994 | 0.983 | 0.920 | 0.921 |
Cumulative % variation of species data | 2.03 | 3.98 | 5.41 | 6.64 |
Cumulative % variation of species–environment relations | 23.83 | 46.74 | 63.54 | 78.01 |
Significance (permutation test on all axes): pseudo-F = 2.5, P = 0.002 |
As shown by the ‘sample-environment’ biplot diagram (Fig. 3), the first axis defined a temperature gradient varying from warm winters and cool summers to cool winters and warm summers. The second axis, however, was a moisture gradient with the wet end at the bottom and the dry end at the top of the biplot.

A CCA biplot defined by the first two axes showing the relationship between the 13 spruce alliances and environmental variables (climatic variables and geographical location).
Along the temperature axis, most of the spruce alliances were scattered on the left half of the biplot and only three occurred on the right, suggesting that the habitats characterized by warm winters and cool summers in high elevation areas tend to be favored by most of the spruce communities. Along the moisture axis, most of the spruce alliances occurred in habitats with moderate precipitation (the mid-part of the biplot). Alli. P. schrenkiana and Alli. P. obovata were two exceptions that were scattered on the dry end, while Alli. P. morrisonicola occurred on the wet end. In general, a cool and wet climate was favored by most of the spruce alliances.
Partitioning the relative importance of climate and geographical location for species composition
At the given degree of freedom (i.e. df = 3), the total explainable variation explained by the geographical coordinates (32.01%) was significantly higher than that explained by the climatic (27.76%) and topographic variables (23.32%). In addition, the combined effects of the three types of environmental variables were fairly low, ranging between 7.12% and 1.10% (Fig. 4, Supplementary Appendix S4).

Partitioning of the variation in species composition explained by the three groups of environmental variables (G, C and T) and their combined effects. The differences between groups were significant (P < 0.001) for all comparisons. Gs, Cs and Ts represent the variation explained by G, C and T, respectively. G∩C is the variation explained by G and C. G∩C∩T is the variation explained by G, C and T. G∩T is the variation explained by G and T, and C∩T is the variation explained by C and T. G = geographical coordinates; C = climatic variables; T = topographic variables.
Relationship between the morphological and phylogenetic divergence of the dominant species and climatic variables
In the same ordination space as that presented in Fig. 3, sect. Omorica occurred mainly in the wet habitats, and sect. Casicta tended to favor climates with moderate moisture, warm winters, and cool summers, yet sect. Picea covered a wide range in the mid-to-dry habitats (Fig. 5). Specifically, the mean annual precipitation for the habitats of sect. Omorica was ca. 229 mm and 426 mm higher than that for sect. Casicta and sect. Picea (Fig. 7b, P < 0.01). Differences in the temperature variables among the three sections were less distinct (Fig. 7d and f).

A CCA biplot defined by the first two axes showing the relationship between the three sections (Zheng and Fu 1978) and environmental variables (climatic variables and geographical location). The drawings inserted in the biplot indicate the leaf cross-section of the three sections.
When grouped into phylogenetically distinct clades (Ran et al. 2006), spruce species from Clade-II A and Clade-II B (two older sub-clades) occurred in habitats with warm winters and cool summers and in areas with mid to higher elevation, while those from Clade-III (a younger clade) occurred in habitats with cold winters and warm summers in mid to low elevation areas, extending across a wide moisture gradient (Fig. 6). Specifically, for the habitats of the two older sub-clades, the mean temperature in the coldest month was ca. 8–10°C higher, yet the accumulated temperature during the growing season (≥5°C) was ca. 297–438°C lower than the habitats of Clade-III (Fig. 7c and e, P < 0.01). Differences in MAP among the clades were less distinct (Fig. 7a).

A CCA biplot defined by the first two axes showing the relationship between the two phylogenetically distinct clades (Ran et al. 2006) and environmental variables (climatic variables and geographical location). A simplified cladogram of Picea was inserted in the biplot (Ran et al. 2006).

Box plots showing the differences in climatic variables among the phylogenetic clades (Ran et al. 2006) and among the three morphological sections of spruce species in China. The central box in each box plot indicates the interquartile range and median, while whiskers show the 10th and 90th percentiles. Mean values marked with different letters indicate a significant difference at P < 0.01.
DISCUSSION
Species composition of spruce forests is closely related to climate and geographical location
Species composition of a community is assumed to be largely determined by the size of the local species pool, which in turn is determined by the regional pool, while climate determines various species pools (Lintulaakso et al. 2019; Pärtel et al. 1996; Zobel et al. 1998). It is therefore reasonable to predict that those plant communities that are geographically closely located, and hence climatically similar, should share more species in common. This prediction is supported by the results from the present study.
A hierarchical vegetation classification resulting from TWINSPAN is ultimately determined by the dissimilarity of species composition between plant communities (Hill and Šmilauer 2005). Our results indicated that the dissimilarity of species composition between spruce forests was strongly related to their geographical locations. For example, Alli. P. morrisonicola, whose natural distribution is in the subalpine area of Taiwan, which is geographically isolated from the Chinese mainland by the Taiwan Strait, was identified at the first division. For the remaining 12 spruce forests distributed on the main land, the alliances that were closely geographically located tended to share more species in common compared with those that were geographically distant. The results from the CCA further showed that the 13 spruce alliances were climatically distinct from one another, which suggests that species composition between different spruce forests was closely related to their climatic differentiation. Specifically, the higher the similarity in species composition, the closer the distance between spruce forests is in the climatically determined ordination space.
While the close relationship between climate or geographical location and the distribution of spruce forests is rather obvious, the variation in species composition explained by climatic and geographical variables was fairly low. This finding can be explained by the following. First, the total number of species in our dataset was relatively high (>1000), suggesting high intraspecific variation in plant characteristics, which are assumed to affect plant competition, propagation and dispersal and will certainly influence the species composition of a community (Zobel et al. 1998, 2005). Second, the measured abiotic variables in this study might account for a small proportion of the total environmental variables that determine species composition (Økland 1999). Specifically, edaphic factors and geographical barriers, which were not used as explanatory variables in this analysis, might play important roles in determining species composition (Augusto et al. 2003; Just 1947; Qian et al. 2003).
Geographical location tends to be a more powerful variable that affects species composition compared with climatic variables
While climatic variables and geographical location are mutually correlated, the latter may have a stronger influence on plant distribution (Huo et al. 2014). Geographical location can be used as a proxy for local or regional climatic conditions and can indicate the geographical constraints that influence post-glacial plant expansion and dispersal (Svenning and Skov 2004). It is well known that plant distribution or species composition in refugia is likely influenced by the vicissitudes of past geology and climate but has little to do with current climate (Araújo et al. 2005; Chou et al. 2011; Stewart et al. 2010). In eastern Asia, the uplift of the Tibetan Plateau during the late Tertiary and the Quaternary resulted in many post-glacial refugia (Hampe and Petit 2010; Wu and Wang 1980; Xing et al. 2017). In addition, the Taiwan Strait, which is a major geographical barrier, might have prevented the exchange of flora between Taiwan and the mainland, leading to a relatively high proportion of endemic taxa in the flora of Taiwan (Wu and Wang 1980; Yang et al. 2018). Thus, it is understandable that geographical locations would have greater potential than climatic variables for explaining plant distribution, which holds true especially for the areas where post-glacial refugia are fairly common.
In this case, the present distribution of P. morrisonicola may be defined as a refugium of the post-glacial period because of the decline of its natural range as a result of climate warming during the post-glacial period rather than cooling during the glacial period (Tsukada 1966; Xu et al. 1973, 1980). Taiwan spruce has been identified as a vulnerable species by the ‘IUCN Red List of Threatened Species’ because of its declining population and reduced occupancy area in high elevation zones (Zhang et al. 2013). In our field survey, we found that Taiwan spruce forests grow mainly in relatively closed microhabitats in high elevation zones under cool but humid conditions. However, these types of habitats are relatively scarce on the island.
Morphological divergence among the spruce species is largely determined by moisture conditions
The shape of leaf cross-sections and the position of stomates on the leaf surface are two prominent and key characteristics that have been used to divide spruce species into subdivisions at the first split in Flora of North America (Taylor 1993), and in Flora of China (Fu et al 1999; Zheng and Fu 1978). Generally, spruce with flattened or broad triangular leaves whose stomates occur on the adaxial leaf surface tend to have a narrow natural range, while those with quadrangular leaves whose stomates occur more or less equally on all leaf surfaces are naturally widespread (Fu et al 1999; Taylor 1993; Zheng and Fu 1978). Our results indicate that spruce with quadrangular leaves (sect. Picea) tend to be distributed in drier habitats, while those with flattened or broad triangular leaves (sect. Omorica and sect. Casicta) occur mainly in wetter climates, which is consistent with a previous study on the morphological divergence of spruce species worldwide (Wang et al. 2017).
Taiwan spruce (i.e. P. morrisonicola) tends to be an exception to the above-mentioned pattern. Specifically, Taiwan spruce with quadrangular leaves belongs to sect. Picea but is currently distributed in subtropical areas with high rainfall. This incongruence may arise from a more complex process, that is, the frequent expansion and retreat of spruce in Asia that was associated with severe climatic oscillation during the glacial cycles in the late Pliocene and the Pleistocene (Xu et al. 1980). The last glacial maximum (1.8–2.0 myr BP) during the Quaternary might have a particular influence on these processes (Sommer and Zachos 2009; Willis and Niklas 2004).
It is well known that while tropical and subtropical forests survived in refugia or became extinct during the Quaternary glaciations, temperate and boreal forests expanded to low latitude or altitude areas (Hampe and Petit 2010; Willis and Niklas 2004). Thus, a higher proportion of pollen from boreal components, such as spruce in specific sediments, is generally assumed to indicate a cold period, and vice versa (Xu et al. 1980). The earliest spruce pollen fossil in Asia was found on the Tibetan Plateau in the late Oligocene to the early Miocene (Wu et al. 2007). Since then, spruce pollen has frequently been found in sediments of the late Pliocene and the Pleistocene in northern, eastern and southwestern China (Shi 1996; Xu et al. 1973, 1980), and on Taiwan island (Tsukada 1966). In addition, the proportion of spruce pollen in the sediments varied substantially with their geological ages, suggesting that spruce was widespread across eastern Asia but underwent frequent expansion and retreat during the glacial cycles in the past. Meanwhile, Taiwan island was intermittently connected to the mainland during the cold periods of the last glacial maxima due to the decrease in sea level, which ensured the exchange of flora between the mainland and Taiwan island (Chou et al. 2011; Tsukada 1966).
However, as the coldest period during the last glacial maximum, the Tali glacial period (50–60 kyr BP) likely resulted in an 8.0–11.0°C reduction in mean annual temperature in Taiwan, which was further intensified by a cold coastal current. Consequently, the climate must have been more continental with a strong winter monsoon (Tsukada 1966). Such climatic conditions may be more suitable for spruce from sect. Picea, which can withstand cold and dry habitats, but less so for spruce species with flattened leaves, which tend to be adapted to warm and humid climates. We therefore deduce that spruce with flattened leaves might not have been dispersed across Taiwan during the glacial maxima. In addition, the appearance of the Taiwan Strait in the post-glacial period might have prevented the exchange of flora between Taiwan and the mainland, which did not allow the post-glacial expansion of spruce into Taiwan. Consequently, spruce species with flattened leaves are totally absent from spruce forests in Taiwan.
Despite the incongruence between the climate and the morphological characteristics of P. morrisonicola, this species should be among the most warm-humid adapted species in sect. Picea. Interestingly, P. wilsonii, which is a phylogenetically closely related species to P. morrisonicola in sect. Picea (Lockwood et al. 2013; Ran et al. 2006), is also a warm-humid adapted species that is distributed across the Chinese mainland.
Phylogenetic divergence among spruce species is mainly influenced by temperature conditions
Molecular phylogenetic studies show that spruce species worldwide can be grouped into three clades, two of which are naturally distributed in China (Lockwood et al. 2013; Ran et al. 2006). Although the origin of spruce remains controversial, a consensus on the phylogenetic divergence among the spruce species in China has been reached. For example (Supplementary Appendix S3), Clade-II and Clade-III in Ran et al. (2006) are equivalent to Clade-III and Clade I in Lockwood et al. (2013). This provides a good base for further consideration of the factors driving their divergence.
The results of this study show that temperature tends to be the first determinant that accounts for the divergence among phylogenetically defined clades which is in line with Wang et al. (2017). These findings suggest that temperature could play a critical role in the context of spruce evolution and speciation, although moisture conditions are not negligible. This inference is consistent with the fact that temperature was the most prominent variable influenced by the glacial cycles (Xu et al. 1980). Specifically, the divergence among the three spruce clades might have occurred in a period of severe oscillation of temperature during the late Oligocene or the early Miocene (cf. Lockwood et al. 2013).
Although the origin of spruce dates back to the early Tertiary or the late Cretaceous, most spruce species in Asia may have arisen more recently (Ran et al. 2006; Sigurgeirsson and Szmidt 1993). The Himalayan-Hengduan Mountains, that is, the area that is overlapped by the three sections and the clades in the CCA ordination diagrams (Figs 3, 5, and 6), is likely the first destination where spruce migrated from high latitude areas when the world became colder since the Miocene (Ran et al. 2006). Our results clearly show that this region is a center of spruce diversity and spruce morphological and phylogenetic divergence. Given the incongruence between the climate and the morphological characteristics of P. morrisonicola, its present distribution is likely a refugium of the post-glacial period.
Morphological convergence among phylogenetically distinct spruce species might arise from parallel evolution, that is, the repeated appearance of similar characteristics that occur among distantly related species, which is fairly common in nature (Hoekstra and Price 2004; Orr 2005; Schluter et al. 2004; Went 1971). In this case, P. schrenkiana (Clade-II A) and P. obovata (Clade-III), which are two phylogenetically distinct taxa with quadrangular leaves, were likely selected by a relatively dry climate, while the flattened leaves of P. brachytyla (Clade-II B) and P. jezoensis (Clade-III) were possibly selected by a relatively wet climate.
While both morphological sections and phylogenetic clades of the spruce species in question carry strong climatic signals, their divergence is subject to different selective pressures.
Supplementary Material
Supplementary material is available at Journal of Plant Ecology online.
Appendix S1: The morphological traits and species for the three sections of spruce in China identified in Zheng and Fu (1978).
Appendix S2: Supplementary methods.
Appendix S3: The phylogenetic clades of spruce in China identified in Ran et al. (2006) and in Lockwood et al. (2013).
Appendix S4: A summary of the 35 analyses for the partitioning of the variation in species composition explained by the three groups of environmental variables (G, C and T) as well as their joint effects.
Funding
This work was supported by National Natural Science Foundation of China (41571045).
Acknowledgements
We thank Lijiang Zhou, Hongchun Wang, Zhi Ma, Ziying Chen and Tiancai Chen for field assistance.
Conflict of interest statement. None declared.
REFERENCES
Author notes
These authors contributed equally to this work.