Statistical analysis
To examine if tree seedling density (number of seedlings per m2) was higher overall on nurse logs than on the forest floor, we used a Wilcoxon rank test because the data did not meet parametric model assumptions of normality and homogeneity of variance across groups. We used seedling density per area of the plot because area of the plot varied with nurse log diameter. We then examined whether substrate type including forest floor and all three decay classes of nurse logs influenced tree seedling density, bryophyte depth, percent canopy cover, percent cover of the three most abundant bryophyte species (Rhizomnium glabrescens , Hylocomium splendens ,Antitrichia curtipendula ) using 1-way ANOVA or Kruskal-Wallis tests if parametric model assumptions couldn’t be met. Bryophyte depth was missing from one plot and percent cover of H. splendens had to be square-root transformed to meet parametric model assumptions.
To determine the factors associated with seedling density, we used a generalized linear model (GLM) with a Poisson error distribution. Explanatory variables included substrate type (nurse log or forest floor), decay class of nurse log, bryophyte depth, % canopy cover, and bryophyte species with a total percent cover of >300 across all plots and that were found in at least 10 plots (Table S1). The bryophyte species included as predictors were: Hylocomium splendens , Rhizomnium glabrescens , Antitrichia curtipendula , Rhytidiadelphus loreus , Sphagnum girgensohnii , and Kindbergia praelonga . We used the dredge function in the MuMIn package to determine the model with the lowest Akaike Information Criterion (AIC) value (Barton, 2009). We calculated Spearman rank correlation coefficients for all pairs of explanatory variables using the rcorr function in the Hmisc package in R (Harrell & with contributions from Charles Dupont and many others, 2020). We included forest floor as decay class 5 in our model. We also calculated variance inflation factors (vif) to test for collinearity among our predictor variables using the vif function in the car package; any vif < 5 was deemed not correlated.
To examine the influence of substrate (forest floor and nurse logs characterized by decay class) on bryophyte community composition, we conducted a non-metric multi-dimensional scaling ordination (NMDS) using the metaMDS function in the vegan package in R (Oksanen et al., 2010). We used a Bray-Curtis distance matrix on bryophyte percent cover in each plot. Only plots that had at least two bryophyte species were included (n = 98). We used ellipses to denote variances in bryophyte species composition among substrate types. We fit bryophyte depth (BD), tree seedling density (seedling), and percent canopy cover (CC) as environmental vectors to the NMDS ordination using the envfit function in the vegan package in R with 1000 permutations (Oksanen et al., 2010). To test for significant differences in species composition among substrate types, we used the adonis function in the vegan package, which is essentially a multivariate analysis of variance. We then used the pairwiseAdonis function in the vegan package to compute all pairwise comparisons among substrate types with adjusted p-values (Martinez Arbizu, 2020).
To test the light levels under and beside Hylocomium splendensmats, we used a paired t-test.
RESULTS