Environmental factors related to community-level functional traits in limestone hill forests along an altitudinal gradient: a case study in northern Thailand


 This study investigated the environmental factors affecting functional traits, which have been shown to be important for species assembly in diverse forest stands on limestone hills in northern Thailand. We established 54 plots of 400 m2 in three forest sites (lower, middle, and upper) established along an altitudinal gradient on a limestone hill. The functional traits were assessed and then linked to environmental factors governing forest composition. Results indicated that elevation, rocky outcroppings, and sunlight were important factors affecting functional trait diversity at the study site. Areas with high values of these three factors exhibited increased community-level leaf size, specific leaf area, and leaf thickness, all of which are associated with light-demanding species. However, in areas with low values of these three factors, we observed increased community-level wood density and maximum plant height, which are characteristic of shade-tolerant species. Elevation also positively affected functional dispersion and functional richness values, indicating a wide functional trait space in higher elevation areas, but lower areas exhibited a narrower functional trait space. We suggest that combining a trait-based approach with environmental factors can reveal patterns of species composition in limestone forests.


Introduction
Functional diversity quantifies the diversity of species' traits in biological communities, and this metric is regarded as key to understanding ecosystem processes as well as responses to environmental stress or disturbance (Laliberté & Legendre, 2010). Higher functional diversity signifies great-er differences among species' trait values, more distinct ecological functions, and thus potentially enhanced functional stability against perturbations caused by anthropogenic or environmental stresses (Cardinale et al., 2012). The ability of plants to respond phenotypically to changing environmental conditions provides a reliable tool for the assessment of environment-induced changes in species community composition. Environmental factors can determine which species within the regional pool can thrive in a given community, based on their suitability to local environmental conditions (Cornwell et al., 2009). Previous studies have suggested that directional changes in functional diversity along environmental gradients can be interpreted as evidence of habitat-filtering processes in trait space (e.g., Rolo et al., 2016;Laughlin et al., 2018;Mitchell et al., 2018). However, once species have passed through environmental filters, biotic interactions may also shape the functional structure of communities (Bello et al., 2015). Community traits may converge towards certain values if highly competitive species tend to displace weak competitors (Mayfield & Levine, 2010). On the other hand, functional diversity may increase if minimizing trait similarity is the dominant process. For example, the species regeneration of early successional stages progressed from functionally narrow (i.e., with similar trait values; trait convergence) towards an increasingly wide functional trait space at late successional stages (i.e., trait divergence; Buzzard et al., 2016). Therefore, communities are expected to be shaped by the joint effects of abiotic filtering and biotic interactions.
Limestone hills are elements of karst landscapes, which are characterized by distinctive topography that develops as water dissolves soluble calcareous rock. Tropical forests on limestone hills are viewed as relict communities, with many endemic and rare species (Lynch & Zomlefer, 2016). In Thailand, many studies have demonstrated (e.g., Latinne et al., 2013;Phutthai & Hughes, 2017) that limestone hills provide a critical habitat for endemic and rare species of plants and animals. However, natural disturbances and human activities have complex negative effects on the biodiversity of these limestone hill ecosystems. Therefore, limestone forest restoration is critical, and it is necessary to consider the environmental factors of limestone hill forests when selecting suitable species. Despite such drastic threats to this unique ecosystem, the ecological processes of forest growth on limestone hills in Thailand remain poorly understood (Clark et al., 2014). Therefore, the goal of the present study is to assess the responses of functional traits to environmental conditions that may considerably differ with respect to elevation. Our primary objective was to elucidate key ecological processes to inform necessary improvements in limestone forest management. We focused on determining the environmental factors that most strongly affect tree functional traits in limestone hill forests of Thailand.

Study site
The study area is located on a limestone hill in Phrae Province, northern Thailand ca. 550 km from Bangkok (17°70'-18°84'N, 99°58'-100°32'E) ( Figure 1). The site spans elevations (Elv) of 320-460 m a.s.l. The mean annual temperature and rainfall are 32 °C and 189.6 mm, respectively. The region experiences two main seasons: (i) the wet season (May-October; mean rainfall and temperature of 159 mm and 32.6 °C, respectively) and (ii) the dry season (November-April; mean rainfall and temperature of 23 mm and 25.6 °C, respectively). The dry season is subdivided into cool-dry (November-January) and hot-dry (February-April) sub-seasons (Meteorological Department, 2019).

