Woody Plant Structural Diversity Changes across an Inverse Elevation-Dependent Warming Gradient in a Subtropical Mountain Forest

: Examining the changes in woody plant structural diversity along an inverse elevation-dependent warming gradient will enhance our mechanistic understanding of how warming affects forest communities because such an inverse elevational gradient reflects a warming trend in a mountain landscape. Here, we investigated the effects of warming on the patterns of species composition and structural diversity in a subtropical broadleaved forest. We calculated a warming index based on elevational difference and modeled the aspect-related potential incident radiation (PDIR) using nonparametric multiplicative regression. We tested the changes in structural diversity of three communities for significant differences along the warming gradient. We associated both the warming index and PDIR with the principal components and tested their relationships for significant differences. We found that trees of different sizes varied in their response to the warming gradient. While a significant decreasing trend was exhibited in both species diversity and size diversity for trees of all sizes and for adult trees along the warming gradient, no significant changes in seedlings were found, and the average basal area value was the highest for the warmest community. Our findings demonstrated that a short-range elevational gradient was adequate to separate the communities in species composition and structural diversity. Patterns of structural diversity along the warming gradient varied in size classes. The community at a higher elevation had more indicator species that were unique in separating the community from others. Principal component analysis showed that the first two principal components were negatively correlated with the warming index, indicating that warming destabilized species composition and community structure. Our study suggests that warming is the major driver of changes in structural diversity and species composition of woody plant communities in a subtropical broadleaved forest and that warming may promote tree productivity at the community level but reduce structural diversity at the quadrat level


Introduction
Ecologists and forest resource managers have long been interested in elevational patterns of forest communities and their structural characteristics [1][2][3][4] because these patterns represent the responses of forest communities to microclimate change [5][6][7].In a mountain landscape, temperature decreases as elevation increases at a scale of 0.65 • C/100 m according to the common environmental lapse rate [8,9].Conversely, an inverse elevational gradient corresponds to a warming gradient and can be used to predict warming effects on the structural attributes of forest communities with the "space-for-time substitution" method [10].
Past studies have reported contrasting elevational patterns of structural diversity in forest communities.Some studies revealed that species diversity was negatively correlated Forests 2024, 15, 1051 2 of 15 with elevation, while other studies indicated that species diversity was higher with increasing elevation [11][12][13].The changes in community structure and species diversity along elevational gradients have long attracted the interest of ecologists and forest scientists.The response of forest community to climate warming can be shown by examining the changes in community structure and species diversity across an inverse elevational gradient [5,10], because in a mountain landscape a warming trend exists with descending elevation according to an environmental temperature lapse rate of 0.65 • C per 100 m [8,9].The subtropical mountain forest ecosystem is an ideal natural experimental environment for studying the variation of plant community structure and diversity along the elevational gradient, because the subtropical forest is rich in species resources [14,15], and the mountains with sloping topography can directly or indirectly affect the distribution of thermal radiation, the key factor controlling the richness, abundance, and dominance of plant species in a community [10,16].Therefore, understanding the changes in composition, structure, and diversity of subtropical forest communities at different elevations can provide valuable information for the forestry sector and the scientific community to cope with climate change.
Here, we focused on the association of a warming trend along an inverse elevational gradient with the changes in the structural diversity of woody plant communities.We calculated a warming index based on the common temperature lapse rate and modeled an energy-related variable, the potential direct incident radiation (PDIR), using nonparametric multiplicative regression (NPMR).We examine the changes in structural diversity metrics across communities along an inverse elevational gradient.We also related the changes in species composition and structure of the communities to the warming index and PDIR.
We aimed to answer the following questions: (1) Does structural diversity differ along an inverse elevation-dependent warming gradient?(2) How patterns of species richness, total abundance, species diversity, size diversity, and species composition are affected by the inverse elevation-dependent warming? and (3) Can the communities be separated by certain structural parameters along a short-range elevational gradient?

