2.1 Study site and experimental design

Experimental plots were established on an area that was previously used as a pasture near La Gamba in the southern region of Costa Rica (8°42′03.78″N, 83°12′06.14″W). The area is located in the humid tropical lowlands with a mean annual temperature of 28.1°C and mean annual precipitation of 5,870 mm (average 1997–2018) measured at the Tropical Research Station La Gamba. Rainfall is <400 per month from December to April, but even the driest month (February) receives 160 mm on average. Soils are loamy clays (Inceptisol—Aeric Endoaquept) and can be flooded in the wet season.

Seeds and seedlings were collected from the local forest and grown for three to six months in a nursery. Between 2012 and 2014 trees were planted in 56 experimental plots of 6 × 6 trees at a spacing of 3.5 m × 4 m with a buffer strip of at least one tree row between plots. Trees were classified as low WD (≤0.5 g/cm³) species, high WD (>0.5 g/cm³) species and legumes, with WD obtained from published data (Zanne et al., 2009) and own local measurements. Trees of these functional groups were planted in combinations of one, two or three functional groups always with nine species per plot, each with four individuals, which were distributed randomly within the plots. Overall, the species planted were selected depending on the availability of seeds or seedlings and species planted in plots were those for which WD data were available, or that were legumes. There was no target to plant a defined number of individuals per species as not all were available in high numbers. The plots were established as a long-term study on the effect of variable functional diversity at constant species diversity on processes and ecosystem functions. We do not expect such effects to become evident before canopy closure and several years of plant growth and here therefore evaluate individual tree growth, but control for potential plot effects in statistical models. During the first three years after planting the understory vegetation approximately one meter around the trees was cleared manually.

2.2 Measurements of growth and traits

Tree height (the distance of the uppermost leaves to the flat ground) was measured in 2013, 2015 and 2016 using a marked pole for trees up to 4 m or a hypsometer (Vertex IV; Haglöf) for taller trees. Diameter was measured at 1.3 m when trees had reached this height. The dataset obtained, included a total of 5,307 height measurements for 2,217 individuals out of 152 different species. For species-specific analyses we here evaluate only 47 species with a minimum of 40 measurements from 12 to 65 (mean: 30) individuals per species.

