Abstract

Aims

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.

Methods

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.

Important Findings

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.
Figure 1:

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).

Table 1:

PCA of climate variables used in the factor analysis, showing the correlations of climate variables with the three components (1, 2 and 3)

Component matrixComponent 1Component 2Component 3
MAP0.0390.8840.428
SM−0.1810.9010.085
MAT0.4690.1260.865
TCM−0.1570.1890.954
TWM0.962−0.122−0.189
GDD50.985−0.0970.060
GDD00.954−0.0380.233
MI0.2960.630−0.534
Α0.7810.424−0.057
Eigenvalues3.1122.2491.834
% of Variance38.90028.10722.919
Cumulative %38.90067.00889.927
Component matrixComponent 1Component 2Component 3
MAP0.0390.8840.428
SM−0.1810.9010.085
MAT0.4690.1260.865
TCM−0.1570.1890.954
TWM0.962−0.122−0.189
GDD50.985−0.0970.060
GDD00.954−0.0380.233
MI0.2960.630−0.534
Α0.7810.424−0.057
Eigenvalues3.1122.2491.834
% of Variance38.90028.10722.919
Cumulative %38.90067.00889.927
Table 1:

PCA of climate variables used in the factor analysis, showing the correlations of climate variables with the three components (1, 2 and 3)

Component matrixComponent 1Component 2Component 3
MAP0.0390.8840.428
SM−0.1810.9010.085
MAT0.4690.1260.865
TCM−0.1570.1890.954
TWM0.962−0.122−0.189
GDD50.985−0.0970.060
GDD00.954−0.0380.233
MI0.2960.630−0.534
Α0.7810.424−0.057
Eigenvalues3.1122.2491.834
% of Variance38.90028.10722.919
Cumulative %38.90067.00889.927
Component matrixComponent 1Component 2Component 3
MAP0.0390.8840.428
SM−0.1810.9010.085
MAT0.4690.1260.865
TCM−0.1570.1890.954
TWM0.962−0.122−0.189
GDD50.985−0.0970.060
GDD00.954−0.0380.233
MI0.2960.630−0.534
Α0.7810.424−0.057
Eigenvalues3.1122.2491.834
% of Variance38.90028.10722.919
Cumulative %38.90067.00889.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.
Figure 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).

Table 2:

CCA of species composition in relation to environmental variables

Axis1234
Eigenvalues0.8780.8440.6190.533
Species–environment correlation (pseudo-canonical correlation)0.9940.9830.9200.921
Cumulative % variation of species data2.033.985.416.64
Cumulative % variation of species–environment relations 23.8346.7463.5478.01
Significance (permutation test on all axes): pseudo-F = 2.5, P = 0.002
Axis1234
Eigenvalues0.8780.8440.6190.533
Species–environment correlation (pseudo-canonical correlation)0.9940.9830.9200.921
Cumulative % variation of species data2.033.985.416.64
Cumulative % variation of species–environment relations 23.8346.7463.5478.01
Significance (permutation test on all axes): pseudo-F = 2.5, P = 0.002
Table 2:

CCA of species composition in relation to environmental variables

Axis1234
Eigenvalues0.8780.8440.6190.533
Species–environment correlation (pseudo-canonical correlation)0.9940.9830.9200.921
Cumulative % variation of species data2.033.985.416.64
Cumulative % variation of species–environment relations 23.8346.7463.5478.01
Significance (permutation test on all axes): pseudo-F = 2.5, P = 0.002
Axis1234
Eigenvalues0.8780.8440.6190.533
Species–environment correlation (pseudo-canonical correlation)0.9940.9830.9200.921
Cumulative % variation of species data2.033.985.416.64
Cumulative % variation of species–environment relations 23.8346.7463.5478.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).
Figure 3:

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.
Figure 4:

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.
Figure 5:

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).
Figure 6:

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.
Figure 7:

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

Araújo
MB
,
Pearson
RG
,
Rahbek
C
(
2005
)
Equilibrium of species’ distributions with climate
.
Ecography
28
:
693
5
.

