Vascular plant species richness and bioindication predict multi-taxon species richness

1. Plants regulate soils and microclimate, provide substrate for heterotrophic taxa, are easy to observe and identify and have a stable taxonomy, which strongly justi - fies their use as indicators in monitoring and conservation. However, there is no consensus as to whether plants are strong predictors of total multi-taxon species richness. In this study, we investigate if general terrestrial species richness can be predicted by vascular plant richness and bioindication. 2. To answer this question, we collected an extensive dataset on species richness of vascular plants, bryophytes, macrofungi, lichens, plant-galling arthropods, gastro - pods, spiders, carabid beetles, hoverflies


| INTRODUC TI ON
Given the global biodiversity crisis, the need for establishing causes for spatial and temporal variation in biodiversity is acute (Ceballos, Ehrlich, & Dirzo, 2017;Hill et al., 2016).However, the majority of species worldwide are still undescribed and nowhere on Earth are all the locally resident species known.As a result, the causes for spatial variation in biological diversity represent a perpetual challenge for ecological science (Pennisi, 2005), with few generally accepted causal mechanisms and models (e.g., DeMalach, Zaady, & Kadmon, 2017;Grace, Adler, Stanley Harpole, Borer, & Seabloom, 2014;Pärtel, Bennett, & Zobel, 2016).Moreover, we also lack costeffective, validated methods for assessing biodiversity.
Vascular plants are the dominant primary producers of terrestrial ecosystems and plants are quite accurate indicators of the abiotic environment, in which they grow.Here, we test whether plant community composition and richness may be used to predict the overall biodiversity through bioindication of abiotic conditions and human impact as well as biotic diversification of organic matter into separate pools of live biomass as well as litter, wood and dung (ecospace sensu Brunbjerg, Bruun, Moeslund, et al., 2017).
The intractability of total species surveys, has motivated the use of surrogate species in conservation planning (Margules & Pressey, 2000;Sarkar & Margules, 2002), with the underlying assumption that species richness correlate between taxonomic groups (Gaston, 1996).Surrogate species are assumed to reflect the distribution of other species or taxonomic groups, but also to indicate the occurrence of habitats and species of high conservation value (Pearman & Weber, 2007).Much research has focused on testing surrogacy and selecting the best taxa (reviewed in Rodrigues & Brooks, 2007).It has generally been found that correlations in species richness across taxa vary depending on spatial scale (grain and extent), geographic location (Hess et al., 2006), and taxonomic focus (e.g., Wolters, Bengtsson, & Zaitsev, 2006).Overall, biodiversity surrogacy studies have shown only weak predictive power (Rodrigues & Brooks, 2007;Su, Debinski, Jakubauskas, & Kindscher, 2004).Similarity in community composition shows more convincing results than species richness.This may be because species composition exhibits a stronger relationship to environmental gradients than species richness (Prober et al., 2015;Su et al., 2004).For selecting areas of conservation interest, multi-taxon surrogacy has been proposed as a more robust measure of biodiversity than single-taxon surrogacy (Smith-Patten & Patten, 2015).Environmental surrogates have also been tested, but found to be less useful for prediction than cross-taxon surrogates (Rodrigues & Brooks, 2007).
Plants are included in biodiversity monitoring programmes for several good reasons: Plants are sessile and reflect conditions at the place of observation, plants are less seasonal and their detection is less dependent on weather conditions than are fungi and arthropods, plants occur in most ecosystems, and skilled field botanists are generally available.Despite their wide use, the evidence for using plants as surrogates for multi-taxon biodiversity is equivocal (Myšák & Horsák, 2014;Saetersdal et al., 2004;Westgate, Barton, Lane, & Lindenmayer, 2014;Wolters et al., 2006).Complex metrics representing habitat quality based on weighted measures of vegetation structure (e.g., native plant species richness, number of trees with hollows, and total length of fallen logs), plant species richness and functional diversity, have also been suggested to work as surrogates for overall biodiversity, but with limited success (Hanford, Crowther, & Hochuli, 2017;Kwok, Eldridge, & Oliver, 2011).Despite the moderate support, plant-based monitoring programmes and conservation guidelines remain a common practice, even at supranational levels.For example, in the EU Habitats Directive (1992), plants are implicitly assumed to work as indicators for both habitat types (socalled Annex 1 habitats) and their conservation status.Moreover, averaging plant indicator values (e.g., Ellenberg Indicator Values, Ellenberg et al., 1991) is commonly used in vegetation studies to assess local conditions (e.g., Diekmann, 2003).The validity of plantbased bioindication has been confirmed by direct measurement of the environmental conditions and by plant growth experiments (e.g., Bartelheimer & Poschlod, 2016;Schaffers & Sýkora, 2000).
Our approach to bioindication follows the ecospace framework (Brunbjerg, Bruun, Moeslund, et al., 2017).Since plants can be used as indicators of the abiotic environment, they can describe the ecospace position of a site.With regard to ecospace expansion, i.e., the differentiation of organic matter, each different plant species constitute a potential substrate for specialised insects and fungi (Basset et al., 2012;Brunbjerg, Bruun, Moeslund, et al., 2017;Strong, Lawton, & Southwood, 1984;Zhang et al., 2016).While the species richness responses to ecospace position along environmental gradients may vary among taxonomic groups, we generally expect ecospace expansion by plant species richness to have a positive effect, at least on the richness of heterotrophic taxa.
Plants are highly responsive to land-use change, which usually involves replacement of natural vegetation by crops and weeds, a process generally considered a major cause of biodiversity loss (Lehsten et al., 2015;Pimm et al., 2014).Effects may be detectable in plant communities for decades or even longer (Gustavsson, Lennartsson, & Emanuelsson, 2007;Hermy & Verheyen, 2007).We expect a plant-derived land-use intensity indicator to be useful for prediction of multi-taxon species richness, especially when used in combination with a plant-derived indicator of abiotic conditions and plant species richness.
Here, we put the value of a plant species list to a test.We use a comprehensive dataset of 130 sites, each 40 × 40 m, sampled for richness of plants and a range of other taxa including DNA derived OTUs (operational taxonomic units), collectively spanning the major environmental variation in terrestrial habitats within Denmark (Figure 1), situated in the Northern temperate part of Europe.These data allow us to investigate the following questions: 1. Can plant species richness be used to predict species richness of other taxa across habitat types?
2. Can we improve this prediction by adding environmental bioindication based on the observed plant species in a multiple regression?
Our study is the first to comprehensively validate the common use of vascular plants in terrestrial conservation planning and monitoring, which has been incorporated-based on anecdotal evidence and tradition-into national and supranational legislation.