Study Sites
Our study sites were located in the Yinpingshan Nature Reserve (22 • 52 ′ -22 • 56 ′ N, 114 • 05 ′ -114 • 15 ′ E), along a slope on Mount Yinping, the core zone of the nature reserve.This area is in the southeastern part of Dongguan, South China's Guangdong Province (Figure 1), and has a lower subtropical humid climate [17,18].The mean annual temperature ranges from 21 • C to 22 • C, with the maximum mean temperature of 27.5 • C in July and the minimum mean temperature of 13.5 • C in January [17].The mean annual precipitation is 1800 m, approximately 80% of which is concentrated from April to September [17].Soil type is clay loamy latosol developed from granite.The peak of Mount Yinping is at an elevation of 898.2 m, the highest peak in the local area of Dongguan.Vegetation in the area is dominated by subtropical evergreen broadleaved forests [17,18].

Sampling Design and Data Collection
We established three 1 ha square plots on a forested slope of Mount Yinping in the nature reserve after topographic surveying with a total station (Nikon DTM-310, Nikon Geotecs Co., Ltd., Tokyo, Japan).We placed the plots in a relatively homogeneous evergreen broadleaved forest and separated them by at least 100 m in a vertical direction.Each plot was divided into 100 square quadrats or 10 × 10 m grid cells for plant and topographic data collection.The quadrats were numbered and marked by placing PVC (polyvinyl chloride) tubes at the 10 m corners.Elevation at each 10 m corner was recorded for calculation of slope and aspect of each quadrat [19].Elevation for each quadrat was calculated by averaging the elevation values at its four corners [20].

Sampling Design and Data Collection
We established three 1 ha square plots on a forested slope of Mount Yinping in the nature reserve after topographic surveying with a total station (Nikon DTM-310, Nikon Geotecs Co., Ltd., Tokyo, Japan).We placed the plots in a relatively homogeneous evergreen broadleaved forest and separated them by at least 100 m in a vertical direction.Each plot was divided into 100 square quadrats or 10 × 10 m grid cells for plant and topographic data collection.The quadrats were numbered and marked by placing PVC (polyvinyl chloride) tubes at the 10 m corners.Elevation at each 10 m corner was recorded for calculation of slope and aspect of each quadrat [19].Elevation for each quadrat was calculated by averaging the elevation values at its four corners [20].
We tallied all trees ≥ 1 cm diameter at breast height (DBH) in each quadrat.All censused individuals were identified to species and recorded with species name and DBH value measured to the nearest 0.1 cm.For exact DBH measurement, trees ≤ 5 cm DBH were measured with a dial caliper, while those > 5 cm DBH were measured with a diameter tape.Plant taxonomy and systematics followed Ye and Peng [21].All the field data were entered in a spreadsheet to construct species-by-quadrat datasets for further calculations.

Statistical Analysis
We calculated species richness (number of species), total abundance (number of stems), the Shannon-Wiener diversity index, and Pielou's evenness for each quadrat using the species-by-quadrat dataset.Shannon-Wiener diversity index and Pielou's evenness were calculated with the following formula [22]: We tallied all trees ≥ 1 cm diameter at breast height (DBH) in each quadrat.All censused individuals were identified to species and recorded with species name and DBH value measured to the nearest 0.1 cm.For exact DBH measurement, trees ≤ 5 cm DBH were measured with a dial caliper, while those > 5 cm DBH were measured with a diameter tape.Plant taxonomy and systematics followed Ye and Peng [21].All the field data were entered in a spreadsheet to construct species-by-quadrat datasets for further calculations.

