Statistical analysis
Trait-environment interactions: I used generalized linear mixed-effects models to assess differences between contiguous forest and fragments in the interactions of traits with different environmental variables in explaining species abundances (Brown et al., 2014; Jamil et al., 2013). As detailed in Jamil et al. (2013), the response variable was the abundance of each species per site, including zeros, modeled using a negative binomial error structure. Site-to-site variation in abundance was modeled in relation to a three-way interaction between traits, climate, and habitat type (fragment/contiguous). Plots (nested within sites) and species were included as random intercepts. A zero-inflation component was included to account for the fact that most species were absent in many plots. Zero-inflation was modeled for the intercept and as a function of the trait in each model. Zero-inflated models performed better than models without zero-inflation (based on AIC). I also decomposed the four traits using Principal Components Analysis to get composite phenotypes using the first three axes that explained ~91% of the interspecific variation in trait values (Figure S1). I used ‘varimax’ rotation to maximize the loadings per axis. The first axis loaded for maximum height (0.94) and seed size (0.66), second axis for wood density (0.99), and third axis for SLA (0.96) and partly seed size (-0.55).
I used Principal Components Analysis to decompose the temperature and precipitation variables from WorldClim into two composite axes that together explained ~90% of variation in climate across sites (Figure S2). Climate axis 1 explained 63.2% variation and primarily captured the precipitation gradient, with positive values corresponding with greater precipitation. Climate axis 2 represented a temperature gradient; more positive values indicated warmer sites. Because environmental variables were correlated (Figure S3), I used the Variance Inflation Factor to remove collinear variables (Fox and Weisberg, 2019). Variables with VIF > 4 were removed from the model. Separate models were run per trait. In all models, Climate axis 1, Climate axis 2, CWD, and soil C:N ratio were the final environmental predictors retained. In addition, because elevation was strongly correlated with CWD, Climate axis 1, and Climate axis 2, and had been previously seen as a key predictor of compositional change in this landscape, I implemented models with elevation as the sole predictor for changes in trait-mediated abundances.
RLQ and fourth-corner analyses: The second question was to assess trait-mediated changes in species composition across environmental gradients. This may occur via a) individual traits changing along individual environmental gradients, b) individual traits changing along combination of environmental conditions, or c) trait combinations (functional syndromes) responding to a combination of environmental conditions. To tease apart these possibilities, I used the RLQ and Fourth-Corner methods. RLQ identifies the links between environmental gradients and trait syndromes as reflected by species abundances, and Fourth-Corner analysis performs formal statistical tests of hypotheses regarding trait-environment linkages (Dray et al., 2014; Peres-Neto, Dray, & ter Braak, 2017). For this procedure, R is the matrix of environmental variables by sample sites, Q has species x trait data, and L is a site x species matrix of abundances or occurrences. The RLQ outputs decompose correlations in the data along orthogonal axes for a global measure of trait-environment linkages (Dray et al., 2014). The Fourth-Corner analysis evaluates multiple bivariate correlations between traits and environmental variables, using permutation tests on the RLQ outputs to assess the relationship between (a) individual functional traits and single or combined environmental gradients (RLQ environment scores), and (b) environmental variables and functional syndromes (RLQ traits scores) (Dray et al., 2014).
Following the procedure detailed by Dray et al. (2014), I first analyzed the L -table using correspondence analysis (CA), followed by standardized principal components analysis (PCA) of the R andQ -tables. Then, I used the RLQ method to evaluate the proportions of variation in the environmental ordination and the trait ordination from the RLQ axes. I tested the global significance of trait-environment relationships (RLQ axes) and bivariate correlations between traits and environmental variables with the ”randtest.rlq” and ”fourthcorner.rlq” functions, respectively. To control false positive error rates in multiple comparisons, I used the Benjamini-Hochberg false discovery rate (BH-FDR) procedure (Benjamini & Hochberg, 1995) with 9,999 permutations to estimate P-values.
I ran separate RLQ analyses for fragments and contiguous forest to compare their trait-environment linkages. SLA and seed size were log transformed and all traits scaled prior to analysis. All calculations and graphs were made using package ”ade4” in R version 4.0.1 (R Core Team, 2016). For this analysis, I included all environmental variables together since the multivariate analysis accounted for collinearity.
Trait covariance: To assess how environmental gradients influenced trait covariance, for each plot I calculated the covariance of each pairwise combination of the four traits using methods and R codes in Dwyer and Laughlin (2017). Trait covariance may change in direct response to environmental conditions or because biotic characteristics of sites, such as the diversity of species or functions, vary with environment. As site-level biotic characteristics, I used rarefied species richness, two metrics of multi-trait functional diversity (functional dispersion and functional evenness) and the abundance-weighted variance of each trait in the pair of traits used to assess trait covariance. Rarefied richness provides a measure of species diversity after accounting for differences in stem density. Functional dispersion quantifies the divergence of species in the trait space and incorporates changes in both functional richness and functional divergence (Mason et al., 2013; Mouchet et al., 2010). Functional evenness reveals the uniformity of trait distribution among constituent species. Variance of individual traits provides a measure of the functional space represented by each trait.
The relationships between biotic characteristics the environment can differ between fragments and contiguous forest. So, I first examined how each biotic characteristic varied in relation to an interaction of habitat type (fragment and contiguous forest) and environmental gradients. For this, I used simple linear models. Preliminary analysis with mixed-effects models showed minimal effects of site-level variance on model coefficients and variance explained, and higher AIC than linear models. None of the biotic characteristics varied between fragment and contiguous forest (Table S1) or with the environmental variables considered here (Tables S2 – S4), so I used them along with environmental variables as predictors for trait covariance.
I used linear models to model plot-level covariance of each trait pair as a function of the interactions of habitat type (fragment/contiguous) with environmental variables and biotic characteristics of a plot. As with trait-abundance associations, I used VIF to remove correlated predictors and reduce collinearity. Four climate variables were consistently retained in all models: CWD, Climate axis 1, Climate axis 2, and soil C:N ratio. Because elevation was strongly correlated with the first three of these variables, I repeated all models with elevation as the sole predictor for changes in site-level trait covariance.