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.