Statistical Analysis
We calculated species richness (number of species), total abundance (number of stems), the Shannon-Wiener diversity index, and Pielou's evenness for each quadrat using the species-by-quadrat dataset.Shannon-Wiener diversity index and Pielou's evenness were calculated with the following formula [22]: where H ′ is the Shannon-Wiener diversity index, P i is the relative abundance of the i-th species in proportion to the total abundance of species in a quadrat; and E is the Pielou's evenness, S is the number of species in a quadrat.Size diversity was represented by variability in basal area (BA), size diversity index, and evenness [15].We defined the size diversity as the variability in basal areas across quadrats.The size diversity index and evenness were calculated following Su, Jia, Zhang and Chen [15], which are based on the Shannon-Wiener index formula and Pielou's evenness, respectively.
where D size is the size diversity index, P BAi is the relative basal area of the i-th species in proportion to the total basal area of species in a quadrat.
To test the effect of the inverse elevation-dependent warming on structural diversity of woody plant communities, we first grouped all tree individuals into three size classes, i.e., seedlings, saplings, and adults, as defined by DBH values using the following criteria [16]: seedlings, 1 cm ≤ DBH < 2.5 cm; saplings: 2.5 cm ≤ DBH < 12.5 cm; adults: DBH ≥ 12.5 cm.Therefore, we calculated the species diversity, size diversity, and evenness indices by four size groups, i.e., all sizes, seedlings, saplings, and adults, for each quadrat.We then used Kruskal-Wallis test, a nonparametric alternative to analysis of variance (ANOVA) that is appropriate for analyzing field ecology data [16], to test for significant differences in these diversity metrics across the three communities along an inverse elevational gradient.
We calculate an inverse elevation-dependent warming index for each quadrat using the following formula [10]: where Elev max is the elevation of the highest quadrat, while Elev quadrat is the elevation of a particular quadrat including the highest one.Therefore, the elevational difference of the highest quadrat, or hence the warming index, equals zero [10].We performed principal component analysis (PCA) on the species-by-quadrat abundance dataset to reduce the dimensions to fewer orthogonal axes.PCA is a conventional ordination method that allows the analysis of datasets of non-independent variables, such as species composition data [24].We extracted the scores or eigenvectors of the first two PCA axes and assessed the effects of warming on tree species composition by relating the axis scores to warming index and PDIR by axis overlays.
To evaluate the differences in species composition across communities of an inverse elevation gradient, we performed multi-response permutation procedures (MRPP).MRPP is a nonparametric multivariate method for testing the hypothesis of no difference between two or more groups of entities [25], such as plant communities and forest stands with compositional data.The output of MRPP includes an agreement statistic A, which describes within-group homogeneity, and the test statistic T, which describes the separation between the groups [25].A p-value is calculated based on randomizations in MRPP [25].
To assess the contributions of some key species to the resemblance as well as separation of groups (communities), we plot the quadrat occurrence, species abundance, and total basal area of the shared species with all presence in the three communities.We also performed indicator species analysis (ISA) to identify species indicative of a community of a particular elevation.We calculated indicator values for each species in each community using Dufrêne and Legendre's method [25], which combines information on the concentration of species abundance in a particular community and the faithfulness of occurrence of a species in that community [25].An indicator value for a species in a particular community was expressed as a percentage.The highest indicator value of a species across communities was regarded as the overall indicator value of that species for a particular community [18,25].The cut-off of the indicator value as a criterion for detecting indicator species for a particular group is usually arbitrary.An indicator value of 25 may, on average, correspond to a relative abundance of 50% and a relative frequency of 50%, because the indicator value is based on the product of relative abundance and relative frequency [25].
The calculation of diversity metrics, PCA ordination and axis overlays, MRPP, and ISA were carried out with PC-ORD: Multivariate analysis of ecological data, version 7 (MjM Software, Gleneden Beach, OR, USA).Descriptive statistics and Kruskal-Wallis tests were executed with Statistica version 8 (Statsoft, Inc.Tulsa, OK, USA).