Leaf and wood traits were obtained from sun-exposed branches of 3–14 (mean 4.3) individuals of the planted trees per species. WD was determined for branches approximately 1-cm thick and 10-cm long. We removed bark and pith, measured volume by the water displacement method and calculated WD (in g/cm³) as dry weight per fresh volume. Wood samples were dried at 80°C for at least 3 days. Theoretical hydraulic conductivity (Kh) was estimated from xylem vessel size and densities. For this, branch wood sections (10–30 µm thick) were cut with a microtome, stained with astral blue and safranin before embedding in Euparal to produce permanent microscopic slides. Digital images were taken with a DM5500B light microscope (Leica, Germany) and all vessels (n = 23–235, mean: 98) within a radial section extending from pith to bark were colour-coded with Photoshop CS6 (Adobe Systems) and subsequently measured automatically with ImageJ 1.49o (http://imagej.nih.gov/ij).

Vessel size and density (VD per mm²) were used to calculate theoretical hydraulic conductivity as

where ρw is the density of water at 20°C (998.2 kg/m³), η is the viscosity of water at 20°C (1.002 × 10^–3 Pa s) and Dh is the hydraulically weighted vessel diameter (m, Sterck, Zweifel, Sass-Klaassen, & Chowdhury, 2008).

Specific leaf area was obtained by scanning the lamina (excluding the petiole and the rachis) with a desktop scanner and dividing the area by the lamina dry weight. LDMC was calculated by dividing the leaf dry mass by its water-saturated fresh-mass. Leaves were re-hydrated by placing them in a humid and dark container with the petioles submerged in water for 12–24 hr (Garnier, Shipley, Roumet, & Laurent, 2001). Leaf nitrogen concentration (N) was measured in 1–2 mg finely ground leaf powder with an elemental analyzer (EA 1110; Thermo Fisher) coupled to an isotope ratio mass spectrometer (Delta^PLUS; Finnigan MAT) by a ConFlo III interface (Thermo Fisher) at the SILVER laboratory of the University of Vienna.

2.3 Habitat specialization

Habitat specialization was determined based on a multinomial classification method using a conservative supermajority threshold (K = 2/3) and a high significance level (p = 0.005; see Chazdon et al., 2011). In addition to data on species abundance compiled by Letcher et al. (2015), we included local data from twelve 500 m² plots in secondary forests (Oberleitner, 2016) and from 15 OG and five secondary forest plots of 1 ha (Morera-Beita et al., 2019). Species were classified as OG specialists, SG specialists, generalists and those which were too rare to classify (Figure S1). As an alternative to this categorical classification of habitat specialization (HS_cat), we calculated a continuous variable of habitat specialization (HS_con), which is the number of individuals recorded in OG habitats divided by the total number of individuals for each species and ranges between 0 (only recorded in SG) and 1 (only from OG forests). This is similar to the successional index (Peña-Claros, 2003), but this index includes the age of forests, which is often unknown and thus available for fewer species.

2.4 Data analysis

To compare tree growth limited by little or no competition by neighbours, we measured growth during the first years after planting and before canopy closure would shade smaller individuals. Because some trees died and others were planted later, the mean observation period was 662 (±276) days. Since young trees particularly of slow-growing species were <1.3 m tall for the first measurements, diameter at breast height (DBH) data were not available for these trees, which would bias results against slow-growing trees and species if DBH-growth data were used. We thus evaluated height growth, which is commonly used in similar settings (Gustafsson et al., 2016; Martínez-Garza, Bongers, & Poorter, 2013; van Breugel et al., 2011). Height and diameter growth were strongly correlated (Figure S2). Because height (or any other measure of) growth is often not constant over time or plant size, comparing different groups of individuals can be problematic. We therefore fitted species-wise functions (linear, three-parameter logistic and Gompertz) to the height: age (the time of planting being age = 0) relationship (Paine et al., 2012, 2015) and used the best-fitting function (Figure S3), selected by the lowest AIC (Akaike information criterion), to calculate size-standardized height increment for an individual of 1 m height. We choose 1 m because some species had barely grown >1 m and using a greater height would have meant extrapolating the fitted growth functions beyond the data range in these species. Mean growth rates per species over the entire size range were mostly very similar to size-standardized growth (Table S1).

We used linear models (LM) to test for correlations among functional traits (WD, SLA, LDMC, leaf N and Kh), growth and habitat specialization (HS_con and HS_cat) using species-mean values. As trees were growing in individual plots and age might have an effect on growth rates, we also calculated linear mixed effect (LME) models with species, plot and the median age between two height measurements of individual trees as random variables. Size would not be suitable as a random variable because size at any given age will be strongly related to growth rate. To simplify comparisons, particularly of the variation explained, with other studies, which mostly use species-mean values to analyse the relationship between growth and traits, we primarily present LMs, and show LME results in Supporting Information.

To test if correlations were caused by relatedness between species, we calculated regressions with phylogenetically independent contrasts (PICs; Felsenstein, 1985) using the phylogeny GBOTB.extended in the r package V.PhyloMaker, which is based on an angiosperm mega-phylogen (Zanne et al., 2014). PICs were calculated with r package picante with regressions forced through the origin.

Height growth rates and Kh were log transformed to homogenize variances, and HS_con, which ranges between 0 and 1, was arcsine transformed (Quinn & Keough, 2002). For LMEs, trait data were scaled by dividing by the root mean square because variables strongly differing in magnitude can bias LMEs. LMEs were calculated with the function lmer of r package lmerTest, which includes significances for terms using the Satterthwaite approximations to estimate degrees of freedom. The r package MuMIn was used to calculate pseudo r² of LMEs as the variance explained by fixed factors (Nakagawa & Schielzeth, 2013). For the relationship between HS_con and traits LMEs were not appropriate because HS_con is not related to plots or individual trees, so only LMs were used here. We calculated principal component analyses (PCA) to evaluate associations among traits and between traits, growth and habitat specialization. One PCA was calculated with the leaf and wood traits (WD, Kh, SLA, LDMC and N), a second PCA included standardized growth. We tested if the first principal component, i.e. the major axis of the trait-space, with or without growth, was correlated with HS_con. Since the classification of habitat specialization is based on widely differing numbers of records, which will affect the accuracy of the classification, all regressions with HS_con were weighted by the logarithm of the number of records per species.

We explored the causal relationship between functional traits, growth and habitat specialization with structural equation models (SEM). For the initial model standardized growth was explained by wood and leaf traits and HS_con by traits plus growth. The model was reduced stepwise by removing non-significant parameters, which also reduced AIC. Including relationships among wood or leaf traits (such as WD ~ Kp) never improved model fit for growth or HS_con and always increased AIC, and was therefore omitted from the models. SEMs were calculated with R library piecewiseSEM. All analyses were performed in R version 3.5.1 (R Core Team, 2015).

3.1 Differences among successional habitat preferences

Based on the available records, most of the species could be assigned to one of the three habitat specialization types: OG specialists (10 species), generalists (16) and SG specialists (11). Five species were too rare (n < 15) in the available records to be classified and five species had no habitat records at all (Table S1).

Specific leaf area was positively correlated with leaf N (r = 0.451, p = 0.003), and negatively with LMDC (r = −0.377, p = 0.012). WD was negatively correlated with Kh (r = −0.407, p < 0.011) and correlations between SLA and WD as well as Kh were marginally significant (p < 0.1; Table 1).

Of all wood and leaf traits, only LMDC differed significantly among the three classes of habitat specialization (HS_cat; Table 2). If only OG and SG were compared, OG specialists had significantly lower SLA and higher LDMC values. Kh was significantly and LMDC marginally significantly (p < 0.1) correlated with the continuous measure of habitat specialization (HS_con). HS_cat and HS_con explained similar proportions in the variation of traits, but overall the r² of ANOVAS using HS_cat was somewhat higher than r² of linear regressions with HS_con. Note that for the linear correlations also species with very few records in OG or SG forests were included (and weighted by the log of the number of records), whereas five species with too few records were not classified for the ANOVA with HS_cat. Because the species were weighted by the number of records in linear regressions, omitting the five species with few records had hardly any effect on the trait: HS_con relationships (not shown).

Standardized height growth ranged from 0.25 m per year in Minquartia guianensis to 4.96 m in Schizolobium parahyba (Table S1) and was nearly twice as high in SG specialists (mean 2.00 ± 0.31 SE) than in OG specialists (1.18 ± 0.14). Pairwise comparisons of growth rates showed significant differences between OG and SG specialists (p = 0.01) and marginally significant differences (p = 0.095) between SG specialists and generalists, but none between OG specialists and generalists. Growth rates were significantly related to habitat specialization, tested as HS_cat or HS_con, which explained 21% and 15%, respectively, of variation in standardized growth rates (Figure 1). Omitting five species that had too few records to reliably classify as HS_cat slightly improved the correlation between HS_con and growth (r² = 0.19, p = 0.0067).

3.2 What explains growth and habitat?

Wood traits tested were significantly correlated with standardized height growth (r² = 0.183, p = 0.005 for WD and r² = 0.156, p = 0.012 for Kh), while SLA was marginally related (p = 0.075) and LDMC and N were not related to growth (Table 3). PIC correlations were very similar to species-wise correlations (Table S2), so these correlations were not mainly caused by the relatedness between species. Statistically, habitat specialization (Table 2) explained growth better than any functional trait except WD, though HS per se is not a trait that can functionally explain growth. When a potential plot effect was accounted for by treating plot as an error term in mixed effect models, again WD, Kh and HS_con were significantly correlated with tree growth (Table S3).

In the PCA with traits including growth (Figure 2) WD scales opposite to standardized growth with leaf traits roughly orthogonally to the WD–growth axis. HS_con is significantly (p = 0.007, r² = 0.25) correlated with PC1, which explains 40.3% of the trait variation. If growth is not included in the PCA (Figure S4), the correlation between HS_con and PC1 is lower, but still significant (p = 0.012, r² = 0.21). While there is an overlap, SG specialists, generalists and OG specialists also separate mainly along axis 1 of the PCA (Figure 2).

The SEM that best explained the relationship between traits, standardized growth and HS was the model where growth was explained by WD, and HS by growth and N (Figure 3). The final model had a lower AIC than more complex models and explained 23% of the variation in growth and 42% of the variation in HS_con. A model where HS_con was explained only by growth and not by leaf N had still lower AIC, but was rejected as it explained substantially less of HS_con (22%) and the effect of N was highly significant in the selected model. The starting full model is shown in Figure S5 and model fit criteria and parameters are shown in Table S3.

4.1 Habitat and functional traits

Based on data of species from a Costa Rican lowland rainforest we found habitat specialization being surprisingly poorly correlated with most functional traits that we investigated. Binning species into three groups (SG and OG specialists and generalists) yielded significant differences only for LDMC, whereas a classification using the proportion of species recorded in OG forests, which is a continuous variable, found significant or marginally significant effects for Kh and LDMC (Table 2). The likelihood of finding significant effects certainly depends on the number of species included and the way species are classified. However, using the same classification method as we did and substantially more species (66 SG and 133 OG specialists; generalists were omitted here), a study of Costa Rican canopy trees found significant differences only in WD, but none (p > 0.5!) in SLA and LDMC (Lasky et al., 2014), although interestingly trait diversity in SLA and LDMC increased with stand age. Also, a large compilation of functional traits and successional habitat for neotropical forests comparing pairs of OG and SG forests in various locations found no significant differences in five out of five pairs for seed mass, in four out of five pairs for seed length, in three out of eleven pairs for SLA, in one out of five pairs for LDMC and in four out of ten pairs for WD (Letcher et al., 2015). Thus, while the trends in our results do agree with the pattern reported by others (significant or not, WD and LDMC were higher in OG, and leaf N, SLA and Kh higher in SG specialists), in many cases the relationship between functional traits and the successional status of tree species in tropical forests is weak.

The lack of strong relationships is somewhat surprising, since specialization along a successional gradient requires specific functional capacities to establish, grow and survive under the given environmental conditions. Light-demanding species are typically characterized by functional traits that permit rapid growth and fast adjustment to changing light conditions, such as high SLA, high leaf N, low WD and low LDMC (Bazzaz & Pickett, 1980; Chazdon, 2014; Kitajima, 1996; Valladares & Niinemets, 2008). While many SG species are more light demanding and OG species shade-tolerant (Strauss-Debenedetti, 1996), shade tolerance and habitat selection should not be treated as equivalent. In a study in subtropical China shade tolerance, quantified via the abundance of saplings in low-light microsites, explained only 26% of the successional class, assessed via the abundance in forests of different ages, and other measures of shade tolerance explained even less (Feng et al., 2018). Environmental filtering may act more strongly to exclude light-demanding pioneer species from establishing in shaded understorey than to exclude shade-tolerant OG specialists from establishing early in succession (Poorter & Arets, 2003).

4.2 Growth and functional traits

To compare growth rates among tree species and to relate species-wise growth with traits, growth has been quantified by various measures including diameter growth, basal area growth, RGR or height growth. One consideration in the choice of the measure is that the quantification should ideally be independent of plant age or size, but over a longer time-span or size gradient, this is true for no available measure of growth. For instance RGR in woody plants tends to decline strongly with size (Iida et al., 2014), as biomass is increasingly allocated to non-productive tissues. As species have different growth trajectories (Groenendijk, Bongers, & Zuidema, 2017) a measure of growth that is constant for one species may thus be increasing or decreasing with age in others. To account for this, we used size-standardized growth (Paine et al., 2012) or accounted for a possible age effect in LMEs. Still, if the growth functions are nonlinear the standardized growth rates depend on the standard size chosen (1 m in our case because some species had not grown much taller), so using a different size might still affect the outcome.

A second issue that may confound the comparison of growth rates is that most studies in natural forests compare individuals growing under widely different light conditions. Given that light-demanding species tend to be found in open canopy microsites whereas in shade-tolerant species at least the juveniles grow under low light, this may produce a biased comparison of the species-specific growth rates. This bias can be partially overcome by selecting the top quantile of growth rates and assuming that this represent individuals at least not growing under strong light limitation (Wright et al., 2010), by comparing individuals growing in a similar light environment (Poorter & Arets, 2003), or by accounting for light in statistical models (Rüger et al., 2012). The second and third approach require a reliable quantification of the light environment for individual trees inside forests, which is challenging to obtain. This problem is avoided when all species are grown under uniformly high light, but such settings are uncommon. Reforestation projects typically operate with a limited selection of trees and also prefer species that grow fast in open fields, sometimes planting shade-tolerant species when the first generation of trees provides shelter (Lamb, 2011). The reforestation at Finca Amable started with a large number of species across the successional gradient, which provided an ideal setting to test the intrinsic growth of species related to their traits and habitat specialization.

Wood density was negatively and Kh was positively associated with potential growth rates, while the effect of leaf traits was much lower. The dominant effect of WD on growth is consistent with previous studies (Poorter et al., 2008; Wright et al., 2010), which mostly also found lower correlations with leaf traits.

Another difference between our study and many others is that we evaluate juvenile growth, thus differences between our and other studies might be explained by architectural differences and ontogenetic shifts in carbon allocation. A meta-analysis (Gibert et al., 2016) confirmed model predictions that the effect of SLA on growth decreases with plant size whereas the effect of WD remains strong during ontogeny, though that study was largely based on seedlings. Changes in light availability are related to functional trait plasticity, which tends to be particularly high for leaf traits because of their close linkage to the photosynthetic capacity and carbon balance of the plant (Rozendaal et al., 2006). In contrast, WD also changes among ontogenetic stages (Hietz, Valencia, & Wright, 2013) but does not respond to changes in irradiance to the same extent as leaf traits, and therefore displays less intraspecific variation (Siefert et al., 2015). If important functional traits associated with growth display life stage specific variation, this may result in additional variation in the trait-growth relationships for juvenile and adult trees.

4.3 Growth and habitat affinity

Standardized height growth was the best predictor of habitat specialization (HS_con, Table 2) and juveniles of SG specialists grew nearly twice as fast as OG specialists under conditions of little competition. This confirms previous studies that observed significant differences in growth rates related to habitat specialization, which was seen in seedlings under high and low light (Poorter, 2005) and in one study also in saplings under high light (Martínez-Garza et al., 2013). Most other studies compared growth of SG and OG specialists in natural habitats, where OG specialists would grow under restricted light conditions and differences in growth rate could be partially explained by differences in resource (light) availability. A large reforestation study in Panama using many species did not define successional status (van Breugel et al., 2011), but comparing height growth from that study with HS_con from the dataset we used gave very similar results to our study with HS_con explaining 26% of height growth (Figure S6).

Our study was conducted in a very humid climate (annual rainfall of c. 6,000 mm). Whereas early successional stages in tropical wet forests are dominated by species with more acquisitive traits including low WD and high growth rates, in tropical dry forests early successional species tend to be conservative with high wood densities and drought avoidance strategies (Lohbeck et al., 2015; Poorter et al., 2019). We are not aware of any study testing potential growth rates of species along a tropical dry forest succession, but we do not expect that the relationships found in our and other studies of wet forests would also be seen in strongly water limited tropical forests.

Looking for causal relationships among traits (including growth) and habitat affinity showed WD to be the factor that most explains growth and leaf N the only other trait related to HS_con. Thus the frequently reported link between WD and successional habitat in tropical wet forests (Chazdon, 2014; Poorter et al., 2019) appears to be indirect via its effect on potential growth (Figure 3). With the exception of leaf N, pairwise correlations between leaf or wood traits and HS_con (Table 2) are mainly indirect, if they are significant at all. Ecological specialization as well as growth of species must be explained by some traits, but even the best model explained only part of the variation of HS_con or growth. This means that either the quantification of habitat specialization is not ideal or that, apart from growth and leaf N, the traits commonly associated with the successional gradient are not the relevant ones.

4.4 Scaling the successional gradient

The ecological niche is an important concept to explain species distribution and diversity (Johnson, Condit, Hubbell, & Comita, 2017). In forests it is partly defined by habitat specialization along a successional gradient, which is often used to classify species in tropical forests. Shade tolerance is frequently used as a proxy for successional habitat specialization. It is mostly defined as the ability to grow in shade, and can be quantified by the frequency of juveniles found in different light environments in the forest (Poorter, Bongers, Sterck, & Wöli, 2005), demographic rates (the growth–mortality trade-off; Wright et al., 2010) or responses to light gradients such as the light-compensation point (Valladares & Niinemets, 2008). Habitat specialization, by contrast, is more directly assessed by the occurrence of species in SG versus OG forests (Chazdon et al., 2011) or in forests of different ages, which was found to be the best indirect measure of shade tolerance (Feng et al., 2018). Since the age of many secondary forests is difficult to obtain, we used the classification system of Chazdon et al. (2011), to classify trees into SG and OG specialists and generalists. Comparing the predictive power of either these three classes (HS_cat) or a continuous gradient (HS_con) shows that a simple measure based on the relative number of records from OG forests is about as good in explaining growth or functional traits (Table 2; Figure 1). Clearly, the classification of species with few records is not very reliable, but this can be easily dealt with in statistical models and including these (small black symbols in Figure 1) does not much affect the outcome of our analysis. Including species with few records, available data (Letcher et al., 2015) thus at least tentatively classify more than 1,500 neotropical species. This dataset can easily be amended with local data and similar data can be compiled for other regions. Neither the categories nor the continuous classification used here accounts for the number of regions or plots the species were recorded from, which might introduce some bias or over-confidence in the classification of species recorded in high numbers from a single site. Compiling a similar dataset with information on forest age would provide a still more accurate classification of habitat specialization. Given the only moderate correlation with measures of shade tolerance (Feng et al., 2018), we caution against using shade tolerance and successional habitat specialization as equivalent and studies should specify clearly what they are testing. Comparing standardized and quantitative measures of habitat specialization, shade tolerance, growth and traits would help to better understand the variation in ecological strategies and niche specialization of the hyperdiverse tree communities in tropical forests.

The successional niche is defined by various environmental factors including resource availability (light, nutrients, water), abiotic stress and competition. In addition, reproductive traits can differ between SG and OG species (Chazdon, Careaga, Webb, & Vargas, 2003), dispersal can shape successional forests (González-Castro, Yang, & Carlo, 2019) and seed removal was found to decrease with stand age (Peña-Claros & De Boo, 2002), factors that are not reflected in the traits we studied. While the traits most commonly studied are easy to measure and available for many species, they are imperfectly correlated with performance traits (including growth) that may be more relevant. While correlations between habitat and leaf or wood functional traits have provided many insights, our study showed that these can be indirect and thus may invite to misleading interpretations. We suggest that analysing traits that are more directly functional, including growth rates, shade tolerance or drought resistance in the case of the successional gradient, might be more informative and advance our understanding of what ecological specialization is based on, although such data are often more laborious to obtain.