| M ATE R I A L S A N D M E TH O D S
We selected 130 study sites (40 × 40 m) distributed across five geographic areas in Denmark (Figure 1).Within each area, sites were placed in three clusters for logistical reasons, but with a minimum distance of 500 m between sites to reduce spatial covariance.Site selection was stratified according to primary environmental gradients.We allocated 30 sites to cultivated habitats and 100 sites to natural habitats.The cultivated subset was stratified according to major land use type and the natural subset was selected amongst uncultivated habitats and stratified according to gradients in soil fertility, soil moisture, and successional stage.We deliberately excluded saline and aquatic habitats, but included temporarily inundated depressions as well as mires and fens.The final set of 24 habitat strata consisted of the following six cultivated habitat types: Three types of fields (rotational, grass leys, set aside) and three types of plantations (beech, oak, spruce).The remaining 18 strata were natural habitats, constituting all factorial combinations of: Fertile and infertile; dry, moist and wet; open, tall herb/scrub and forest.These 24 strata were replicated in each of the five geographical areas.We further included a subset of 10 perceived hot spots for biodiversity in Denmark, selected subjectively by public voting among active naturalists in the Danish conservation and management societies, but restricted so that each region held two hot spots.See Brunbjerg, Bruun, Brøndum, et al. (2017) for more details on site selection and stratification.

| Collection of biodiversity data
The field inventory aimed at an unbiased and representative assessment of the multi-taxon species richness in each of the 130 sites.We collected data on vascular plants, bryophytes, lichens, macrofungi, arthropods and gastropods.Due to the limited size of each site, we excluded vertebrates from our study, as we did not expect the richness of these to be reliably covered.We conducted a complete inventory of vascular plants and bryophytes.For the remaining taxa, which are more demanding to find, catch, and identify, we aimed at collecting a reproducible and unbiased sample through a standardised effort.We carefully examined each site for gastropods, lichens, plant-galling and mining arthropods/fungi (plant-galling arthropods hereafter) and macrofungi, actively searching contrasted microhabitats and substrates (soil, herbaceous vegetation and debris, wood, stone surfaces, and bark of trees up to 2 m).We used a standard set of passive traps to survey insects (pitfall traps, meat-baited and dung-baited traps, yellow pan traps and Malaise traps) during periods of standard length and timing.Biodiversity survey methods are detailed in (Brunbjerg, Bruun, Brøndum, et al., 2017).