Augusto
L
,
Dupouey
JL
,
Ranger
J
(
2003
)
Effects of tree species on understory vegetation and environmental conditions in temperate forests
.
Ann For Sci
60
:
823
31
.

Bailey
SF
,
Rodrigue
N
,
Kassen
R
(
2015
)
The effect of selection environment on the probability of parallel evolution
.
Mol Biol Evol
32
:
1436
48
.

ter Braak
CJF
(
1986
)
Canonical correspondence analysis: a new eigenvector technique for multivariate direct gradient analysis
.
Ecology
67
:
1167
79
.

ter Braak
CJF
,
Smilauer
P
(
2012
)
CANOCO Reference Manual and User’s Guide: Software for Ordination, Version 5
.
New York, NY
:
Microcomputer Power
.

Chou
YW
,
Thomas
PI
,
Ge
XJ
, et al. (
2011
)
Refugia and phylogeography of “Taiwania” in East Asia
.
J Biogeogr
38
:
1992
2005
.

Farjón
A
(
2001
)
World Checklist and Bibliography of Conifers, 2nd edition
.
London, UK
:
Royal Botanical Garden, Kew
.

Fu
L
,
Li
N
,
Mill
RR
(
1999
)
Picea.
In
Wu
ZY
,
Raven
PHE
(eds).
Flora of China
.
Beijing & St. Louis
:
Science Press & Missouri Botanical Garden Press
,
25
32
.

Gauthier
S
,
Grandpré
LD
,
Bergeron
Y
(
2000
)
Differences in forest composition in two boreal forest ecoregions of Quebec
.
J Veg Sci
11
:
781
90
.

Hampe
A
,
Petit
RJ
(
2010
)
Cryptic forest refugia on the ‘Roof of the World’
.
New Phytol
185
:
5
7
.

Haxeltine
A
,
Prentice
IC
(
1996
)
BIOME3: an equilibrium terrestrial biosphere model based on ecophysiological constraints, resource availability, and competition among plant functional types
.
Glob Biogeochem Cycle
10
:
693
709
.

Hewitt
G
(
2000
)
The genetic legacy of the Quaternary ice ages
.
Nature
405
:
907
13
.

Hill
MO
,
Jr
HGG
(
1980
)
Detrended correspondence analysis: an improved ordination technique
.
Vegetatio
42
:
47
58
.

Hill
MO
,
Šmilauer
P
(
2005
)
TWINSPAN for Windows Version 2.3
.
Huntingdon & Ceske Budejovice
:
Centre for Ecology and Hydrology & University of South Bohemia
.

Hoekstra
HE
,
Price
T
(
2004
)
Parallel evolution is in the genes
.
Science
303
:
1779
81
.

Huo
H
,
Feng
Q
,
Su
YH
(
2014
)
The influences of canopy species and topographic variables on understory species diversity and composition in coniferous forests
.
Sci World J
2014
:
252489
.

Hutchinson
MF
(
1989
)
A new objective method for spatial interpolation of meteorological variables from irregular networks applied to the estimation of monthly mean solar radiation, temperature, precipitation and windrun
.
CSIRO Division of Water Resources Tech. Memo
89
:
95
104
.

Just
T
(
1947
)
Geology and plant distribution
.
Ecol Monogr
17
:
127
37
.

LePage
BA
(
2001
)
New species of Picea A. Dietrich (Pinaceae) from the middle Eocene of Axel Heiberg Island, Arctic Canada
.
Biol J Linnean Soc
135
:
137
67
.

Liang
Q
,
Xu
X
,
Mao
K
, et al. (
2018
)
Shifts in plant distributions in response to climate warming in a biodiversity hotspot, the Hengduan Mountains
.
J Biogeogr
45
:
1334
44
.

Lintulaakso
K
,
Polly
PD
,
Eronen
JT
(
2019
)
Land mammals form eight functionally and climatically distinct faunas in North America but only one in Europe
.
J Biogeogr
46
:
185
95
.

Lockwood
JD
,
Aleksić
JM
,
Zou
J
, et al. (
2013
)
A new phylogeny for the genus Picea from plastid, mitochondrial, and nuclear sequences
.
Mol Phylogenet Evol
69
:
717
27
.

Økland
RH
(
1999
)
On the variation explained by ordination and constrained ordination axes
.
J Veg Sci
10
:
131
6
.

Økland
RH
,
Rydgren
K
,
Økland
T
(
2003
)
Plant species composition of boreal spruce swamp forests: closed doors and windows of opportunity
.
Ecology
84
:
1909
19
.

Orr
HA
(
2005
)
The probability of parallel evolution
.
Evolution
59
:
216
20
.

Pärtel
M
,
Zobel
M
,
Zobel
K
, et al. (
1996
)
The species pool and its relation to species richness: evidence from Estonian plant communities
.
Oikos
75
:
111
7
.

Qian
H
,
Klinka
K
,
Økland
RH
, et al. (
2003
)
Understorey vegetation in boreal Picea mariana and Populus tremuloides stands in British Columbia
.
J Veg Sci
14
:
173
84
.

Ran
JH
,
Wei
XX
,
Wang
XQ
(
2006
)
Molecular phylogeny and biogeography of Picea (Pinaceae): implications for phylogeographical studies using cytoplasmic haplotypes
.
Mol Phylogenet Evol
41
:
405
19
.

Schluter
D
,
Clifford
E
,
Nemethy
M
, et al. (
2004
)
Parallel evolution and inheritance of quantitative traits
.
Am Nat
163
:
809
22
.

Shi
N
(
1996
)
Development of spruce and fir in north China during the Pliocnen and the early Pleistonene: palaeoclimatic implications
.
Quaternary Sci
16
:
319
28
.

Sigurgeirsson
A
,
Szmidt
AE
(
1993
)
Phylogenetic and biogeographic implications of chloroplast DNA variation in Picea
.
Nordic J Bot
13
:
233
46
.

Sommer
RS
,
Zachos
FE
(
2009
)
Fossil evidence and phylogeography of temperate species: ‘glacial refugia’ and post-glacial recolonization
.
J Biogeogr
36
:
2013
20
.

Stewart
JR
,
Lister
AM
,
Barnes
I
, et al. (
2010
)
Refugia revisited: individualistic responses of species in space and time
.
Proc R Soc B-Biol Sci
277
:
661
71
.

Sun
H
,
Li
ZM
(
2003
)
Qinghai-Tibetan Plateau uplift and its impact on Tethys flora
.
Advance in Earth Sciences
18
:
852
62
.

Svenning
JC
,
Skov
F
(
2004
)
Limited filling of the potential range in European tree species
.
Ecol Lett
7
:
565
73
.

Taylor
RJ
(
1993
)
Picea.
In
Flora of North America Editorial Committee
(ed).
Flora of North America North of Mexico
.
New York, NY
:
Oxford University Press
.

Toby
S
,
Chytrý
M
(
2002
)
Vegetation surveys in the circumboreal coniferous forests: a review
.
Folia Geobot
37
:
365
82
.

Transeau
EN
(
1909
)
The relation of the climatic factors to vegetation
.
Am Nat
43
:
487
93
.

Tseng
YS
(
1991
)
Study on the vegetation ecology of Salixianxi watershed in central Taiwan
. Master Thesis.
China
:
Taiwan University
.

Tsukada
M
(
1966
)
Late Pleistocene vegetation and climate in Taiwan (Formosa)
.
Proc Natl Acad Sci USA
55
:
543
8
.

Walker
D
(
1986
)
Late Pleistocene-early Holocene vegetational and climatic changes in Yunnan Province, southwest China
.
J Biogeogr
13
:
477
86
.

Wang
GH
,
Li
H
,
Zhao
HW
, et al. (
2017
)
Detecting climatically driven phylogenetic and morphological divergence among spruce (Picea) species worldwide
.
Biogeosci
14
:
2307
19
.