Field measurements
Field measurements were conducted from October 2015 until December 2016. Three forest sites along an altitudinal gradient on a limestone hill were selected for study thereafter referred as: a lower forest site (LFS), a middle forest site (MFS), and an upper forest site (UFS). The forest sites were differentiated by the Elv and geomorphology along the limestone hill (Table 1). We established a 20 × 20 m permanent plot in each forest site, thereby generating a total of 54 plots for a total area of 2.16 ha (Table 1). We measured the diameters at breast height (DBH at 1.30 m above ground level) for all trees with DBH ≥ 4.5 cm in each established plot. All trees were tagged and identified to the species level. In addition, leaf specimens were collected to verify species identification by comparing them to specimens in the herbarium at the Department of National Park, Wildlife and Plant Conservation. Nomenclature followed the Flora of Thailand (Smitinand, 2014).

Measurements of environmental factors
In addition to elevation we also estimated the percentage of rocky outcropping cover (%RC), percentages of nitrogen (%N), phosphorus (%P), and potassium (%K) in the soil and measured photosynthetically active radiation (PAR), and soil moisture content (SMC) for each of the 54 plots.
To measure SMC and soil nutrients (N, P, K), we collected 100 cm 3 soil samples from the topsoil layer (0-15 cm) in October 2016 using a soil core sampler by taking sub-samples from the center and at each of the four corners of each plot (five points per plot). Two sets of soil samples per plot were collected: the first for analysis of SMC (%), determined as the ratio of fresh weight to dry weight, and the second for soil nutrient analyses (%N, %P, and %K). Measurements of N, P and K were determined using the Kjeldahl, Blay II, and Flame Pho-tometric methods, respectively, and analyses were conducted at the soil laboratory of Chiang Mai University. All samples were collected on the same day, 10 days after the last rainfall and close to the onset of the dry season, when rain was infrequent.
Photosynthetically active radiation (µmol photons m -2 s -1 ) is a measure of the photon flux within the 400-700 nm spectral band of solar radiation. PAR was measured at the center and at each of the four corners of each plot using a Spectrum 3415FX LightScout External Sensor (Spectrum Technologies, Plainfield, IL, USA). PAR measurements were taken at 1.3 m above ground level on a sunny day (08:00-10:00) in November 2016. The average PAR for each plot was calculated as the arithmetic average of five measurements. The average Elv (m) was measured at the center and at each of the four corners of each plot (five points per plot) using a Sun Company Altimeter 202 (Sun Company, Wheat Ridge, CO, USA). The %RC was estimated by counting the number of grid cells wherein rock appeared, which was obtained by dividing each 20 × 20 m plot into 100 2 × 2 m grid cells. These environmental factor measurements were used to represent environmental conditions in our analyses of woody tree abundance and functional traits at each site.