| OTU richness from DNA metabarcoding
DNA metabarcoding of e.g., a soil sample is increasingly used to assess biological communities (Taberlet, Coissac, Pompanon, Brochmann, & Willerslev, 2012).Sequences are traditionally grouped into ecologically meaningful OTUs, which are then used as richness estimates (Bálint et al., 2016;Frøslev et al., 2017).In order to extend our biodiversity assessment to organisms assumed to be poorly represented in traditional surveys, we included OTU richness of fungi and general eukaryotes from soil samples and of arthropods from Malaise traps (called fungal OTUs, eukaryote OTUs, and Malaise OTUs, respectively).
We collected soil from all sites and subjected it to metabarcoding through DNA extraction, PCR amplification of genetic marker regions (DNA barcoding regions) and massive parallel sequencing on the Illumina platform as described in (Brunbjerg, Bruun, Brøndum, et al., 2017).The soil sampling scheme included the mixing of 81 soil cores from each site in an attempt to get a representative sample.
For this study, we used sequencing data from genes amplified with primers targeting fungi and eukaryotes.For eukaryotes, we used the primers 18S_allshorts (Guardiola et al., 2015(Guardiola et al., , 2016) ) with a slight modification of the forward primer (TTTGTCTGGTTAATTCCG) to exclude fungi.For fungi, we amplified the ITS2 region with primers gITS7 (Ihrmark et al., 2012) and ITS4 (White, Bruns, Lee, & Taylor, 1990).Furthermore, we extracted DNA from the ethanol of Malaise traps and subjected it to sequencing.For this, we amplified a region of the CO1 gene with primers ZBJ-ArtF1c and ZBJ-ArtR2c targeting primarily arthropods (Zeale, Butlin, Barker, Lees, & Jones, 2011), and subjected it to a sequencing approach similar to the other markers.
We extracted, amplified, and sequenced DNA from two independent trapping periods for each site and combined these to obtain a single richness measure.
The bioinformatic processing of the sequence data followed the strategy outlined in (Brunbjerg, Bruun, Brøndum, et al., 2017) including post-clustering curation with the LULU algorithm (Frøslev et al., 2017) in order to obtain reliable alpha diversity estimates for OTU data.Although, it is widely acknowledged (e.g., Bálint et al., 2016) that species richness is difficult to estimate from sequencing data of environmental DNA, Frøslev et al. (2017) showed that careful bioinformatics processing can produce richness measures based on OTU data with strong correlation to richness metrics based on survey data.For this study, a simple OTU count was used as a DNA-based richness metric, after ensuring that variation in sequencing depth between samples only had a minor impact (results not shown).

| Conservation Index
A metric of conservation value was produced to test if plants are surrogates of species of conservation concern.For vascular plants, macrofungi, lichens, gastropods, spiders, and arthropods we used red list for Denmark (Wind & Pihl, 2004).For taxonomic groups lacking a national red list (bryophytes and plant-galling arthropods), an expert-based red listing was created for this project using the same criteria as the official red lists (bryophyte expert: Irina Goldberg, plant-galling arthropod expert: Hans Henrik Bruun).Each red listed species contributed to a weighted score of threatened species per site (the Conservation Index) as follows: red list status RE (regionally extinct) and CR (critically endangered) = 4 points, red list status EN (moderately endangered) = 3 points, red list status VU (vulnerable) = 2 points, and red list status NT (near threatened) and DD (data deficient) = 1 point.

