Dark diversity reveals importance of biotic resources and competition for plant diversity across habitats

Abstract Species richness is the most commonly used metric to quantify biodiversity. However, examining dark diversity, the group of missing species which can potentially inhabit a site, can provide a more thorough understanding of the processes influencing observed biodiversity and help evaluate the restoration potential of local habitats. So far, dark diversity has mainly been studied for specific habitats or large‐scale landscapes, while less attention has been given to variation across broad environmental gradients or as a result of local conditions and biotic interactions. In this study, we investigate the importance of local environmental conditions in determining dark diversity and observed richness in plant communities across broad environmental gradients. Using the ecospace concept, we investigate how these biodiversity measures relate to abiotic gradients (defined as position), availability of biotic resources (defined as expansion), spatiotemporal extent of habitats (defined as continuity), and species interactions through competition. Position variables were important for both observed diversity and dark diversity, some with quadratic relationships, for example, plant richness showing a unimodal response to soil fertility corresponding to the intermediate productivity hypothesis. Interspecific competition represented by community mean Grime C had a negative effect on plant species richness. Besides position‐related variables, organic carbon was the most important variable for dark diversity, indicating that in late‐succession habitats such as forests and shrubs, dark diversity is generally low. The importance of highly competitive species indicates that intermediate disturbance, such as grazing, may facilitate higher species richness and lower dark diversity.


| INTRODUC TI ON
The global biodiversity crisis represents one of the most critical challenges in the 21st century (Butchart et al., 2010;Díaz et al., 2019;Tittensor et al., 2014). Achieving conservation goals and prioritizing efforts requires appropriate metrics to quantify biodiversity and identify the factors driving the declines. The most commonly used measure is observed species richness which traditionally depends on field surveys to count the individual species. Although observed diversity can provide valuable insights into the richness of species within a given site, it does not account for the absent part of the species pool that could potentially inhabit that site considering suitable environmental conditions and biogeographic history, that is, the dark diversity (Pärtel, Szava-Kovats, & Zobel, 2011). Identifying this part of the biodiversity can provide a more thorough understanding of the processes influencing biodiversity and help evaluate the restoration potential of local habitats (Lewis, Szava-Kovats, Pärtel, & Evolution, 2016).
In contrast to observed diversity, dark diversity focuses on the portion of diversity potentially able to occur in a particular habitat type, but which is currently missing. The ultimate potential of biodiversity at a given site is mainly determined by large-scale biogeographic and evolutionary processes (i.e., species diversification and historic migration patterns) which create the set of species that can theoretically inhabit a site, defined as the species pool (Cornell & Harrison, 2014;Pärtel, Zobel, Zobel, van der Maarel, & Partel, 1996;Zobel, 2016). Although this species pool influences biodiversity patterns, the observed species are filtered by local processes such as biotic interactions, population dynamics, dispersal, anthropogenic disturbance, and stochastic events (Cornell & Harrison, 2014;Pärtel, Szava-Kovats, & Zobel, 2013;Ronk, Szava-Kovats, & Pärtel, 2015;Zobel, 2016). Dark diversity therefore reconciles the role of local (biotic interactions, abiotic filters, dispersal, stochastic events) and large-scale processes (species diversification and historic migration patterns) underlying biodiversity patterns and biological communities (Pärtel, 2014;Pärtel et al., 2011). This metric can provide insight into the determinants of missing species by helping us understand what characterizes species that are often missing or sites missing many species. Quantifying dark diversity patterns, in combination with observed diversity patterns, can allow researchers to better understand the mechanisms and processes acting on individual populations or entire communities .
The potential value of dark diversity to guide conservation and restoration planning has been demonstrated for mammals (Estrada, Márcia Barbosa, & Real, 2018), sharks (Boussarie et al., 2018), and fungi , but most studies have considered plants (Bennett et al., 2016;Moeslund et al., 2017;Ronk, De Bello, Fibich, & Pärtel, 2016;Ronk et al., 2015). Dark diversity has also proven valuable in understanding plant diversity patterns, such as determining that vascular plant dark diversity across Europe follows a latitudinal gradient (Ronk et al., 2015). So far, the plant traits likely to increase a species' probability of being part of the dark diversity include stress intolerance, height, adaptation to low light and nutrient levels, and producing fewer and heavier seeds Riibak et al., 2015;Riibak, Ronk, Kattge, & Pärtel, 2017). Moreover, understanding the ecological processes governing plant dark diversity contributes to understanding biodiversity in general, that is, plants are bioindicators of their abiotic environment and anthropogenic impact (Bartelheimer & Poschlod, 2016), and they form the living and dead organic carbon pools and biotic surfaces that are the niche space for all other taxonomic groups DeAngelis, 2012), and vascular plants predict biodiversity across environmental gradients and broad taxonomic realms, and are related to the occurrence of regionally red-listed species of other taxa .
Nevertheless, as a relatively new concept, more research is required to establish its full potential and to understand the ecological processes governing dark diversity across plant communities.
Most dark diversity research has ignored the variability between types of habitats and have mostly been restricted to narrow sets of variables and specific habitats (Riibak et al., 2015) or large-scale landscapes (Ronk et al., 2015(Ronk et al., , 2016, with no studies examining how dark diversity varies across large environmental gradients or the importance of local conditions and biotic interactions. Applying dark diversity within one habitat type may produce adequate results, for example, as seen for grasslands (Riibak et al., 2015), but biodiversity varies greatly across ecosystems and is highly dependent on the habitat and region of interest . Dark diversity can be used to derive community completeness, a relativized biodiversity index, which has been proposed as a valuable tool for facilitating biodiversity comparisons irrespective of regions, ecosystems, and taxonomic groups (Pärtel et al., 2013). The community completeness index can be defined, in general terms, as the proportion of species from the regional species pool which have dispersed to and established at a site after abiotic and biotic filtering (Pärtel et al., 2013). Since patterns in observed species richness may mimic patterns in dark diversity (e.g., exhibit a strong latitudinal gradient; Aning, 2017;Pärtel et al., 2011;Ronk et al., 2015;Zobel, 1997), community completeness can provide a different aspect of biodiversity as it accounts for the variation in species pool size and expresses biodiversity on a relative scale (Pärtel et al., 2013). For instance, completeness exhibited no relationships to latitudinal gradients, but strong relations to anthropogenic disturbance (higher completeness in areas with lower disturbance) for fungi , plants (Ronk et al., 2015(Ronk et al., , 2016, and birds (Cam, Nichols, Sauer, Hines, & Flather, 2000). Comparing the environmental processes influencing these biodiversity measurements can provide valuable information for better prioritization of resources and understanding patterns of biodiversity. However, despite observed species richness and its determining factors being relatively established, dark diversity and its completeness counterpart are new methodologies, and as such, have not been well investigated, and the factors influencing them are not fully understood.
Determining the set of species that can theoretically inhabit a site, the species pool, is typically estimated using species co-occurrence patterns with Beal's smoothing which assumes species with shared ecological requirements and biogeographic history will have similar likelihoods of being present at a particular site (Beals, 1984; de Bello et al., 2012;Lewis et al., 2016;McCune, 1994;Münzbergová & Herben, 2004). Since particular pairs of species will be less likely to co-occur if they compete for the same resources, this approach also assumes to account for competitive interactions which are a major factor influencing species occurrence patterns (Cornell & Harrison, 2014;de Bello et al., 2012), especially in plant communities (Riibak et al., 2015). Although the Beals co-occurrence matrix approach is expected to account for biotic interactions in the regional pool estimates, competitive exclusion can lead to lack of species at sites, thereby increasing dark diversity and decreasing species richness. To examine this further, we used Grime's plant life strategies (competitor/stress-tolerator/ruderal, C-S-R (Grime, 1979)) to quantify the occurrence of competitive species in plant communities (Ejrnaes & Bruun, 2000); that is, when a site has a high Grime C value, it is dominated by few, but abundant, competitive species.
One way to consider the roles these factors play in dark diversity measurements can be provided with the recently developed ecospace framework . Ecospace divides the environmental causes of variation in species diversity into three main domains: (1) the position in environmental hyperspace (position), (2) the availability and variation of biotic resources (expansion), and (3) the spatiotemporal extent of habitats (continuity). This framework can be used to quantify and examine the roles of, for example, environmental filtering (position), as well as succession and human disturbance (continuity) on dark diversity. Ecospace also recognizes the role of vegetation structure and diversification of organic matter (expansion), as a contributing factor of biodiversity, bringing to light often-ignored trophic interactions that exist between taxa . The importance of expansion on biodiversity has been illustrated by a recent study where vegetation structure was found to influence biodiversity across trophic groups of plants, animals, fungi, and bacteria (Penone et al., 2019). This framework unites theories such as niche theory, island biogeography theory, and a suite of community assembly theories into one framework for further development of a general theory of terrestrial biodiversity .
In this study, we investigate the importance of local environmental conditions and competition for dark diversity and completeness in plant communities across habitats and compare with results for observed plant richness. We ask what the relative importance of the ecospace dimensions is across habitats, and when habitat differences are accounted for in the diversity measure. We discuss how dark diversity can contribute to new aspects for informed conservation and management.

| Study sites
Our data stem from Biowide (www.biowi de.dk), a nationwide survey of biodiversity in Denmark (Brunbjerg, Bruun, Broendum, et al., 2019). A total of 130 study sites (40 m × 40 m) were evenly distributed across five geographic regions in Denmark with a minimum distance of 500 m between sites (Figure 1a). Each site is sampled in four 5-m circle plots ( Figure 1b). Sampling was designed to evaluate the ecospace framework, stating that biodiversity varies according to abiotic conditions, buildup and diversification of organic resources, and spatiotemporal continuity . The sites were stratified to represent the main variation in major environmental gradients. Thirty sites were cultivated habitats and 100 sites natural and semi-natural habitats. The cultivated subset was stratified to represent major land-use types, and the natural subset was stratified to represent major gradients in soil fertility, soil moisture, and successional stage. Saline and fully aquatic habitats were deliberately excluded, but temporarily inundated depressions, as well as wet mires and fens, were included. The final set of 24 environmental strata consisted of the following six cultivated habitat types: three types of fields (rotational, grass leys, and set aside) and three types

| Plant species richness
Vascular plant species richness was inventoried in 5-m circular plots at each site (four plots at each site) by trained botanists during the summer 2014 and spring 2015 to account for variations in phenology. We removed all subspecies, hybrids, variations, and neophytes (i.e., species that are not considered a natural part of the vegetation given their history and dispersal ability; see appendix tables 6-8 in Buchwald et al., 2013). The nomenclature in this study follows the species checklist Denmark database from https://allea rter-datab asen.dk. The data contained a total of 580 species of vascular plants in the 115 sites spanning from open habitats to shrubs and late-succession forests.

| Explanatory variables
According to ecospace, the position variables included in the model were soil moisture index (SMI) and soil fertility index (SFI). For each site, SFI represents the predicted value from the best linear model (of all sites) of site mean Ellenberg N (plant-based bioindication of nutrient status (Ellenberg et al., 1991)) as a function of soil Ca, leaf N, leaf NP, and soil type. We calculated a soil moisture index for each site using the predicted values from the best linear model soil organic carbon as % soil C in 0-10 cm soil layer (g/m 2 average of four soil samples taken at each site) as described in Brunbjerg, Bruun, Broendum, et al. (2019), and (e) vegetation heterogeneity derived from Lidar. The latter included both canopy height variation (variance of the 90th percentile for points >3 m within the site) and shrub layer height variation (variance of the 90th percentile for returns between 30 cm and 3 m [Brunbjerg, Bruun, Dalby, et al., 2019]).
Spatial continuity included landscape characteristics: (a) fraction of intensive fields within a 500-m buffer of the site and (b) the fraction of natural habitats in a 1 km × 1 km grid from a national mapping, interpolated using Spline in ArcGIS 10.2.2 (weight 0.5, number of points 9 [Ejrnaes et al., 2014]). Temporal continuity was estimated directly for each study site as time since major land-use change: For each site, a temporal sequence of aerial photographs and historical maps was inspected starting with the most recent photographs (photographs from 2014, 2012, 2010, 2008, 2006, 2004, 2002, 1999, 1995, 1968, 1954, 1945) and ending with historical maps reflecting land use in the period 1842-1945. Temporal continuity (the year in which a change could be identified) was reclassified into a numeric 4-level variable: 1: 1-14 years, 2: 15-44 years, 3: 45-135 years, and 4: >135 years (Brunbjerg, Bruun, Broendum, et al., 2019).
Lastly, to examine the assumption that competitive exclusion can lead to lack of co-occurrence and hence increased dark diversity, we used Grime's plant life strategies to quantify the importance of local interspecific competition for establishment of species. Grime's plant strategy theory (competitor/stress-tolerator/ruderal, C-S-R) states a three-way trade-off between life strategies that facilitate competition for resources (competitive strategy), survival in stressful environments, for example, with high salinity, flooding, and drought (stress tolerance), and survival in disturbed environments (ruderalism; Grime, 1979). The original C-S-R species strategies were recoded as numeric values for each plant species where a total of 12 points were distributed to the different strategies as described in Ejrnaes and Bruun (2000). To

| Regional pool, dark diversity, and completeness
All statistical analyses were performed in R version 3.5.3 (R Core Team, 2019). To calculate regional pools, we used vegetation data from the 5-m circular Biowide plots (four at each site) as well as an additional dataset of plant inventories in 5-m circles from the national monitoring program (Danish Nature Agency, 2016). This dataset includes 52,362 plots with more than five species recorded by trained botanists and is added to increase the co-occurrence matrix for robustness in the Beal's calculations of regional pools (see below). We did not include species-poor plots, that is, those with less than five observed species, resulting in 448 plots from Biowide and 52,362 plots from the additional dataset. The regional pool was calculated using the Beals index (Beals, 1984), as recommended by Lewis et al. (2016).
The Beals index represents the probability that a particular species will occur within a given site based on the assemblage of co-occurring species (Beals, 1984;McCune, 1994;Münzbergová & Herben, 2004).
We calculated the Beals index using the "beals" function in the "vegan" package (Oksanen et al., 2017). The threshold for including a particular species in the regional species pool is recommended to be the 5th percentile of the Beals index value for the species in question (Gijbels, Adriaens, & Honnay, 2012;Ronk et al., 2015). Preceding the calculation of each threshold, the lowest Beals index value among plots with occurrence of the species in question was identified, and all plots having values below that minimum were not considered.
Analyses were done at the site level (n = 115) by creating a site regional pool combining the four plot regional pools at each site. In addition to the four plots, plants had been inventories in the whole site, and therefore observed species in the site, but not included in the regional pools (n = 2) were added to the regional pools to ensure that site regional pool included all observed species. Then, dark diversity was calculated for each site as the difference between the regional pool and the observed species richness (Pärtel et al., 2011) and completeness was calculated following Pärtel et al. (2013) using the formula ln(observed richness/dark diversity). Because dark diversity is relative it may not be suitable for comparison across habitats (Scott, Alofs, & Edwards, 2011), and completeness is suggested as a possible alternative (Pärtel et al., 2013). In this study, we found that completeness was highly correlated with observed species richness ( Figure 2), so instead, we propose a dark diversity measure corrected for habitat types using the residuals of a model of dark diversity as a function of the habitat type. To represent habitat type, we conducted a supervised classification of nine a priori defined training classes and using ordination gradients from an NMDS ordination based on all identified species, that is, plants, fungi, and insects from the Biowide project. The classification resulted in eight classes spanning gradients in succession (early, late), moisture (wet, dry), and nutrients (rich, poor). The habitat types explained 51% of the variation in dark diversity.

| Statistical analyses
Soil pH, litter mass, organic carbon, organic matter, and vegetation variation were log-transformed, and all explanatory variables were standardized. We conducted model selection by first testing for multicollinearity using the variance inflation factor (VIF; Zuur, Ieno, & Elphick, 2010). We removed canopy height variation and organic matter resulting in a maximum VIF of 2.9. The remaining variables were used as explanatory variables (linear terms) in a set of generalized linear models (GLMs) with Poisson distribution for the count data for dark diversity and plant richness, and normal distribution for residual dark diversity, but as models for plant richness were overdispersed, we chose negative binomial models for this response instead. To avoid spurious correlation in the models, we excluded variables with no hypothesized relationship and constrained the response direction and shape to ecologically plausible responses (Burnham & Anderson, 2002;Zuur et al., 2010). In general, more resources, more diverse resources, more environmental variation, and increasing temporal and spatial continuity are all hypothesized to increase plant species richness. Exceptions to this are litter mass, intensive fields, and competitive ability, which are expected to decrease plant species richness. Generally, we expect the opposite relationships with dark diversity. Following this approach, we excluded negative effects of expansion, continuity, and heterogeneity variables on plant species richness. We did not include interaction terms because of lacking plausible ecological hypotheses. To allow for nonlinear relationships for position variables corresponding to the intermediate disturbance hypothesis (Connell, 1978;Townsend, Scarsbrook, & Dolédec, 1997) and intermediate productivity hypothesis (Fraser et al., 2015), we used AIC (Burnham & Anderson, 2002) to evaluate whether inclusion of quadratic terms for the variables SMI, SFI, light, soil pH, and bare soil improved the model fit (using delta AIC < 2 as a rule), in which case relevant quadratic terms were added. Subsequently, we checked again the responses' direction and excluded ecologically implausible responses (Burnham & Anderson, 2002). We then dropped remaining variables sequentially based on AIC using the function drop1() in the lme4 package (Bates, Sarkar, Bates, & Matrix, 2007). The final regression models were tested for overdispersion and evaluated by visual inspection of residual plots and for spatial autocorrelation using Moran's I in the ape package (Paradis, Claude, & Strimmer, 2004).

| RE SULTS
The species richness per site ranged from 8 to 127 species, and dark diversity ranged from 84 to 243 species. Completeness and plant species richness were highly positively correlated (r s = .95, Figure 2), and completeness was therefore excluded from further analysis. Dark diversity was less correlated with plant species richness (r s = .17, Figure 2). The final regression models explained between 14% and 65% of the variation in dark diversity, residual dark diversity, and species richness (Table 1). We found position variables to be important for dark diversity and plant species richness (Figures 3 and 5). Soil moisture invoked a unimodal response in dark diversity but a bimodal response in species richness with a low at inter-

| D ISCUSS I ON
In this study, we found position variables to be important for both plant richness and dark diversity. Once abiotic and successional differences between habitat types were accounted for, that is, residual dark diversity, the position variables were no longer significant indicating that the habitat classification used here is largely driven by abiotic gradients. Plant richness was highest at intermediate conditions of soil fertility and pH, corresponding to the intermediate productivity hypothesis, which states that few species can tolerate the environmental stresses at low productivity and a few highly competitive species dominate at high productivity (Fraser et al., 2015). Species richness increased with pH, possibly corresponding to the generally large regional species pool in calcareous habitats, and aligns with previous research indicating that plant diversity has a strong positive association with soil pH in temperate and boreal regions (Pärtel, 2002;Pärtel, Helm, Ingerpuu, Reier, & Tuvi, 2004). The unimodal relationship between dark diversity and soil moisture may reflect that hydrological gradients can be strong environmental filters determining plant communities (Fraaije et al., 2015;Silvertown, Araya, & Gowing, 2015;Valdez, Hartig, Fennel, & Poschlod, 2019) with more distinct communities at the extremes than at intermediate soil moisture, that is, specific adaptations for waterlogged and very dry soil are required (Ernst, 1990). Therefore, fewer coincidental species may appear in these extreme habitats compared with habitats of intermediate soil moisture, resulting in a low estimated regional pool and lower dark diversity at the extremes.
Expansion of ecospace seems more important for plant species richness and less important for dark diversity as only organic F I G U R E 2 Spearman rank correlations and dot plots of site regional pool, dark diversity, residual dark diversity, plant species richness, and completeness. The red line in the plot shows a loess smoothing line Completeness carbon is significant for dark diversity. This is somewhat unexpected as many of the hypotheses for associations between expansion and richness would also apply inversely to dark diversity.
For example, variation in shrub height was positively correlated with species richness corresponding to a general positive effect of vegetation heterogeneity on species richness by increasing available niche space and providing grazing refuges (Stein, Gerstner, & Kreft, 2014). In our experience, shrub height variation also reflects land-use history, that is, intensive, mechanical land use causes homogenization and general low vegetation that concurs with low species richness. Concordantly, shrub height variation should be negatively related to dark diversity; however, it is insignificant.
Species richness was negatively related to litter mass corresponding to a general negative effect of litter on germination, establishment, richness, and biomass of plants (Xiong & Nilsson, 1999). This effect, however, does not seem to apply to dark diversity, where litter mass is insignificant. We found a negative relationship between bare soil and species richness, indicating that although bare soil is expected to increase the survival at germination (e.g., Roth, Seeger, Poschlod, Pfadenhauer, & Succow, 1999) and thereby affect community composition and increase plant richness, it may also reflect a trade-off where bare soil means fewer growing individuals and species. The only expansion variable to influence dark diversity, residual dark diversity, and species richness was soil organic carbon, influencing all three diversity measures negatively.
While organic carbon may not be a biotic resource for plants in accordance with ecospace, it is related to long-term stability in late-successional stages as soil organic carbon accumulates over time (Luyssaert et al., 2008). While long-term stability and hence time for establishment imply lower dark diversity, high organic carbon also coincides with species-poor habitats, such as forests, mires, bogs, and heaths.
For the first time, we attempted to investigate the effect of competition on dark diversity. The importance of competition for dark diversity hints that the regional pool estimate using Beals does not fully account for biotic interactions although it is expected from the co-occurrence approach inherent in the Beals index (Cornell & Harrison, 2014;de Bello et al., 2012;Riibak et al., 2015). However, scale may also play a role. The plant communities used in this study cover broad environmental gradients TA B L E 1 Regression models of dark diversity (DD), residual dark diversity (Residual DD), and plant species richness (PlantRich) (including the additional dataset of 52,362 plots), and when estimating the regional pool from the co-occurrence matrix, abiotic filtering may therefore be more prevalent over biotic filtering (Kraft et al., 2015). High competitive exclusion as a biotic filtering of the regional pool may still occur at sites with relatively few but abundant, competitive species and implies that other species are missing. This corresponds to our findings of both decreased species richness and increased dark diversity in communities characterized by competitive species (high community mean GrimeC).
Excluded species from the regional pool can colonize where communities dominated by competitive species are disturbed (e.g., by grazing, erosion, or other). On the contrary, competition was not significant for the residual dark diversity. It is possible that accounting for the variation in habitat types not only accounts for variation in position but also the inherent competitive strategies of the species in the habitat types, for example, late-successional stages with long-term stability are generally dominated by competitive species. Recent review and opinion papers (Cadotte & Tucker, 2017;Kraft et al., 2015) have highlighted that biotic interactions of species may be related to environmental gradients.
Residual dark diversity seems to decrease with the fraction of natural land in the surroundings of the sites possibly by influencing local processes, that is, landscapes with high density of nature are likely to have a higher local pool of species, increased dispersal, increased species survival, metapopulation structures, and less negative edge effects of intensive land use . It is possible that this effect is otherwise masked by the large effect of habitat type and position variables in both dark diversity and observed diversity, making the effect of natural landscapes only visible once this other variation has been accounted for.
This study also shows that different diversity measures contribute with different aspects to better understand drivers of diversity. For example, the negative effect of organic carbon on plant species richness could indicate that carbon storage leads to loss of species richness, whereas the effect on dark diversity actually indicates that fewer species are missing when carbon storage is high. Here, we compared dark diversity and completeness, thought to be less dependent on habitat types (Pärtel et al., 2013). However, we found that completeness was highly correlated with observed species richness and therefore added no new information to our study of plant community diversity aspects.
We therefore corrected for differences in regional pool through a residual dark diversity to provide a measure comparable across habitats.  With global biodiversity rapidly decreasing, it is vital to understand the drivers of biodiversity to prioritize conservation and make management more efficient. In this study, besides the ecospace position variables, competition seems to be the strongest predictor of plant richness. Conservation management focusing on intermediate disturbance such as grazing can disturb competitive communities making room for more species and thereby decrease dark diversity.
Besides ecospace position, organic carbon was the most important variable for both dark diversity measures indicating that longterm stability in late-successional habitats decreases dark diversity.
Examining the influencing factors of different measures of biodiversity can lead to better decision-making in the future conservation of the world's biodiversity.

ACK N OWLED G M ENTS
We sincerely thank Aage V. Jensen Nature Fund for financial support to CF, AKB, JM, LD, KC, and JV for the project "Dark Diversity in nature management." The Biowide project and REJ were supported by a grant from the Villum Foundation (VKR-023343). MP has been supported by the Estonian Ministry of Education and Research (IUT20-29), and the European Regional Development Fund (Centre of Excellence EcolChange). The authors declare no conflict of interest.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
Plant species lists from 5-m plots and environmental variables: Dryad https://doi.org/10.5061/dryad.9kd51 c5d5.