1. Introduction

Research on species diversity patterns along mountain elevational gradients has increased considerably over the past three decades (Colwell et al., 2004, Currie and Kerr, 2008, Manish et al., 2017, Rahbek, 1995). However, our understanding of the underlying mechanisms that drive patterns of biodiversity remains inadequate (Colwell et al., 2016; Wu and Lei, 2013; Zhang and Sun, 2009). One of the main drivers of biodiversity along elevational gradients is spatial (Stiegel et al., 2011; Wu and Lei, 2013; Zhang and Sun, 2009). For example, the mid-domain effect, which has been widely investigated in several regions and taxa (Brehm et al., 2007; Colwell et al., 2016; Colwell and Lees, 2000; McCain, 2004, 2009, 2010; Rahbek, 1995, 2005; Storch et al., 2006), predicts that species richness along an elevational gradient will resemble a unimodal curve with the highest species richness at a middle elevation (Currie and Kerr, 2008). Another spatial factor that affects species richness is area. The area of different elevations varies, which likely drives patterns of species richness; for example, larger areas may harbor more species (Tang et al., 2009; Zhang and Sun, 2009).

It is widely accepted that species richness alone is insufficient to reveal the underlying mechanisms that shape patterns of biodiversity along an elevational gradient (Dı́az and Cabido, 2001). An alternative approach is to consider functional traits and functional diversity (Violle et al., 2007). Functional diversity is based on the functional traits and could provide us a new door to discover biodiversity patterns (Villéger et al., 2008). Functional traits can determine species niches and reflect the influence of resource acquisition and environmental tolerances (Bässler et al., 2016; Dı́az and Cabido, 2001). Although the main (and often only) measure of diversity along elevational gradients is still species richness (Currie and Kerr, 2008; Rahbek, 1995), increasing number of studies have used functional diversity to measure biodiversity (Hanz et al., 2019; Hao et al., 2018; Lopez et al., 2018; Matsuoka et al., 2019). Although the two related facets of species diversity and functional diversity may provide complementary insight into our overall understanding of biodiversity patterns along elevational gradients (Dı́az and Cabido, 2001), patterns of biodiversity are scale-dependent, regardless of whether species diversity or functional diversity is measured (Bhatta et al., 2018; Calba et al., 2014; Jarzyna and Jetz, 2018). Different phenomena act at different scales (or due to a researcher's subjective choice of scales) (Epstein et al., 2018; Higgins, 2010; Pegg and Taylor, 2007). Due to data collection restrictions and field survey boundedness, few studies have compiled data on traits from publications at regional scale (Epstein et al., 2018; Shiono et al., 2015; Šímová et al., 2018; Swenson et al., 2012) and combined species richness pattern with functional diversity pattern (Bässler et al., 2016, Higgins, 2010, Kessler et al., 2009).

The Hengduan Mountains system is one of the world's biodiversity hotspots (Myers et al., 2000) and an important center of speciation and distribution of alpine taxa. In addition, the Hengduan Mountains are a natural laboratory for studying the factors that drive patterns of biodiversity (Ding et al., 2020; Xing and Ree, 2017). In the Hengduan Mountains many climatic variables change faster (e.g., temperature) along elevational gradients in a relatively short horizontal distance than in flat low land (Zhang et al., 1997; Zhang, 1998). The alpine plants of the Hengduan Mountains have evolved a variety of special adaptations (such as cushion plants) to their extreme environment. Studies on alpine taxa are crucial for understanding biodiversity (Yang et al., 2019). In this region, studies have examined patterns of elevational species richness in animals (Fu et al., 2006, 2007; Wu et al., 2013a, 2013b) and plants (Feng et al., 2006a, 2006b; Shen et al., 2004; Wang et al., 2007; Yuan et al., 2008; Zhang et al., 2009a, 2009b). However, few studies have examined species richness patterns of a specific plant taxon in the Hengduan Mountains (Zhang et al., 2009b). Furthermore, few studies have analyzed species traits and functional diversity while focusing on a specific taxon at a large scale.

Delphinium species are valuable medicinal and economic angiosperms (Chen et al., 2020; Wang, 1981; Yin et al., 2020). Exploring large-scale patterns of species richness in Ranunculaceae may provide data and theoretical support for species diversity protection policies.

In this study, we mined data from the literature on species from the genus Delphinium to: 1) analyze elevational patterns of species richness, functional diversity and trait occurrence at a large scale in the Hengduan Mountains; 2) identify relationships between species diversity and functional diversity across the same large-scale elevational gradient.