| Abiotic factors and environmental calibration
We used field-measured abiotic variables to validate the plant-based environmental calibration.Environmental recordings and estimates included soil pH, soil C, N and P, soil moisture, leaf N and P concentrations, air temperature, light intensity, soil surface temperature, and humidity, number of trees >40 dbh, dead wood volume, and vapour pressure deficit (VPD).For further details on methods for collection of the abiotic data see Brunbjerg, Bruun, Brøndum, et al. (2017).
Mean Ellenberg Indicator Values were calculated for light conditions, soil nutrient status, soil pH, and soil moisture based on the plant lists for each site and the species' abiotic optima (Ellenberg et al., 1991).

| Natural habitat index
To supplement the plant-based environmental calibration, we calculated a plant-based natural habitat index reflecting land use intensity using 115,071 vegetation quadrats from the national monitoring of terrestrial biodiversity (Nielsen et al., 2012;Svendsen, Bijl, Boutrup, & Norup, 2005).Vegetation data were grouped in quadrats from Annex 1 habitats (A1) (EU Habitats Directive 1992) of conservation value (excluding nitrophilous tall herb fringes), other quadrats sampled in natural areas (Na), and quadrats sampled in agricultural (Ag) landscapes (road verges, hedges, soil banks etc.).For each species in the dataset, a natural habitat score was calculated as: where f () = frequency of species in the mentioned habitat category.
The species level score is thus a number between 0 and 1, where 0 implies that the species only occurs in agricultural biotopes and 1 implies that the species only occurs in habitats of conservation concern.
The natural habitat index was calculated for each site as the mean of species scores and reflects land use intensity and land use history.

| Analyses
We used Spearman rank correlation to test for correlation between species richness of vascular plants and the richness of other taxonomic groups including macrofungi, lichens, bryophytes, gastropods, plant-galling arthropods, carabid beetles, hoverflies, spiders, summed richness (all non-vascular plant species), fungal OTUs, eukaryote OTUs and Malaise OTUs.To validate the plant-based environmental calibration, Spearman correlations were calculated between Ellenberg Indicator Values, the natural habitat index, and measured environmental variables (soil moisture, soil C, N, and P, leaf N, and P, soil pH, surface, and air temperature, light intensity, number of trees >40 dbh, dead wood volume, and VPD).
We also grouped the 130 study sites into five different land-use intensity categories from protected Annex 1 habitats, over other uncultivated areas, plantation forest, and extensively farmed habitats to intensive farmland.AnovA followed by Tukey's post hoc tests was used to test for differences in mean natural habitat index value between the five habitat types.
To assess the efficiency of plants as predictors for multi-taxon richness, we divided the statistical analyses in a bivariate regression testing the predictive value of plant richness without taking other environmental co-variates into account and a multiple regression testing the effect of plant richness given environment, and estimating the total predictive potential of vascular plants (combining plant richness and bioindication).We produced models for a range of diversity metrics including species richness of above-mentioned groups as well as the Conservation Index.
Data exploration was carried out following the protocol described in Zuur, Ieno, and Elphick (2010).Collinearity was assessed using variance inflation factors (VIF) sequentially disregarding variables showing VIF values >3 from the VIF calculations (Zuur et al., 2010).Ellenberg nutrient status was found to be correlated with Ellenberg pH and our Conservation Index (VIF > 3) and was excluded from all models.If generalized additive model (GAM) smoothers fitted to the residuals of the models were conservatively significant (p < 0.01) for any of the predictors, we included polynomials in the final model.To account for the geographically nested design, we used generalized linear mixed modelling (GLMM) with a Poisson distribution and a log link function and with cluster (n = 15) as random intercept.If overdispersion was detected, we used GLMM with Negative Binomial error distribution instead of Poisson error distribution (Hilbe, 2011), given that the deviance information criterion (DIC) indicated a better fit based on a ΔDIC < 2, criterion.Model assumptions were verified by plotting residuals versus fitted values and versus each covariate in the model.We assessed the residuals for spatial dependency.Modelling was performed in r version 3.4.3(R Core Team, 2017), models were fitted by approximate Bayesian inference using the inlA package (Rue, Martino, & Chopin, 2009).
Explanatory variables were scaled prior to model implementation.
Marginal posterior distributions were summarised by 95% Bayesian credible intervals (BCI) corresponding to the 0.025 and 0.975 quantiles of the posterior distribution (Zuur, Ieno, & Seveliev, 2017).We chose not to perform model selection, but covariates can be considered statistically important if the 95% BCI do not overlap zero (Zuur et al., 2017).As an aid in the interpretation and comparison of model-based predictions, we calculated Pearson's product moment correlations between fitted values (only for fixed variables) and observed species richness of the response variables (pseudo R 2 ).