Went
FW
(
1971
)
Parallel evolution
.
Taxon
20
:
197
226
.

Willis
KJ
,
Niklas
KJ
(
2004
)
The role of Quaternary environmental change in plant macroevolution: the exception or the rule?
Philos Trans R Soc Lond Ser B-Biol Sci
359
:
159
72
.

Woodward
FI
,
Lomas
MR
,
Kelly
CK
(
2004
)
Global climate and the distribution of plant biomes
.
Philos Trans R Soc Lond Ser B-Biol Sci
359
:
1465
76
.

Woodward
FI
,
Williams
BG
(
1987
)
Climate and plant distribution at global and local scales
.
Vegetatio
69
:
189
97
.

Wu
ZY
,
Wang
HS
(
1980
)
Floristic characteristics of vegetation of China.
In
Wu
ZY
(ed).
Vegetation of China
.
Beijing
:
Science Press
,
82
140
.

Wu
Z
,
Wu
Z
,
Hu
D
, et al. (
2007
)
Geological evidences for the Tibetan Plateau uplifted in late Oligocene
.
Acta Geologica Sinica
81
:
577
87
.

Xing
Y
,
Ree
RH
(
2017
)
Uplift-driven diversification in the Hengduan Mountains, a temperate biodiversity hotspot
.
Proc Natl Acad Sci USA
114
:
E3444
51
.

Xiong
Y
,
Li
QK
(
1987
)
Soil of China
.
Beijing
:
Science Press
.

Xu
R
,
Kong
ZC
,
Du
NQ
(
1980
)
Plant assemblages of Picea and Abies in the Pleistocene and implications for Quarternary study
.
Quaternary Sci
5
:
48
56
.

Xu
R
,
Tao
JR
,
Sun
XJ
(
1973
)
On the discovery of a Quercus semicarpifolia bed in Mount Shisha Pangma and its significance in botany and geology
.
Acta Botanica Sinica
15
:
103
19
.

Yang
Y
,
Chen
L
,
Peng
H
(
2018
)
Clarification of the taxonomic confusion between Ilex formosana and I. tetramera (Aquifoliaceae), with the description of a new species I. shukunii
.
Phytotaxa
382
:
182
92
.

Yang
GZ
,
Chen
YG
,
Zhao
WC
, et al. (
2002
)
Plant Resource Investigation in Nantzuhsien Creek Watershed in Yushan of National Park
.
Nantou
:
Yushan of National Park administrative office
.

Zhang
XS
(
2007
)
Vegetation Map of the People’s Republic of China (1: 1 000 000)
.
Beijing, China
:
Geological Publishing House
.

Zhang
D
,
Katsuki
T
,
Rushforth
K
(
2013
)
Picea morrisonicola
. In
The IUCN Red List of Threatened Species 2013: e.T34383A2852220
. (25 October 2018, date last accessed).

Zhao
H
,
Guo
K
,
Yang
Y
, et al. (
2018
)
Stipa steppes in scantily explored regions of the Tibetan Plateau: classification, community characteristics and climatic distribution patterns
.
J Plant Ecol
11
:
585
94
.

Zheng
W
,
Fu
L
(
1978
)
Flora of China
.
Beijing, China
:
Science Press
,
123
67
.

Zhou
YL
,
Yang
QX
(
1980
)
Cold temperate coniferous forests.
In
Wu
ZY
(ed).
Vegetation of China
.
Beijing, China
:
Science Press
,
191
203
.

Zobel
M
,
van der Maarel
E
,
Dupré
C
(
1998
)
Species pool: the concept, its determination and significance for community restoration
.
Appl Veg Sci
1
:
55
66
.

Zobel
M
,
Otsus
M
,
Rünk
K
, et al. (
2005
)
Can long-distance dispersal shape the local and regional species pool?
Folia Geobot
40
:
35
44
.

Author notes

These authors contributed equally to this work.

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
Handling Editor: Christian Schöb
Christian Schöb
Handling Editor
Search for other works by this author on: