Diversity metrics
Species diversity
We measured community richness using: (1) species richness (SR), (2) the Shannon-Wiener diversity index (H’) (Shannon and Weaver 1949), and (3) Smith and Wilson’s (1996) index of evenness (Evar). Species presence/absence data was collected for all sub-plots, so to avoid pseudo-replication we calculated a single species richness measure for each plot by averaging the species richness of the nine sub-plots. In addition to calculating the overall species richness of each plot, we calculated the species richness of shrubs (25 species total), herbs (102 species total), and graminoids (52 species total) separately. Our other two measures of community richness required species abundance data, which was only collected for one sub-plot in each plot, so H’ and Evar measurements at the plot-level were taken from a single sub-plot.
Phylogenetic diversity
Phylogeny construction
To study the effect of reindeer grazing on phylogenetic diversity, we constructed a regional vascular plant phylogeny. After reducing all taxonomy assignments to the species-level, which involved re-assigning subspecies, hybrids, and synonymous species, we identified 145 unique species from our samples to include in the regional phylogeny. DNA sequences for these species for the genes matK and rbcL were collected from GenBank (Supplementary material Appendix 1, Table A1). If sequence information was not available for one or both of matK and rbcL for a species, we used the closest-relative within the same genus with available matK and/or rbcL sequence information. Sequences for the identified species or close relative were available for 138/145 species for rbcL and 140/145 for matK (Supplementary material Appendix 1, Table A1). Two species did not have available sequence information for either gene and did not have a close relative that could be substituted:Cerastium cerastoides and Carex parallela . C. cerastoides was manually inserted at the base of the Cerastiumclade (Scheen et al. 2004) and C. parallela was manually inserted into a clade with Carex dioica (Lipnerová et al. 2013).
To ensure the correct reading frame prior to sequence alignment, all sequences were translated to amino acids in ExPASy (Gasteiger et al. 2003). Alignments were completed in MEGA7 (Kumar et al. 2016) using MUSCLE with the nucleotides for the coding regions of plant plastids and default settings. Non-informative gaps were manually removed from the matK alignments and excess lengths were trimmed from the ends of both genes. We identified GTR+G+I as the best model for rbcL and GTR+G as the best model for matK using the Model Selection tool in MEGA7.
Bayesian inference of trees was performed using MrBayes (Huelsenbeck and Ronquist 2001) with two partitions: one for matK and one for rbcL. We ran our model for 100 000 000 generations with 4 chains at a temperature of 0.2 and a stop value of 0.01 for convergence. The sample frequency was set to every 1000 generations with a burn-in fraction of 0.25 and the trees were dated by constraining seven nodes (Table 1; Bell et al. 2010). The analysis converged at 10 304 000 generations. We illustrate the posterior with the majority rule consensus tree, created by collapsing clades with posterior probabilities less than 50% to polytomies (Supplementary material Appendix 1, Figure A1).
Phylogenetic structure
We used mean pairwise distance (MPD) as our measure of phylogenetic structure. MPD is more sensitive to changes in distantly related taxa than is the mean nearest taxon distance (Webb et al. 2002). To prevent bias resulting from calculations due to distantly related species (e.g.Lycopodium species), we used only the angiosperms (133/145 sampled species) to calculate MPD. We calculated a standardized measure of MPD with the aid of the function ses.mpd , in the R package “picante” by comparing the observed phylogenetic community structure to a specified null model with a randomized community structure (Kembel et al. 2010, R Core Team 2019). Using taxa that were identified to the species-level, we calculated both (1) a presence/absence-based measure of MPD (calculated at the sub-plot level and then averaged within each plot) and (2) an abundance-based measure (calculated for a single sub-plot in each plot), by weighting the pairwise distances by the product of the relative abundance of each species in each pair. Both MPD metrics were calculated for every tree in the posterior sample (created by merging the two runs and removing the 25% burn-in: final n = 15,458 trees) and averaged to produce a single measure for each plot. For our null model we used the independentswap algorithm (with 1000 iterations per run and 999 runs), which randomizes the community data matrix while maintaining species occurrence frequency and sample species richness (Gotelli 2000). We used the p-value, or quantile, of observed MPD vs. the MPD of null communities as our standardized response variable, as this metric is less biased than the more common Net-Relatedness Index (Vamosi et al. 2014). Our measure describes the rank of the observed phylogenetic dispersion relative to the distribution produced by the null model. A value of 0 corresponds to a community that is more clumped than any of the null communities, a value of 0.5 corresponds to a community that has a median dispersion relative to the null communities, and a value of 1 corresponds to a community that is more dispersed than any of the null communities. We refer to these metrics as phylogenetic dispersion (from the species presence/absence data) and abundance-weighted phylogenetic dispersion (from the relative species abundance data).