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

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

Figure 1

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

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³ 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 Photometric 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⁻² s⁻¹) 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.

We selected functional traits that were reflective of different ecological strategies for 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² g⁻¹), leaf size (LS; cm²), leaf thickness (LT; mm), wood density (WD; g cm⁻³), 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² g⁻¹) 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⁻³) 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).

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⁻¹), stem basal area (m² ha⁻¹), 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 (1)H′=∑i=1npi × lnpiH' = \sum\nolimits_{i = 1}^n {pi\, \times \,{ln}pi} 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 (PCA) performed in PC-ORD software (Version 6, McCune & Mefford, 2011). PCA is a statistical procedure that uses an orthogonal transformation to convert a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables (McCune & Mefford, 2011). The complete data matrix used in this study was 21 tree species and 5 variables for functional trait (SLA, LS, LT, WD, and H_max). Eigenvalues and eigenvectors were extracted from the covariance matrix of original variables to determine the relationship between the variables.

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) in R version 3.4.1 (R Core Team, 2017). For single traits, we calculated a CWM to describe the dominant trait value in each site. CWM was calculated as (2)CWM=∑i=1npi ×tri{\rm{CWM}} = \sum\limits_{i = 1}^n {p_i \, \times tr_i } 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 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).

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² ha⁻¹) and stem density (2175 stems ha⁻¹) were highest in the UFS, followed by the MFS and LFS (Table 3). The five most important value indices of tree species within each forest site are shown in Table 4. Streblus taxoides, Ficus macleilandii, and Dracaena cochinchinensis were the most prevalent species at all elevations. Ficus macleilandii and Dracaena cochinchinensis were prevalent in the MFS and UFS, while Alphonsea siamensis and Sterculia pexa were pervasive in the MFS and LFS. The most common species in the MFS were intermediate between those in the LFS and UFS. Parishia insignis prevailed in all three forest sites. Streblus taxoides and Aglaia edulis were extensive in the LFS only, whereas Sterculia lanceolata and Meliosma pinnata were widespread in the UFS only (Table 4).

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). Several dominant species displayed distinguishing characteristics, such as a high H_max in Mallotus peltatus, Firmiana colorata, and Diospyros castanea, thicker and larger leaves in Dracaena cochinchinensis, Streblus taxoides, and Sterculia villosa, and a high SLA and dense wood in Ficus macleilandii, Ficus curtipes, and Brucea javanica.

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

For the multivariate traits, the FDis and FRic values were positively related to Elv (p < 0.01) (Table 6). For single traits, the CWMs of all functional traits were negatively and/or positively related to the %RC and Elv (Table 6). PAR was positively related to the CWMs of all three leaf traits (LS, SLA, and LT; Table 6). The CWM of LS was significantly related to all environmental factors, exhibiting positive relationships with PAR, %RC, Elv, and %P, and negative relationships with SMC, %N, and %K. The CWM of SLA was positively affected by PAR, %RC, Elv, and %P and negatively affected by %N. The CWM of LT was positively, significantly related to PAR, %RC, and Elv. The CWM of H_max exhibited a positive relationship with SMC but a negative relationship with %RC and Elv (Table 6).

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

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.