Community Structure
The three forest communities we investigated were relatively homogeneous in the way that they had similar dominant species (Table 1).The total number of taxa summed up from all the communities, or the γ-diversity, was 152 species from 105 genera and 54 families, representing 11,149 individual trees (Table 1).The total basal area in the aggregate for each community showed a general decreasing trend with decreasing elevation, which was consistent with the change in total abundance (Table 1).For the number of taxa, the community at the highest elevation is greater than the community at the lowest elevation; however, the biggest number of species, genera, and families were found in the community at the mid-elevation (Table 1).Twenty-eight species were simultaneously shared by the three communities at different elevations (Figure 2).These species with all presence in the communities represented 7247 individual trees, approximately 2/3 of the total number of individuals found in the communities; however, they differed in quadrat occurrence, species abundance, and total basal area across communities (Figure 2).As shown in Figure 2a, the five most common species with an overall quadrat occurrence were Schima superba Gardner et Champ.(248), Itea chinensis Hook.et Arn.(205), Machilus chinensis (Benth.)Hemsl.(158), Schefflera octophylla (Lour.)Harms (98), and Psychotria rubra (Lour.)Poir.(89).Species varied in their distribution among the communities.Both Schima superba and Itea chinensis had approximately 90% quadrat frequency in Community 1 as well as in Community 3; however, Schima superba had 66% quadrat frequency (66 out of the 100 quadrats) in Community 2, whereas Itea chinensis had only 24% quadrat frequency in Community 2.
As shown in Figure 2b, the five most abundant species summed up from all three communities were Schima superba (2182 stems), Itea chinensis Hook.et Arn.(1765 stems), Machilus chinensis (756 stems), Litsea rotundifolia Hemsl.(314 stems), and Psychotria rubra (286 stems).Seventeen of the shared species had an abundance of less than 100 stems.The most abundant species, Schima superba, had 841 stems distributed in Community 1, 609 in Community 2, and 732 in Community 3, while other species were distributed in a more aggregated pattern.For example, the second most abundant species, Itea chinensis, had 758 stems distributed in Community 1 and 918 stems in Community 3, but only 89 stems in Community 2. Eight of the less abundant species had only one stem distributed in one of the communities and several to dozens of stems in two other communities (Figure 2b).The distribution of basal area of the species with all presence (Figure 2c) shared similar patterns with that of abundance (Figure 2b), as well as with quadrat occurrence among the communities, with only a few exceptions, such as for the two most abundant species, Schima superba and Itea chinensis.As shown in Figure 2b, the five most abundant species summed up from all three communities were Schima superba (2182 stems), Itea chinensis Hook.et Arn.(1765 stems), Machilus chinensis (756 stems), Litsea rotundifolia Hemsl.(314 stems), and Psychotria rubra (286 stems).Seventeen of the shared species had an abundance of less than 100 stems.The most abundant species, Schima superba, had 841 stems distributed in Community 1, 609 in Community 2, and 732 in Community 3, while other species were distributed in a more aggregated pattern.For example, the second most abundant species, Itea chinensis, had 758 stems distributed in Community 1 and 918 stems in Community 3, but only 89 stems in Community 2. Eight of the less abundant species had only one stem distributed in one of the communities and several to dozens of stems in two other communities (Figure 2b).The distribution of basal area of the species with all presence (Figure 2c) shared similar patterns with that of abundance (Figure 2b), as well as with quadrat occurrence among the communities, with only a few exceptions, such as for the two most abundant species, Schima superba and Itea chinensis.

Species Diversity
Trees of different sizes responded differently to the inverse elevation-dependent warming gradient (Figure 3).For all sizes of trees (Figure 3a,e), as well as for adult trees (Figure 3d,h), quadrat-level species richness and total abundance significantly decreased across communities along an inverse elevation-dependent warming gradient; but for seedlings, no significant difference was found across a warming gradient (Figure 3b,f).The trends of change in both species richness and total abundance were consistent for trees of different sizes.For saplings, total abundance decreased with warming (Figure 3c), though this trend was not significant, while species richness significantly decreased with warming (Figure 3g).Within-community variability in both total abundance (Figure 3c) and species richness (Figure 3g) was the highest for saplings in Community 1 and lowest in Community 3, which was at the warmest elevation; but for adult trees, the warmest community (Community 3) exhibited broader min-max ranges and had more outliers and extreme values (Figure 3d,h).
trees of different sizes.For saplings, total abundance decreased with warming (Figure 3c), though this trend was not significant, while species richness significantly decreased with warming (Figure 3g).Within-community variability in both total abundance (Figure 3c) and species richness (Figure 3g) was the highest for saplings in Community 1 and lowest in Community 3, which was at the warmest elevation; but for adult trees, the warmest community (Community 3) exhibited broader min-max ranges and had more outliers and extreme values (Figure 3d,h).
The Shannon-Wiener diversity index and Pielou's evenness (Figure 4) showed similar trends as compared with total abundance and species richness (Figure 3), except that the Shannon-Wiener index of seedlings exhibited an increasing trend with warming (Figure 4b), though not significant.The Shannon-Wiener diversity index and Pielou's evenness (Figure 4) showed similar trends as compared with total abundance and species richness (Figure 3), except that the Shannon-Wiener index of seedlings exhibited an increasing trend with warming (Figure 4b), though not significant.

Size Diversity
Changes in size diversity across communities along a warming gradient (Figure 5) were different from those of species diversity (Figure 4).Changes in basal area across communities along a warming gradient showed a significant decreasing trend for trees of all sizes (Figure 5a, p < 0.00001) and adult trees (Figure 5d, p < 0.00001).A weak decreasing trend was found for saplings with marginal significance (Figure 5c, p = 0.065).Size diversity significantly decreased along a warming gradient for trees of all sizes (Figure 5e, p < 0.00001), saplings (Figure 5g, p < 0.0001), and adult trees (Figure 5h, p < 0.0001).No significant difference was found in the basal area (Figure 5b, p = 0.313) and size diversity index (Figure 5f, p = 0.176) for seedlings across communities.Changes in the basal area (Figure 5a-d) and size diversity index (Figure 5e-h) were similar to those in total abundance and species richness (Figure 3), respectively, indicating that the basal area was determined more by total abundance than by species richness and size diversity was determined more by species richness than by total abundance.

Size Diversity
Changes in size diversity across communities along a warming gradient (Figure 5) were different from those of species diversity (Figure 4).Changes in basal area across communities along a warming gradient showed a significant decreasing trend for trees of all sizes (Figure 5a, p < 0.00001) and adult trees (Figure 5d, p < 0.00001).A weak decreasing trend was found for saplings with marginal significance (Figure 5c, p = 0.065).Size diversity significantly decreased along a warming gradient for trees of all sizes (Figure 5e, p < 0.00001), saplings (Figure 5g, p < 0.0001), and adult trees (Figure 5h, p < 0.0001).No significant difference was found in the basal area (Figure 5b, p = 0.313) and size diversity index (Figure 5f, p = 0.176) for seedlings across communities.Changes in the basal area (Figure 5a-d) and size diversity index (Figure 5e-h) were similar to those in total abundance and species richness (Figure 3), respectively, indicating that the basal area was determined more by total abundance than by species richness and size diversity was determined more by species richness than by total abundance.No significant difference was found in the size diversity-based evenness for seedlings, saplings, and trees of all sizes across communities, but for adult trees, the evenness significantly decreased across communities along a warming gradient (Figure 5m, p < 0.00001).

Species Composition
Multi-response permutation procedures (MRPP) indicated that the species composition of the communities at different elevations extremely differed by overall comparison as well as by pairwise comparison (Table 2).The species composition of Community 1 vs. Community 2 (MRPP, A = 0.10) was more homogeneous than either Community 1 vs. Community 3 (MRPP, A = 0.08) or Community 2 vs. Community 3 (MRPP, A = 0.09).The between-group difference was greater for Community 1 vs. Community 2, rather than for Community 1 vs. Community 3, or for Community 2 vs. Community 3.

Species Composition
Multi-response permutation procedures (MRPP) indicated that the species composition of the communities at different elevations extremely differed by overall comparison as well as by pairwise comparison (Table 2).The species composition of Community 1 vs. Community 2 (MRPP, A = 0.10) was more homogeneous than either Community 1 vs. Community 3 (MRPP, A = 0.08) or Community 2 vs. Community 3 (MRPP, A = 0.09).The between-group difference was greater for Community 1 vs. Community 2, rather than for Community 1 vs. Community 3, or for Community 2 vs. Community 3. Principal component analysis (PCA) well separated the three communities by the ordination of sample units in species space (Figures 6 and 7).Two warming-related variables, warming index and potential direct incident radiation (PDIR), were significantly correlated with the first two principal components from PCA.As shown in Figures 6 and 7, the warming index was negatively correlated with PCA axis 1 (Figure 6; r = −0.659,p = 0.0001) and Axis 2 (Figure 6; r = −0.359,p = 0.001), while PDIR was positively correlated with PCA axis 1 (Figure 7; r = 0.390, p = 0.001) and negatively correlated with PCA axis 2 (Figure 7; r = −0.338,p = 0.001).
Table 2. Multi-response permutation procedures (MRPP) for the species composition across communities.Explanation of the MRPP statistics: A is a measure of effect size and represents the within-group homogeneity; T is a statistic describing the separation between the groups; p is the p-value for significance test of homogeneity using randomization.Community 1, Community 2, and Community 3 are each represented by 1 ha plot at 576.9 m, 433.2 m, and 239.5 m a.s.l., respectively.Principal component analysis (PCA) well separated the three communities by the ordination of sample units in species space (Figures 6 and 7).Two warming-related variables, warming index and potential direct incident radiation (PDIR), were significantly correlated with the first two principal components from PCA.As shown in Figures 6 and 7, the warming index was negatively correlated with PCA axis 1 (Figure 6; r = −0.659,p = 0.0001) and Axis 2 (Figure 6; r = −0.359,p = 0.001), while PDIR was positively correlated with PCA axis 1 (Figure 7; r = 0.390, p = 0.001) and negatively correlated with PCA axis 2 (Figure 7; r = −0.338,p = 0.001).The inverse elevation-dependent warming index was stronger than the aspect-related PDIR in their relationships with the first two PCA ordination axes (Figures 6 and 7).The warming index also had greater strength in separating the communities than PDIR (Figures 6 and 7).

Indicator Species
Twenty-two tree species were detected as bio-indicators of a particular community with a significant indicator value ≥ 25 (Table 3).These indicator species contributed to the separation of the communities.Nine out of these indicator species were found with all presence in the three communities (Figure 2 and Table 3).These were Acronychia pedunculata (L.) Miq., Adina pilulifera (Lam.)Franch.ex Drake, Ardisia quinquegona Blume, Evodia meliaefolia Benth., Itea chinensis Hook.et Arn., Machilus chinensis (Benth.)Hemsl., Psychotria rubra (Lour.)Poir., Sapium discolor (Champ.)Muell.-Arg., and Schima superba Gardner et Champ.These species only had low to medium indicator values (Table 3).The number of indicator species that separated a community from others with a significant indicator value was the highest for Community 2 (10 species), lowest for Community 3 (3 species), and medium for Community 1 (9 species), practically a decreasing trend along a warming gradient.The species with the highest indicator value (IV = 77.0)for indicating Community 1, Rhodoleia championii, had a low indicator value averaged across the three communities by randomization (Table 3).By contrast, Schima superba had an indicator value of 35.8 for indicating Community 1, but it had the highest mean indicator value (IV = 30.6)for all the communities from randomization (Table 3).This result suggested that Schima superba was not uniquely indicative of a particular community.Although Schima superba was the most abundant species, its relative abundance and relative frequency in Community 1 as calculated among the three communities was not as high as other species with less abundance but in aggregated distribution, such as Adina pilulifera, Machilus chinensis, and Psychotria rubra.The inverse elevation-dependent warming index was stronger than the aspect-related PDIR in their relationships with the first two PCA ordination axes (Figures 6 and 7).The warming index also had greater strength in separating the communities than PDIR (Figures 6 and 7).

Indicator Species
Twenty-two tree species were detected as bio-indicators of a particular community with a significant indicator value ≥ 25 (Table 3).These indicator species contributed to the separation of the communities.Nine out of these indicator species were found with all presence in the three communities (Figure 2 and Table 3).These were Acronychia pedunculata (L.) Miq., Adina pilulifera (Lam.)Franch.ex Drake, Ardisia quinquegona Blume, Evodia meliaefolia Benth., Itea chinensis Hook.et Arn., Machilus chinensis (Benth.)Hemsl., Psychotria rubra (Lour.)Poir., Sapium discolor (Champ.)Muell.-Arg., and Schima superba Gardner et Champ.These species only had low to medium indicator values (Table 3).The number of indicator species that separated a community from others with a significant indicator value was the highest for Community 2 (10 species), lowest for Community 3 (3 species), and medium for Community 1 (9 species), practically a decreasing trend along a warming gradient.
The species with the highest indicator value (IV = 77.0)for indicating Community 1, Rhodoleia championii, had a low indicator value averaged across the three communities by randomization (Table 3).By contrast, Schima superba had an indicator value of 35.8 for indicating Community 1, but it had the highest mean indicator value (IV = 30.6)for all the communities from randomization (Table 3).This result suggested that Schima superba was not uniquely indicative of a particular community.Although Schima superba was the most abundant species, its relative abundance and relative frequency in Community 1 as calculated among the three communities was not as high as other species with less abundance

Discussion
One of the most important spatial patterns that attracted ecologists and forest managers is the elevational or altitudinal pattern of plant communities [26][27][28][29].In a broad-range elevational gradient, structural diversity is controlled predominantly by the coupling of water and heat [30,31], thus the changes in diversity patterns are similar to those in latitudinal patterns.In contrast, in a short-to-medium-range elevational gradient, changes in diversity patterns will be affected predominantly by temperature because of similar water availability from rainfalls.A small change in temperature due to the temperature lapse rate in mountain areas will exert an enormous influence on structural diversity [15,31,32].
Structural diversity is mainly composed of three components, i.e., spatial distribution, species diversity, and tree size diversity [15,33].Changes in structural diversity are often affected by a variety of factors, including climate change [34,35], land use history [36], habitat heterogeneity [37], and inter-and intraspecific competitions [38,39].Compared with other indicators, tree size is easy to measure using forestry tools [15,40].Several biodiversity indices have been used in previous studies to evaluate the response of diversity to elevational gradients, however, the richness-or abundance-based diversity indices may not be so sensitive in response to warming as tree size-based diversity indices [16].For species diversity, the most common indices are the Shannon-Wiener index, Simpson index, and Pielou's evenness [40,41].The Shannon-Wiener index reflects the species richness and evenness in a community, which is estimated by measuring the uncertainty of the category to which the randomly selected individuals belong, while the Simpson index represents the likelihood that two randomly chosen individuals in a community will be different species, and Pielou's evenness index can reflect the uniformity of the number of individuals in the community [42].However, traditional diversity indices cannot fully reflect the structural diversity of a community, because the diversity of a forest stand means a wide range of changes in tree species, size, and age class [16,40].Therefore, an extension of the well-known Shannon-Wiener diversity index with DBH or basal area to diversify the diversity structure is more reasonable.This new index based on size structure not only has a similar diversity calculation method to the Shannon-Wiener index but also does not need to divide the diameter and basal area data into different classes.As shown in this study, using the size data to calculate the Shannon-Wiener index will undoubtedly provide a new indicator for the relationship between diversity and elevational gradient.
Previous reports have shown that patterns of species diversity along an elevational gradient varied with studies.Species diversity may decrease or increase [11][12][13], or it will exhibit a unimodal pattern or just unrelated [43] with elevation.These patterns might be site-specific so more studies should be carried out for adequate information to reach a sound conclusion, especially when these studies are taken into account the relationships with climate change.A review report by Chinese scholars documented various patterns of species diversity along an elevational gradient, with the prevailing conclusion that species diversity decreases with the increase in elevation [43].Obviously, the decrease in species diversity is reasonable with a broad-scale elevation increase.Vegetation changes with a broad-scale elevation increase, such that the species diversity of scrublands, coniferous forests, or mountaintop dwarf forests that often occur at high elevations is lower than that of evergreen broadleaved forests at low elevations.In contrast, our results are completely different, that is, species diversity increases with elevation, or in other words, species diversity decreases along the inverse elevational warming gradient.This might be due to the fact that our study specifically examined the change in structural diversity of broadleaved forest communities along a short-range elevational change.It also indicated that only a short-range elevation-dependent warming gradient can make appreciable changes in the structural diversity of the subtropical broadleaved forests.
Although we set up a short-range inverse elevational warming gradient using three evergreen broadleaved forest plots at different elevations and tested the effect of warming on forest biodiversity, a broad spatial gradient, as well as a temporal gradient in combination, should be considered for study in the subtropical area.A broad spatial gradient combined with varied forest types may show differential patterns of changes in structural diversity, which calls for comparative studies.Moreover, whether the global warming trend and extreme climatic events, such as massive floods and cyclones that frequently occur in this area, will alter the response of forest structural diversity to warming should also be explored.