Functional trait measurements
We selected functional traits that were reflective of different ecological strategies for  14 plants and were known to be sensitive to community assembly processes, following Pérez-Harguindeguy et al. (2013). Five functional traits of woody plants were examined and included in analyses: specific leaf area (SLA; cm 2 g -1 ), leaf size (LS; cm 2 ), leaf thickness (LT; mm), wood density (WD; g cm -3 ), and maximum plant height (H max ; m). One individual of each target species was randomly selected per plot to quantify the functional structure of the sampled community following Rolo et al. (2016). Mature sun-leaf samples were taken for assessment of SLA, LS, and LT in October 2016. To determine LS, fresh leaves were scanned using ImageJ software (http:// rsbweb.nih.gov/ij/). Thereafter, LS was measured using the obtained images. Leaf mass was recorded before and after drying at 60 °C for 48 h to a constant weight. SLA (cm 2 g -1 ) was calculated as the ratio of fresh leaf area per oven-dried mass. LT (mm) was estimated as the mean leaf blade thickness taken from five leaf samples, measured using a thickness gauge (China YH-1; Zhejiang, China). Increment cores for WD determination were collected at breast height (1.3 m) from selected individual trees using 5 mm increment borers. WD (g cm -3 ) was calculated from a wood sample (diameter = 0.5 cm, length = 5 cm) as the ratio of the oven-dried mass of the wood sample divided by its fresh volume. The H max (m) of each species was determined using secondary data; obtained from the Flora of Thailand (Smitinand, 2014).

Data analyses
We compared environmental factors among sites (LFS, MFS, and UFS) using an analysis of variance (ANOVA). We calculated the means of measurements made within plots and identified significant differences among sites. For each site, we calculated the following: species richness, stem density (stem ha -1 ), stem basal area (m 2 ha -1 ), species frequency (%, the ratio of the number of plots in which species occurs to the total of all sampling plots), relative stem density (%), relative basal area (%), and relative species frequency (%). From these results, we derived the importance value index (IVI) as the sum of the relative stem densities, relative basal areas, and relative species frequencies for each species in each forest site. We used the IVI to identify the prevalent species in each site. IVI is an excellent indication of the vegetational importance of a species within a site, since it is sensitive to such variables as apparent contagion or exceptional basal area (Curtis & McIntosh, 1951). In addition, the Shannon-Wiener index (H´) was calculated as where pi is the relative abundance of species i, and n is the total number of all species, as a measure of tree species diversity in each site (Hill, 1973). For species-trait relationships, we analyzed all tree species and their functional traits using a principal component analysis ( Community-level means were assessed within each plot to analyze the functional trait diversity of woody trees. We calculated community-level weighted mean (CWM) values for single traits and determined functional dispersion (FDis) and functional richness (FRic) values for multivariate traits (Mouchet et al., 2010). All CWM, FDis, and FRic values were obtained using the 'FD' package (Laliberté et al., 2014) where pi and tr i are the relative abundance and trait value for species i, respectively, and n is the total number of species. The CWM is an estimate of the most probable attribute that a species drawn at random from a community would display (Mouchet et al., 2010). For multivariate traits, we used two indices to quantify each facet of functional diversity for each site, with species distributed in a multidimensional functional space, to describe the variation of functional diversity in each site. We computed the FDis by considering all traits collectively in each site to describe how the species were distributed in trait space (Laliberté & Legendre, 2010). We also computed FRic values by considering all traits collectively in each site to describe the volume of the functional space occupied by the community (Villeger et al., 2008). We then used ANOVA to compare single traits (CWM) and multivariate traits (FDis and FRic) values among sites.
To determine which environmental variables (PAR, %RC, Elv, SMC, %N, %P, and %K) were related to functional trait diversity (single and multivariate traits) of trees in each site, we applied generalized linear mixed models, using the lme4 package (Bates et al., 2015) in R version 3.4.1 software (R Core Team, 2017). The model with the lowest Akaike's Information Criterion (AIC) value was selected as the best-supported model for each species (Ogasawara, 2016).

Environmental factor gradient
Environmental variables (PAR, %RC, Elv, SMC, %N, %P, and %K) differed significantly among the LFS, MFS, and UFS along the elevational gradient (Table 2). Elv significantly differed among the three stands (p < 0.001). The LFS had the highest SMC, followed by the MFS and UFS (p < 0.01). PAR and %RC were higher in the UFS (p < 0.01); however, neither variable differed between the MFS and LFS. For soil nutrients, only K was significantly different among the sites (p < 0.05), with the LFS exhibiting higher K concentration relative to the MFS and UFS (Table 2).

Tree species composition
We identified a total of 21 species of woody trees (Appendix A), from 17 genera and 15 families. The LFS and MFS exhibited higher species richness and diversity than UFS. Both stand basal area (18.9 m 2 ha -1 ) and stem density (2175 stems ha -1 ) were highest in the UFS, followed by the MFS and LFS (

Variation of functional traits along sites
The PCA applied on the functional traits among species explained 63.48% of the total variation. The PCA ordination plot graphically depicted the relationship between tree species and their functional traits across the three forest sites ( Figure  2  The single traits of woody trees differed in the CWM among forest sites (Table 5). The CWM of LS was higher in the MFS and UFS than in the LFS (p < 0.01). SLA was higher in the MFS compared with the LFS and UFS (p < 0.01), and LT was highest in the UFS, followed by the MFS and LFS (p < 0.001). The CWMs of WD and H max were highest in the LFS (p < 0.001 and p < 0.01, respectively), followed by the MFS and UFS. For the multivariate traits, the FRic and FDis values were highest in the MFS and UFS than in the LFS (p < 0.01, and p < 0.001, respectively) ( Table 5). Table 5. Community-level weighted means (CWM) of leaf size (LS), specific leaf area (SLA), leaf thickness (LT), wood density (WD), maximum plant height (H max ), functional richness (FRic), and functional dispersion (FDis) of woody tree species measured in the different forest sites, i.e., lower forest site (LFS), middle forest site (MFS), and upper forest site (UFS), growing along an altitudinal gradient on a limestone hill. Different lowercase letters within a row indicate significantly different means (ANOVA test, p < 0.05; n = 54).

Functional traits of prevalent species
Plant functional traits can drive the growth and establishment potential of species in various forest communities. These traits are often used as proxies to determine whether species utilize different ecological strategies for reproduction and resource capture (McGill et al., 2006). The most abundant species observed in the LFS, Ficus curtipes, had dense wood, suggesting that species in this environment are shade-tolerant and exhibit slow growth rates (Hoeber et al., 2014). The most prevalent species in the UFS, including Excoecari oppositifolia and Dracaena cochinchinensis, tended to exhibit thicker leaves. Species growing in higher, sunnier, drier, and less fertile habitats typically exhibit thicker and longer-lived leaves (Pérez-Harguindeguy et al., 2013), suggesting that these species can balance photosynthetic benefits with the costs of respiration and transpiration and may often contain more photosynthet-ic apparatus per unit area (Onoda et al., 2011). The most widespread species in the MFS, Diospyros castanea and Mallotus peltatus, exhibited a high H max , which suggests that these species are strong competitors for light, as they exhibit short distances between the upper boundaries of the main photosynthetic tissues of the plant and the ground level (Prado-Junior et al., 2017).

Environmental factors affecting functional trait diversity
High functional trait diversity (i.e., high variation in functional traits within a community) is expected to result in enhanced community-level plant productivity and resource-use efficiency (Laureto et al., 2015).
Our results suggested that multivariate functional diversity (FDis and FRis) was highest in the MFS and UFS and lowest in the LFS. In other words, the MFS and UFS exhibited increased functional trait space (i.e., divergent traits), whereas the LFS had a narrower functional trait diversity (i.e., trait convergence) (Laliberté & Legendre, 2010). We found that elevation was a key factor that positively affected multivariate functional diversity within study sites. The higher elevation areas of the limestone hill were characterized by more extensive rocky outcrops, bright light, and dry soil. This unique environment may promote the establishment of specialist species with diverse functional traits. Indeed, limestone outcroppings have been shown to support tropical forest biodiversity at low latitudes around the globe (Felfili et al., 2007). We found that the CWMs of leaf traits (LS, SLA, and LT) were highest in the MFS and UFS; however, the CWMs of WD and H max were highest in the LFS. Facilitative interactions may also occur in woody plant communities of limestone hills at medium and high elevations, as these habitats are associated with faster growing species compared with the species typically found at lower elevation sites (Wright et al., 2010). The increase in %RC with increasing Elv on limestone hills can reinforce canopy gaps and introduce more direct sunlight. These environmental conditions provide a suitable habitat for light-demanding (faster growing) species (Wright et al., 2010). On the other hand, lower elevation areas on limestone hills have few rock outcrop and deeper soils, which are related to higher SMC and nutrient concentrations, with the latter imported from higher elevations. These environmental conditions favor a continuous, tall canopy with low understory light, thus promoting high densities of shade-tolerant (slower growing) tree species (Wright et al., 2010). Our results indicated that three environmental factors, namely Elv, %RC, and PAR, were crucial factors limiting the functional trait diversity within study sites, highlighting that light-demanding (faster growing) species are counterbalanced by shade-tolerant (slower growing) species in tropical limestone forests. Our results are consistent with the concept of environmental filtering of functional diversity. Along the elevational gradient of our limestone hill study sites, we observed a shift from a reduced (convergent) trait space in lower areas toward an increased (divergent) trait space in upper areas, as a result of a higher %RC, Elv, and PAR.

Conclusions
We showed that high Elv enhanced the functional trait space (i.e., divergent traits) in MFS and UFS, while the LFS exhibited narrower functional trait diversity (i.e., trait convergence). Within this limestone hill study area, sites with high Elv, %RC, and sunlight exhibited increased CWMs of LS, SLA, and LT, indicating that these environmental factors were associated with light-demanding (faster growing) tree species. However, in areas with low values of these factors, we observed increased CWMs of WD and H max , indicating associations with shade-tolerant (slower growing) species. Restoration practices within degraded areas on limestone hills must consider the extreme environments created by large areas of bare rock, shallow soil, and high light. To achieve successful restoration in these areas, managers should plant species with large leaves, high SLAs, and/or thicker leaves, as species with these traits are associated with extremely degraded environments.