| Plant-derived environmental bioindication
Plant-derived environmental bioindication was supported by correlation to independent data.Community mean Ellenberg Indicator Values correlated well with corresponding measured abiotic factors (Table S1).The natural habitat index scored highest for EU Habitats Directive Annex 1 habitats, intermediate for other uncultivated habitats, plantation forests and extensively farmed sites, and lowest for sites with intensive cropping (Figure S1).The natural habitat index correlated negatively with Ellenberg nutrient status, indicating that plants with affinity to natural habitats generally occur in infertile environments (Figure S2).Plant species richness correlated positively with plant-derived Ellenberg soil pH and soil nutrient status and with measured soil pH, C, N, and P. In contrast, there was a negative correlation between plant species richness and the number of trees >40 dbh, dead wood volume, and canopy height (Table S1).

| Plant species richness as surrogate for other taxonomic groups
Spearman rank correlation between plant richness and species richness for other taxa revealed no significant correlation for macrofungi, bryophytes, lichens, carabid beetles, and summed richness.The richness of plant-galling arthropods, gastropods, spiders, and hoverflies showed significant positive relationships with plant richness (Table S2).

| Multiple regression based on plant-derived bioindication
Plant species richness was important with positive effects for species richness of all surveyed taxonomic groups except carabid beetles (Table 1 and Figure 2).Plant species richness was also important and positive for richness of fungal OTUs, Malaise OTUs and for the Conservation Index (Table 1 and Figure 2 F I G U R E 2 Relationships between plant species richness and species richness of other taxonomic groups, Conservation Index, and soil fungal, eukaryote, and Malaise operational taxonomic unit (OTU) richness in multiple and bivariate regressions and their 95% Bayesian credible intervals (BCI).Stippled versus full line indicates significance and non-significance at the 0.05 level (parameter estimates whose 95% BCI did not overlap zero), respectively.For model details see Table1 (Figure 3).The corresponding bivariate regression between plant richness and other richness metrics explained below 5% of variation in total richness for gastropods, total richness, bryophytes, fungi and eukaryote OTU richness, 5%-10% explained variation for hoverflies, spiders and Conservation Index and 10%-16% for fungal OTU richness, malaise OTU richness, and plant-galling arthropods (Figure 3).

Ellenberg light was generally important with positive effects
for flying insects such as hoverflies and Malaise OTU richness and otherwise negative or neutral effects on species richness and Conservation Index.Increasing Ellenberg moisture seemed to promote total richness, the richness of spiders, bryophytes, gastropods, lichens, hoverflies, and eukaryote OTU richness, whereas the effect on fungal OTU richness was weak and negative.Ellenberg pH had negative effects on spiders, bryophytes, hoverflies, and total richness, and positive effects on gastropods, and eukaryote OTU richness.The natural habitat index had a positive effect on bryophyte and lichen species richness and the Conservation Index, a negative effect on hoverflies, carabids and spiders and a unimodal effect on macrofungi species richness (Table 1 and Figure S3).

| D I S C U S S I O N
Terrestrial biodiversity of heterotrophic organisms relies on the build-up and diversification of organic matter produced by plants.
Therefore, it may seem reasonable to use the species richness of vascular plants as a proxy of total biodiversity.However, neither plants nor other individual taxa have hitherto been confirmed as reliable surrogates for other taxa.This point was supported by our simple cross-taxon correlations.Although plant species richness did correlate positively with four out of eight surveyed taxonomic groups, the overall performance of plant species richness per se as biodiversity surrogate was poor.A multivariate modelling approach, with both plant-derived bioindication and plant species richness, showed a much stronger-and consistently positive-effect of plant species richness on the species richness of other taxa, as well as on the Conservation Index and OTU richness measures, except soil eukaryote OTU richness, Malaise OTU richness, and carabid species richness (Figure 2).
Monitoring programmes often use plants as general indicators of conservation status of habitats without empirical justification.While plant species richness in itself may be a poor indicator for the richness of other species groups, plant indication may be a cost-effective approach to estimate environmental conditions (e.g., Diekmann, 2003) and possibly also the habitat quality (Andersen, Nygaard, Fredshavn, & Ejrnaes, 2013).We used Ellenberg Indicator Values for light conditions, soil moisture, nutrient conditions, and soil pH (Ellenberg et al., 1991), which are available only for the Central European flora.Our approach may still be applicable in other parts of the world because species scores from ordination of large and representative vegetation datasets typically reflect major environmental gradients (e.g., Ejrnaes, Aude, Nygaard, & Münier, 2002) and may replace Ellenberg Indicator Values in much the same way.While bioindication of environmental conditions is well developed, there is currently no standard approach to estimation of habitat quality by plant lists, despite scientific evidence that plants reflect land-use intensity and landuse history (Hermy & Verheyen, 2007).Plant-based habitat quality scores may be obtained e.g., by expert judgment (Kowarik, 1990) or by empirical evidence (Ejrnaes, Liira, Poulsen, & Nygaard, 2008;Ejrnaes et al., 2002).In this study, we calculated an index reflecting the affinity of plants to protected habitats and found a strong and negative correlation with Ellenberg nutrient status.
An obvious challenge for selecting surrogates that reflect general biodiversity is that different taxonomic or functional groups respond differently to habitat conditions.For example, lichens and bryophytes growing under extremely infertile conditions, often directly on stone and trees, will show a markedly different richness response along a fertility gradient than more competitive vascular plants.Likewise, generalist predatory beetles and spiders may be expected to respond differently than specialist herbivores such as gall wasps or aphids.Therefore, finding that plant species richness, after accounting for the abiotic environment, had a positive effect on species richness of all other taxa except carabids, lends strong support for the idea of biodiversity surrogacy and for vascular plants as useful surrogates.
The taxonomic groups used in this study varied strongly in their ecological dependence on plants.Plant-galling arthropods depend directly on specific plant species as hosts and represent a megadiverse group of phytophagous insects and mites with pronounced host specificity (Jaenike, 1990).Hoverflies are generally less hostdependent as larvae, but utilise plants as sources of pollen and nectar in the adult life stage.On the other hand, we would expect generalist herbivores or predators such as gastropods, carabid beetles and spiders, as well as primary producers such as lichens and bryophytes, to be causally unrelated to specific plant species.
Still, such species may respond to environmental conditions also influencing plant species richness.Macrofungi constitute several functional groups including both generalist decomposers and mycorrhizal symbionts, some of which are specialised on a single plant genus or species.Many decomposer fungi are also specific to certain plant genera or species, while some are necrophagous and highly specialised on arthropods or other fungi.Despite the difference in plant species specificity, the positive effect of plant species richness was consistent across taxonomic groups, pointing to a general applicability of plants as surrogates, even for predatory and decomposer organisms.Plant richness and environmental calibration obtained through bioindication together could account for 48% of the variation in richness of all other surveyed taxa combined.The figures for predicted OTU richness were also supportive with 24%-30% of variation explained.
The amount of explained variation was lowest for species richness of carabid beetles (12%), lichens (18%), and spiders (24%), possibly indicating that vegetation structure or microclimatic properties unrelated to plant community composition may be more important to species in these groups.Mobile generalist predators such as spiders and carabid beetles may rely less on site conditions than sessile species such as plants and fungi.A large proportion of lichens are epilithic or epiphytic on boulders and trees and therefore partly uncoupled from the prevailing environmental site conditions as reflected by vascular plants.
Despite the general usefulness of plants as surrogates, the amount of unexplained variation for specific groups such as lichens, carabids and hoverflies demonstrates that surrogates and indicators should be selected with due reference to spatial scale and the ecology of the target species groups (Kwok et al., 2011;Zurlini & Girardin, 2008).
In order to test the generality of vascular plants as surrogates, we also included three richness metrics derived from DNA metabarcoding-soil fungal OTUs and soil eukaryote OTUs from eDNA and aerial arthropod OTUs from Malaise trap DNA.Despite a thorough sample of 81 regularly spaced soil cores, we have merely covered an approximate 0.01% of the soil surface of each study site.This could pose a bottleneck for getting a representative sample of OTUs in diverse and heterogeneous habitats.We assume that the eukaryotic and the fungal genetic markers are targeting a soil community depending on micro-climate and soil composition, and less on vegetation-at least compared to the organisms recorded above the ground.Furthermore, the inherent problems in getting reliable richness estimates from eDNA sequencing are widely acknowledged (e.g., Bálint et al., 2016).
Nevertheless, the positive effect of plant species richness was reproduced for OTU-richness, albeit insignificant for eukaryotes after taking account of environment, and the explained variance by multiple regressions with plant-derived environmental variables approached 25% for three taxonomically very different OTU taxa.
Rare and threatened species are particularly important to conservation, and we demonstrated that a plant-based model could explain 23% of the variation in our Conservation Index based on occurrence of red-listed species.Our sites were only 40 × 40 m and, thus, too small for a representative sampling of very rare species.
With larger plots, we would expect a higher proportion of explained variation.A general index of site uniqueness could replace the use of rare species for assessment of conservation value of such small sites.
Our natural habitat index was the strongest predictor of variation in the Conservation Index which is in accordance with evidence for the preferences of threatened species for rare natural habitats (Pearman & Weber, 2007 and references therein).
We find it encouraging that the richness of vascular plants is a consistent positive predictor of multiple functional groups comprised by our multi-taxon species richness estimate.However, looking at the direct trophic effects, we see opportunities for further improvement of plant surrogacy.It has long been acknowledged that plants serve as mutualistic partners for other organisms (e.g., Elton, 1949).With respect to the diversification of organic matter, Southwood (1961) and later work by Brändle and Brandl (2001) quantified the richness of phytophagous insects on European trees and showed that the size of their associated biotas vary enormously and predictably, i.e., large, long-lived and omnipresent species may harbour a more diverse pool of insects than small annuals or uncommon species.A thorough examination of reported interactions between plants and associated invertebrates and fungi may be used to create a more powerful surrogate for total biodiversity than the mere number of plant species.
Vascular plants play an important role in conservation management and monitoring.In this study, we demonstrate that plant species are useful surrogates for biodiversity at large, but only when environmental bioindication is taken into account.In line with the ecospace framework for biodiversity (Brunbjerg, Bruun, Moeslund, et al., 2017), future research into the diversification of organic matter (dung, dead wood, litter) and differentiation between vascular plants with varying importance as hosts may further improve the value of plant-related indicators as surrogates of biodiversity in general.

AC K N OW L E D G E M E N T S
R.E., L.B., I.G., T.L., T.G.F., C.F., and A.K.B. were supported by a grant from VILLUM foundation (Biowide, VKR-023343).We thank Aimee Classen and Greg Newman for assistance with the analysis of soil properties, Vagn Alstrup ( †), Ulrik Søchting and Roar Skovlund Poulsen for lichen surveying, Karl-Henrik Larsson for aid in identifying critical corticioid fungi, Leif Örstadius for identifying Psathyrella collections.

F
I G U R E 1 Map of Denmark showing the location of the 130 sites grouped into 15 clusters within five regions (Njut: Northern Jutland; Wjut: Western Jutland; Ejut: Eastern Jutland; FLM: Funen, Lolland, Møn; Zeal: Zealand) ).Multiple regression of species richness of the selected taxa varied in per cent explained variation by fixed variables from 12% for carabid beetle richness to 55% for gastropod species richness TA B L E 1 Model results for the full GLMM models (Poisson or negative binomial) Model results show the effects of plant species richness (plant_rich), Ellenberg light (E_light), Ellenberg moisture (E_moisture), Ellenberg pH (E_pH) and natural habitat index (nat_index, nat_index 2 ), on species richness of various taxonomic groups (total, bryophytes, carabid beetles (carabids), plant-galling arthropods (gallers), gastropods, hoverflies, lichens, macrofungi, and spiders), Conservation index, and OTU richness (eukaryote, fungal, and Malaise).Intercept, parameter estimates (marginal posterior means), and 95% Bayesian credible intervals (BCI, i.e., the 0.025 and 0.975 quantiles of the posterior distribution) are given in parentheses.Parameter estimates with 95% BCI not overlapping zero are shown in bold.GLMM: generalized linear mixed modelling; OTU: operational taxonomic unit.50 Barplot of percentage variance explained by the multiple and bivariate regressions of species richness of the taxonomic groups, Conservation index, and soil fungal, eukaryote and Malaise operational taxonomic unit (OTU) richness.