Conclusions
Our results demonstrated that a short-range elevational gradient was adequate to separate the communities in a homogeneous broadleaved evergreen forest in species composition and structural diversity.Patterns of species richness, total abundance, species diversity, and tree size diversity along an inverse elevational gradient varied in size classes and were scale-dependent.The first two principal components were negatively correlated with the warming index, indicating that warming destabilized species composition and community structure.The community at a higher elevation had more indicator species that were unique in separating the community from others.Our study suggests that warming is the major driver of changes in structural diversity and species composition of woody plant communities in a subtropical broadleaved forest.

ForestsForests 15 Figure 1 .
Figure 1.Location of the study sites in Dongguan, Guangdong Province, China.

Figure 1 .
Figure 1.Location of the study sites in Dongguan, Guangdong Province, China.

Figure 3 .
Figure 3. Boxplots showing changes in total abundance (a-d) and species richness (e-h) across communities at different elevations.The community codes represent the average elevations of a community: 1, 576.9 m; 2, 433.2 m; and 3, 239.5 m.

Figure 3 .
Figure 3. Boxplots showing changes in total abundance (a-d) and species richness (e-h) across communities at different elevations.The community codes represent the average elevations of a community: 1, 576.9 m; 2, 433.2 m; and 3, 239.5 m.

Figure 4 .
Figure 4. Boxplots showing changes in the Shannon-Wiener diversity index (a-d) and Pielou's evenness (e-h) across communities at different elevations.The community codes represent the average elevations of a community: 1, 576.9 m; 2, 433.2 m; and 3, 239.5 m.

Figure 4 .
Figure 4. Boxplots showing changes in the Shannon-Wiener diversity index (a-d) and Pielou's evenness (e-h) across communities at different elevations.The community codes represent the average elevations of a community: 1, 576.9 m; 2, 433.2 m; and 3, 239.5 m.

Figure 6 .
Figure 6.PCA ordination diagram showing changes in sample units with overlays of warming index on the first two axes.

Figure 6 .
Figure 6.PCA ordination diagram showing changes in sample units with overlays of warming index on the first two axes.

Figure 7 .
Figure 7. PCA ordination diagram showing changes in sample units with overlays of potential direct incident radiation (PDIR) on the first two axes.

Figure 7 .
Figure 7. PCA ordination diagram showing changes in sample units with overlays of potential direct incident radiation (PDIR) on the first two axes.

Table 1 .
Summary of community characteristics for the three 1 ha plots at different elevations.

Table 2 .
Multi-response permutation procedures (MRPP) for the species composition across communities.Explanation of the MRPP statistics: A is a measure of effect size and represents the withingroup homogeneity; T is a statistic describing the separation between the groups; p is the p-value for significance test of homogeneity using randomization.Community 1, Community 2, and Community 3 are each represented by 1 ha plot at 576.9 m, 433.2 m, and 239.5 m a.s.l., respectively.

Table 3 .
Woody plant species indicative of a particular community with a significant indicator value ≥ 25.