2. Material and methods 2.1. Study area

The Hengduan Mountains range from 24°40′N to 34°00′N and 96°20′E to 104°30′E (Fig. 1a) (Li, 1987). This mountain region has a complex geological topography, with high vegetation and landscape diversity covering the subtropical to alpine zone (Li, 1989). The Hengduan Mountains also have extremely rich and complex species diversity, harboring an estimated 12,800 species, about 42.5% of the total number of Chinese vascular plants (Sun et al., 2017). The Hengduan Mountains include five administrative regions: western and northwestern Yunnan, western Sichuan, eastern Xizang (Tibet), southeastern Qinghai, and southwestern Gansu. The natural geomorphology of the Hengduan Mountains includes six major rivers that have divided the mountain range into seven major mountain systems from a south-north direction. This river-mountain system contains 92% of the area of all regions and assembles varied habitats (Li, 1989). The natural boundaries of the Hengduan Mountains contain 80 county-level regions (Li, 1987).

2.2. Data sources

We used 4836 county-level records to create a database of Delphinium distribution and trait occurrence across China. For each species, we compiled data on Latin names of species, lower elevational limit, upper elevational limit, elevational range, endemism, county-level distribution, and 15 functional traits.

County-level distribution and elevational records were collected from several online sources, including Flora of China: Volume 6 (Wang and Svetlana, 2001) (FOC, http://www.iplant.cn/foc/vol/6), Global Biodiversity Information Facility (GBIF, https://www.gbif.org/), the National Specimen Information Infrastructure (NSII, http://www.nsii.org.cn/2017/home.php), the Biodiversity of the Hengduan Mountains and adjacent areas of south-central China database (http://hengduan.huh.harvard.edu), the Chinese Virtual Herbarium (CVH, https://www.cvh.ac.cn/), and additional literature. We converted all county synonyms and species synonyms into accepted names. Species names followed the nomenclature of Flora of China and all intraspecies records were merged to species level. To ensure the accuracy of the natural distribution records, we also removed all the cultivation records. County data followed the most recent administrative information from the National Bureau of Statistics of the People's Republic of China (http://www.stats.gov.cn/). To determine whether species were endemic to the Hengduan Mountains, we identified species that have distribution records in Hengduan Mountains and species that are only distributed in or beyond a very small range of the Hengduan Mountains. The area at different elevational bands was calculated using the DEM model at 1 km resolution (https://srtm.csi.cgiar.org/).

2.2.1. Elevational data

In total, we collected elevational data for 102 species of Delphinium in the Hengduan Mountains of China (Table S1). For species that only had information on a single elevation point, we adjusted elevational distribution following previous studies (Brehm et al., 2007; CardelÚS et al., 2006; Wu et al., 2013a). We divided the elevational gradient (1000–5700 m) into 47 100-m-elevation bands; each species occupied at least one elevation band (Li et al., 2014; Zhang et al., 2009b). Here, we refer to these elevation bands as the low (1000–2566 m), middle (2567–4133 m), and high (4134–5700 m) elevational regions.

The most updated source of elevational ranges for Delphinium species is Flora of China. Although Flora of China records nearly the same elevational ranges of species endemic to Hengduan Mountains as other sources, it records a wider range for non-endemic species. The difference is mainly reflected in the low elevation areas outside the study area (Fig. S1). Thus, for non-endemic species, we firstly checked the Vascular Plants of Hengduan Mountains (VPHM). If species were not included in the VPHM, we checked for accurate elevational range data in the Hengduan Mountains in the CVH, NSII, and GBIF databases.

2.2.2. Species trait data

We selected the largest possible number of morphological traits to serve as functional traits for our calculation of functional diversity indices (Raine et al., 2018). Nine quantitative traits were selected: stem height, leaf width, inflorescence length, inflorescence flower number, pedicel length, bracteole length, sepal spur length, sepal spur base in diam, and other sepal lengths. Quantitative traits contained a lower limit value, upper limit value, and mean value (calculated from upper and lower limits). Six qualitative traits were also selected: inflorescence type, sepal color, petal color, staminode color, carpel number, and stem type. Qualitative traits were classified and coded for analysis (Fig. S2, details see Table S2). To demonstrate the patterns of six qualitative traits along the elevational gradient, we transformed them into the ratio of each level in each elevational band. For example, we used the frequency of blue-purple sepal (blue-purple species number/all species number in the same band) in one elevational band to represent the ratio of blue-purple level in sepal color. Correlation analyses were performed to examine the relationship between all morphological traits.

2.2.3. Dividing groups

We analyzed Delphinium species based on endemism and range. To determine patterns of biodiversity among species endemic to the Hengduan Mountains, we first identified 75 endemic species and 27 non-endemic species. Endemic species were defined as species distributed in or only beyond a very small range of the Hengduan Mountains; the remaining species were defined as non-endemic species. MDE null model predicted that larger-ranged species were more sensitive to geometric constrains (Colwell and Lees, 2000). Therefore, we also analyzed the effect of elevational range on patterns of biodiversity. Following the method of previous studies (Wu et al., 2013b), we divided species into small ranges species (100 m), medium ranges species (200–800 m), or large ranges species (900–3400 m). There were 39 species in the small range, 26 in the medium range, and 37 in the large range.

2.2.4. Data analysis

To calculate null species richness and upper/lower 95% confidence interval (CI) for each elevational band, we ran 5000 resampled simulations using the Colwell & Lees Box 5 simulation algorithm (continuous domain analysis) in RangeModel 5 (http://purl.oclc.org/rangemodel) (Chen et al., 2017; Wu et al., 2013a). To estimate the impact of the mid-domain effect model for each species group (i.e., endemic, non-endemic, small, medium, and large range species), we performed a linear regression between empirical richness and mid-domain effect predictions of mean richness values (Colwell, 2008). We also performed polynomial regressions between species richness and elevation. The best model was selected by choosing the lowest Akaike Information Criterion (AIC).

The area of each elevational value was extracted only in the Hengduan Mountains and summed for each elevational band. Following the methods of previous studies (Vetaas and Grytnes, 2002; Wang et al., 2007), we used species richness/log-transformed area to represent species density for total species and each group.

We measured functional diversity using four common indices: functional richness (Fric), functional dispersion (Fdis), functional evenness (Feve), and functional divergence (Fdiv) (Laliberté and Legendre, 2010; Villéger et al., 2008). We used polynomial regressions to analyze the relationships and trends between functional diversity indices and those of elevational gradients and species richness. The best model was selected by choosing the lowest Akaike Information Criterion (AIC). We also used Spearman correlations to detect relationships among the functional diversity indices and species richness. Functional diversity analysis was carried out in the ‘FD’ package (Villéger et al., 2008), and all the other statistical analyses were performed in R 3.5.3 (R Core Team, 2019).

3. Results 3.1. Species richness patterns, MDE model, and area

Species richness patterns of all groups (i.e., endemic, non-endemic, small, medium, and large range) showed a unimodal distribution across an elevational gradient (Fig. 2a, b). Below 3000 m, the richness pattern of species endemic to the Hengduan Mountains is relatively narrow and asymmetrical compared to that of non-endemic species. At higher elevations, the richness pattern of endemic species shifts markedly compared to non-endemic species. The patterns of elevational richness based on range size are also hump-shaped (Fig. 2b). The richness of species with large and small ranges peaked between 3000 and 4000 m, whereas richness for species with medium ranges peaked between 3500 and 4200 m. The relationship between elevational gradient and different groups of species is polynomial (Table S3). The best polynomial regression model is Richness ∼ Elevation + Elevation² + Elevation³, except best model for small-range species is Richness ∼ Elevation + Elevation². Richness patterns for species with large ranges (adj R² = 0.885, p< 0.01) and non-endemic (adj R² = 0.912, p < 0.01) species were the best fit for the MDE model (Fig. 3). We also found that species density patterns are nearly the same as species richness patterns (Fig. S3).

3.2. Functional diversity and species richness

Functional diversity showed a humped curve along the elevational gradient. Fric and Fdiv showed a peak near 4000 m, whereas Fdis and Feve peaked at about 4500 m (Fig. 4). Spearman correlations showed the coefficient of species richness (SR) against functional richness (Fric) is greater than the coefficients of SR against functional dispersion (Fdis), SR against functional divergence (Fdiv), and SR against functional evenness (Feve) (Table S4). Excluding Feve, functional diversity indices were positively correlated with species richness (Fig. 5).

3.3. Species traits along elevational gradients

Several Delphinium traits decreased at higher elevations, including stem height, leaf width, inflorescence length, sepal spur diam., ratio of the blue or purple color of sepal/petal/staminode, ratio of 3 carpel species, and ratio or branched stem species. Other Delphinium traits increased at higher elevations, including the sepal spur length, other sepal lengths, bracteole length, the ratio of the no-blue-or-purple color of sepal, the ratio of other color types of petals, the ratio of the no-blue-or-purple color of staminodes, the ratio of varied carpel number of species, the ratio of simple and both types of branch type species (Fig. 6). Correlation analysis of all traits showed that some traits are strongly correlated with other traits, e.g., stem height, leaf width, inflorescence length, and flower number (Figs. S4 and S5, and Table S5).

4. Discussion 4.1. Species richness pattern along the elevational gradient

In this study, we found that species richness exhibits a unimodal pattern along an elevational gradient. This finding is similar to previous work (Zhang et al., 2009b) even though we use updated, comprehensive data set. We found that species with a large range fit the mid-domain effect better than species with small or medium ranges (Fig. 2b). The classic signature of geometric constraint is that the peak positions of species with larger ranges are more constrained and overlap more within a bounded domain than do those of species with smaller ranges (Colwell et al., 2004; Currie and Kerr, 2008).

Our finding that endemic species richness has a relatively higher elevational peak that of non-endemic species is consistent with previous studies (Kessler, 2000, 2002; Steinbauer et al., 2016; Wu et al., 2013a, 2013b; Zhang et al., 2009b). Most species with small ranges are endemics (Table S6). Endemic species defined by their planimetric distributions may also have a relatively narrower elevational distribution compared with non-endemic species (t-test, p < 0.001). Endemic species are more adapted to high elevational regions in most mountains (Steinbauer et al., 2016). Higher elevational regions contain more complicated, fragmented, and rugged topography than lower elevational regions (island effect). Complex topography and small habitat differentiation restrict the dispersion of local species, leading to more speciation (Kessler, 2002; Steinbauer et al., 2016).

The potential drivers of elevational patterns of species richness include spatial factors explained by the mid-domain effect and area, but also climatic factors, biotic factors, and evolutionary history (Stiegel et al., 2011). The middle elevational peak of species richness may be caused by the combination of all these factors (O'Brien, 2006; Sun et al., 2020).

4.2. Functional diversity pattern along the elevational gradient

Functional diversity, as measured by Fric and Fdiv, reached a peak near 4000 m (Fig. 4a, c), similar to the peak in species richness. Fric is affected by extreme trait values (Villéger et al., 2008). Accordingly, a high number of species at mid-elevations increases various traits values and covers more space in multi-dimensional functional trait space, which causes Fric to increase with species richness. If unweighted presence-absence data are used, Fdiv reaches the highest value when all points are distributed on multi-dimensional global sphere surface in functional trait space. The Fdiv peak at 3000–4000 m indicates that the points were closer to the mean distance sphere at mid-elevation than at lower or higher elevations. The increasing Fdis value toward the high elevational region reveals the inflation of Fric (Fig. 4b). Fdis calculates the mean distance of species to the centroid of the convex hull (Laliberté and Legendre, 2010). Fric reached a peak at near 4000 m, but Fdis reached the highest values at 4500 m. If there are no extreme trait values at mid-elevations to enlarge the space covered by the convex hull, the Fdis should display the same increase as Fric (Fig. 7). At high elevational bands, lower species richness reduced the species number; thus, fewer species with extreme traits were present at these elevations (Fig. 7). Consequently, the mean distance (Fdis) increased, whereas the space covered (Fric) decreased. To measure the distributional evenness, Feve relates the distance between each point to itself (Villéger et al., 2008). Feve dropped at the highest elevational region, which revealed that species are distributed unevenly in functional trait space at high elevation (Fig. 4, Fig. 7).

Previous research on other taxonomic groups has documented those phylogenetic structures are also clustered at high elevational regions (Li et al., 2014). We also found that patterns of functional traits clustered at high elevations. Species have long been known to show convergent (environmental filtering) and divergent (interspecific competition) evolution (Grime, 2006; Pillar et al., 2009) along elevational gradients, although the main driver of functional diversity patterns is environmental filtering (Li et al., 2014; Zhang et al., 2009b).

Elevation is a proxy for examining how compounded environmental variables (e.g., for temperature and precipitation) shape the patterns of functional diversity (Baur et al., 2014). Under the extreme alpine environmental stress, strong niche overlap (similar traits) among species may transform form convergent functional types into contrasting functional types (Dı́az and Cabido, 2001), e.g., the cushion-like form (Fig. 1d), increasing flower number and flower size. Thus, in the alpine zone, environmental filtering causes unevenness and lower species richness.

4.3. Functional diversity and species richness

Although we found that functional diversity increases with species richness, the rate of this increase in functional diversity decreases and plateaus when the species richness continues to increase. Fric saturation has been previously reported at levels of highest species richness (Chen et al., 2016), and here we observe similar saturation of Fdiv and Feve (Fig. 5). This saturation of functional diversity indices indicates that there is more functional redundancy at the mid-elevational regions. Functional redundancy may ensure that if some species are lost, the remaining species can compensate for the loss; with more species performing similar functions, ecosystems have less chance to fail (Dı́az and Cabido, 2001). Thus, functional redundancy is insurance against functional loss, and the middle elevation regions are the most stable zone for species (Chen et al., 2016; Dı́az and Cabido, 2001; Pillar et al., 2013).

Functional diversity is scale-dependent and scale-sensitive, and may show inconsistencies at various scales (Calba et al., 2014; Epstein et al., 2018; Jarzyna and Jetz, 2018). The applicable scale of the theoretical framework of functional diversity (Villéger et al., 2008) is generally a community scale, and most data used to study functional diversity are obtained from field plots. In recent years, many studies on functional diversity have been conducted at large scales based on large administrative units or large grid units, either along a horizontal or vertical gradient, and have been consistent with research conducted at a fine-scale (Epstein et al., 2018; Shiono et al., 2015; Šímová et al., 2018; Swenson et al., 2012). Our study shows that functional diversity indices are significantly correlated with species richness and that functional redundancy appears along elevational gradients. Although our results suggest that these functional diversity indices are interrelated with species richness, the patterns of functional richness are most similar to patterns of species richness (Table S4). Our results also indicate that although conclusions based on studies at a community scale may extend to larger scales (Calba et al., 2014; Epstein et al., 2018; Jarzyna and Jetz, 2018).

4.4. Species traits along elevational gradients

Our results showed that some interspecific traits increase, whereas others decrease, along an increasing elevation gradient. This is consistent with previous findings that functional traits --- both for a single species (intraspecific) (Wang et al., 2012a, 2012b, 2016, 2018) and for various taxa (Mao et al., 2018; Zhu et al., 2010) vary along elevation gradients. Because plant size could considerably effect the heat balance of plants in alpine regions, shrinking stem size, leaf size, reducing flower number, and decreasing branch ratio may reduce surface area and avoid heat losing (Meng et al., 2013) (Fig. 6a, b, d, o). The increase in carpel number balances the decrease in flower number in terms of reproduction strategy (Fig. 6n).

The increase in flower size may attract more pollinators and increase visiting frequency under extreme alpine environment (Ishii and Harder, 2006; Wang et al., 2012a) (Fig. 6a, i). These trends imply that species may invest more energy in reproduction at higher elevations by shrinking body size and increasing flower size. However, further study is needed to examine reproduction strategy related to seed mass and size. In addition, pollinator differences may potentially affect trait patterns along the elevational gradient. Long spurred Delphinium species are pollinated by some long-beak Nemetrinus species at relatively low elevations (Gao et al., 2020; Hansen et al., 2012), whereas short spurred Delphinium species are pollinated by Bombus species (Macior, 1975; Waddington, 1981) at relatively high elevations. This may explain why spur base diameter and the ratio of other colors of sepal increases and spur length decreases towards higher elevations (Fig. 6g, h, k, l). Delphinium species with pale grey, white or translucent sepal and enclosed flowers (Fig. S2c,d, and f) that are found at relatively high elevations are similar to the “greenhouse” structure that appears in Saussurea and Rheum species. This structure may relate to the need for maintaining temperature inside the flower in cold environments. On the other hand, the acceleration of inflorescence development may ensure its rapid growth and reproduction (Yang et al., 2019).

5. Conclusions

In this study, we mined previously published data to explore regional elevational patterns of functional trait diversity and species richness within the genus Delphinium. We found that total species richness is distributed in a unimodal pattern along an elevational gradient in the Hengduan Mountains. Species endemic to the Hengduan Mountains are distributed relatively narrowly on an elevational gradient and may be affected by high speciation rates and low dispersion rates. Species with large ranges and non-endemic species best fit the mid-domain effect. The mid-domain effect may partly explain the unimodal distribution of total species richness along an elevational gradient, but a combination of additional factors is likely involved as well, including evolutionary history, ecological preference, climatic factors, and biotic factors.

We also found that functional diversity in Delphinium is correlated with species richness and exhibits redundancy in species-rich regions. The most functionally stable region for Delphinium along the Hengduan Mountains elevational gradient is the middle elevational region, likely because plants in these regions are capable of tolerating environmental changes. High elevation alpine regions, which are likely the frontier of speciation, contain more newly evolved endemic species. These narrowly distributed endemic species are trapped by extreme terrain and environmental conditions, and have similar traits, forming a functional convergence group. Although the results of our regional study are consistent with those of community or individual studies, a few caveats are worth mentioning. Data from easy-to-access areas (e.g., near good roads) are likely more reliable and accurate than data from sparsely sampled areas (Chen et al., 2013; Qian et al., 2018); in addition, systematic bias due to sampling may arise in specimen elevation records. Our findings prompt us to recommend that environmental variables related to elevation should be considered in future research and the Hengduan Mountains receive continued study.

Author contributions

H.S. and X.K.O. planned and designed the research; L.S.Q. collected and analyzed the data; L.S.Q. wrote and revised the manuscript. H.H.S. helped to revise the manuscript and figures.

Declaration of competing interest

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

Acknowledgments

We thank Jiahui Chen for his kind help to provide the data source of functional traits. We thank Lu Sun for his kind help with R. This study was supported by the Second Tibetan Plateau Scientific Expedition and Research (STEP) program (2019QZKK0502), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA20050203), and the Key Projects of the Joint Fund of the National Natural Science Foundation of China (U1802232).

Appendix A Supplementary data

The following is the supplementary data to this article: