Leaf trait variation in species-rich tropical Andean forests

Study sites and examined tree species

The study was conducted at three sites in the Andes of southern Ecuador along an elevation gradient at ca. 1000 m (Bombuscaro, Podocarpus NP), ca. 2000 m (San Francisco Reserve) and ca. 3000 m elevation (Cajanuma, Podocarpus NP) in the Provinces of Loja and Zamora-Chinchipe. All sites are located in protected forest areas. At each elevation three permanent 1-ha plots were established in 2018, choosing representative portions of old-growth forest without visible signs of human disturbance (Appendix A1).

The forest types at the three sites differ in floristic composition, species richness and structural characteristics⁴⁹: The premontane rain forest (below 1300 m) at the lowermost site reaches 40 m in height with common tree families being Fabaceae, Moraceae, Myristicaceae, Rubiaceae, and Sapotaceae. It is replaced at 1300–2100 m by smaller-statured lower montane rain forest with Euphorbiaceae, Lauraceae, Melastomataceae, and Rubiaceae as characteristic tree families, and above 2100 m by upper montane rain forest with a canopy height that rarely exceeds 8–10 m. Dominant tree families of the latter forest type are Aquifoliaceae, Clusiaceae, Cunoniaceae, and Melastomataceae. Tree species turnover is complete between premontane and upper montane forest, while a few tree species are shared between lower montane and premontane or upper montane forest types.

The climate is tropical humid with a precipitation peak from June to August and a less humid period from September to December. Mean annual temperature decreases with elevation from 20 °C at 1000 m to 9.5 °C at 3000 m, while annual precipitation increases from around 2000 mm at the two lowermost sites to 4500 mm at 3000 m. Typically, there are no arid months with < 100 mm precipitation^50,51.

The study sites are characterized by relatively nutrient-poor soils on metamorphic schists and sandstones (2000 and 3000 m) or granodioritic rocks (1000 m)⁴³. Soils are slightly more fertile on lower slope positions and at lower elevations than on upper slopes and at upper elevations^43,52. The decreasing nutrient availability is reflected by decreasing forest biomass and productivity with increasing elevation^16,43,53, and along the topographical gradient from lower to upper slope position⁸.

We selected 52 tree species in total, 20 tree species each from the species-rich premontane and lower montane forests, and twelve species from the upper montane forest. The species included the most abundant tree species in the permanent 1-ha plots at each elevation and in addition represented the major occurring tree life strategies from each study site, covering fast-growing pioneers to late-successional tree species and understory species to tall canopy trees (Appendix A2). In addition, we used available data on SLA and wood specific gravity (WSG) from previous studies⁵⁴ to select species at the three sites that covered the known range of these two functional traits. The selected species represent 41.0% (at 1000 m), 50.3% (at 2000 m) and 32.0% (at 3000 m) of the total tree basal area (trees ≥ 10 cm dbh) of the three study plots at the respective elevation.

Leaf sampling and trait analyses

We sampled 8–10 replicate trees per species, randomly selected from the known individuals within the three 1-ha plots (for four species only 5–7 appropriate tree individuals were found). In total 421 trees were sampled.

The sampling campaign took place from February to March 2019. We collected 2–3 branches per tree from the top of the crown with as much sun exposure as possible with all attached leaves and took them in sealed polyethylene bags (filled with water-soaked tissues) to the research station, where further processing took place. The branches were placed over night in water to achieve water saturation of the leaf tissue prior to the measurements.

Per tree, 20 young but fully developed sun leaves without signs of herbivory were stripped from the branches and used for the analyses. We took care that all leaves were fully expanded and not senesced and attempted to cover the full range of leaf sizes present. Three additional leaves of average size were taken to quantify leaf thickness and toughness. In case of compound leaves, we measured the morphology of the individual leaflets. In a few species with particularly large leaves (e.g. Cecropia, Graffenrieda, and Pourouma), a smaller number of leaves was investigated due to time constraints during optical leaf area determination.

For determining average leaf size (one-sided leaf area, LA, cm²), the fresh leaves including petioles were scanned in color (Canon LIDE 100, 150 dpi). LA was determined from the scanned leaf silhouettes with the software WinFOLIA 2014a (Régent Instruments, Quebec, QC, Canada). Subsequently, leaf fresh weight was determined and the leaves dried at 60 °C for at least three days to determine leaf dry weight and foliar water content. Leaf dry matter content (DMC, mg g⁻¹) was calculated as the quotient of leaf dry weight to leaf fresh weight. Specific leaf area (SLA, cm² g⁻¹) was calculated by dividing total LA by total leaf dry mass.

Leaf thickness (mm) was measured on three fresh leaves per tree with a digital micrometer (Mitutoyo M293-240-70, Mitutoyo Germany Ltd, Neuss, Germany) at each two locations on both sides of the main leaf vein in the middle between the secondary veins; measurements were subsequently averaged.

Leaf toughness (kN m⁻¹) was estimated as the mean of six punch tests using a digital penetrometer (flat-ended 2.0 mm diameter steel punch, DS-50 N, Imada Inc., Japan) on three fresh leaves (excluding the midrib and other major veins) from each tree⁵⁵.

Leaf dry mass was analyzed for its C and N content with a CN elemental analyzer (Vario EL III, Hanau, Germany), and for its P, Ca, K, Mg and Al content by ICP analysis (Thermo Scientific iCAP 7000 ICP-OES, Thermo Fisher Scientific, Germany) after HNO₃ digestion of the ground leaf material.

Statistical analyses

We compared trait means across the three elevation levels with Tukey’s HSD test.

The phylogenetic relationships of the studied tree species were extracted from the mega-tree of vascular plants ‘GBOTB.extended.tre’ using the R package V.PhyloMaker⁵⁶ (Appendix A3). We used the R package brms⁵⁷ to fit 12 Bayesian phylogenetic multilevel models⁵⁸ with identical structure to describe the effect of the elevational level on the 12 investigated leaf functional traits. All leaf traits were log-transformed in order to handle skewness and heteroscedasticity. The observational units were the tree individuals. The models had the following structure:

where the response variable (Y), a log-transformed leaf trait, is described by a normal distribution with the varying mean (mu) and a standard deviation of (sigma). The predicted value (mu) for observation (i) in plot (j) and of species (k) is described by the intercept (alpha_{0}), a fixed effect for the categorical predictor site with three levels, a random intercept for plot (j) ((alpha_{plotleft[ j right]})) and a random intercept for species (k). Whereas the plot effects were standard random effects as in usual mixed models, the random intercepts for the species were a combination of two components: a phylogenetic effect that incorporated the phylogenetic non-independence of residuals (called ‘phylogeny effect’ in short) and an phylogenetically-independent species effect that accounted for additional variance among species (called ‘species effect’ in short). The random effects were fit as following:

Random plot effects

Random species and phylogeny effects

For the random plot effect, a standard deviation parameter (tau_{plot}) was estimated as in usual mixed models. The two-component random effect on the species level was estimated using two variance parameters, (tau_{phyl}^{2}) for the ‘phylogeny effect’ and (tau_{ind}^{2}) for the ‘species effect’. The overall effect described by a multivariate normal distribution with a covariance of (Sigma_{phyl}). (Sigma_{{phylleft[ {m,n} right]}}) is a matrix with (tau_{phyl}^{2} + tau_{ind}^{2}) on the diagonal and (tau_{phyl}^{2}) times the correlation (rho_{{phylleft[ {m,n} right]}}) derived from the phylogenetic variance–covariance matrix between the two species (m) and (n) elsewhere.

This model structure can be interpreted analogously to a classical PGLS model but offers the advantages of handling several observations per species instead of working with species means, and the possibility to include further random effects and comprehensive Bayesian inference⁵⁸.

Weakly informative normal priors were used for the slope and intercept parameters. All variance components were assigned to weakly informative half-t priors. Models were fit with Hamiltonian Monte Carlo (HMC⁵⁹) via the Stan probabilistic programming language⁶⁰ using R package brms⁴⁹. Sampling was performed for 10,000 iterations after a warmup of 10,000 iterations. The settings used for the HMC algorithm were an adapt_delta value (target acceptance rate) of 0.99 and a maximum tree depth of 15. The posterior distributions of the parameter estimates were summarized by the posterior mean and the 95% highest posterior density intervals (HDI⁶¹). Effects were considered credibly different from zero, when the 95% HDI did not include zero. The contribution of different model components to the total variance in the data was decomposed based on the approach of⁶², extended to a multi-level context analogous to⁶³.

We followed the method proposed by de Bello et al.⁶⁴ for partitioning of quadratic entropy with the Rao index to decompose total community variance into between-species and within-species effects. The R function ‘RaoRel.r’⁶⁴ was used to compute the variance weighting by relative species abundance. Species relative abundance was based on the species’ average contribution to plot basal area in the three permanent 1-ha plots at the respective elevation level. For the analyses, we standardized each trait by dividing the values by the range of possible values for this trait⁶⁴.

A principal components analysis (PCA) was used to summarize the correlation structure of the 12 standardized mean trait values for the 52 tree species at the three study sites.

We applied network analyses to assess the connectivity between the studied traits at each elevation level. For each study site, we calculated a matrix of trait-trait relationships using Pearson correlations. To avoid considering spurious correlations among traits, only statistically significant correlations (p < 0.05) were included^32,33. Then, an adjacency matrix was created by assigning below threshold as 0 and above threshold as the integer of |r|× 10; thus, the adjacency matrix only shows the weighted presence or the absence of connections between pairs of plant traits. Finally, trait networks were visualized and informative parameters calculated with the R package igraph⁶⁵.

In the constructed networks, traits are represented as nodes and their correlations are represented as the edges linking them. We used “degree”, “node strength”, and “betweenness” as indicators of network centrality. The degree is the number of edges that connects a focal node to other nodes; node strength is the sum of the weights of edges linking a focal node to adjacent nodes, and betweenness gives the number of shortest paths from all nodes to all others passing through the focal node.

To characterize the overall networks, we used “edge density”, “average path length” and “average clustering coefficient”. Edge density defines the proportion of present edges among nodes out of all possible edges in the network; average path length is the average “degrees of separation” between all pairs of nodes in the network, and the average clustering coefficient measures the probability that the adjacent nodes of a node are connected.

All statistical analyses were performed using the R v4.0.2 programming environment (https://www.r-project.org/)⁶⁶.

Source: Ecology - nature.com