Modelling species responses with HMSC
We applied HMSC from the Hmsc R package (Tikhonov et al. 2020) to fit a joint species distribution model to plant community data simultaneously including information on traits (genome size), environmental covariates, and phylogenetic relationships in a single model (Pollock et al. 2012; Warton et al. 2015). We included MAP, MAT, altitude, soil N, and soil P as fixed effects (i.e., predictors) and used the sampling site as the sampling unit and the random effect. HMSC estimates species responses (slope parameters) to environmental covariates and uses these responses as the species’ functional niche to estimate the strength and direction in which these niches moderate multiple species responses to environmental filtering. To account for the phylogenetic correlations among all the species, we included a phylogenetic correlation matrix (construction details of the phylogenetic trees see above) in the model’s covariance structure. HMSC model can indicate the strength of the phylogenetic signal based on parameter ρ (varying from 0 to 1, a higher value indicating stronger phylogenetic signal) (Ovaskainen et al. 2017).
We fitted the model to the plant community (with a probit-link for the presence/absence data) with Bayesian inference, using the posterior sampling scheme (Ovaskainen et al. 2016). We used four Markov Chain Monte Carlo (MCMC) chains, each of which consisted of 15,000 iterations, out of which we discarded the first 5,000 as the burn-in and thinned the remaining by 10 to yield in total 1000 posterior samples per chain. We assessed the convergence of the MCMC chains by examining the distribution of the potential scale reduction factor over the parameters that measure the responses of the species to the fixed effects included in the model. The MCMC convergence of the HMSC models was satisfactory as the potential scale reduction factors (psrf) are close to one (the psrf of the 𝛽-parameters that measure species response to environmental covariates (Ovaskainen et al. 2017) were on average 1.16). We examined the explanatory and predictive powers of the HMSC models through species-specific Tjur’s R2 values (Tjur 2009). Twofold cross-validation was performed to assess the predictive power of the model.
Based on the output of HMSC, we first characterized species’ responses to each of the environmental variables as a mean parameter and its 95% posterior probability. Second, we predicted the community weighted mean genome size based on 161 species considered in the model as responses to the included environmental covariates. Note that genome size was log-transformed. Third, we used HMSC to capture the species-to-species associations (positive or negative co-occurrences) with latent variables. We considered statistical support for these associations based on 90% posterior probabilities. Association patterns were classified as positive or negative, with this result used to estimate the frequency of two taxa co-existing (or not) compared to the frequency expected based on shared habitats or random occurrences.