Multi-taxon inventory reveals highly consistent biodiversity responses to ecospace variation

Amidst the global biodiversity crisis, identifying drivers of biodiversity variation remains a key challenge. Scientific consensus is limited to a few macroecological rules, such as species richness increasing with area, which provide limited guidance for conservation. In fact, few agreed ecological principles apply at the scale of sites or reserve management, partly because most community-level studies are restricted to single habitat types and species groups. We used the recently proposed ecospace framework and a comprehensive data set for aggregating environmental variation to predict multi-taxon diversity. We studied richness of plants, fungi, and arthropods in 130 sites representing the major terrestrial habitat types in Denmark. We found the abiotic environment (ecospace position) to be pivotal for the richness of primary producers (vascular plants, mosses, and lichens) and, more surprisingly, little support for ecospace continuity as a driver. A peak in richness at intermediate productivity adds new empirical evidence to a long-standing debate over biodiversity responses to productivity. Finally, we discovered a dominant and positive response of fungi and insect richness to organic matter accumulation and diversification (ecospace expansion). Two simple models of producer and consumer richness accounted for 77 % of the variation in multi-taxon species richness suggesting a significant potential for generalization beyond individual species responses. Our study widens the traditional conservation focus on vegetation and vertebrate populations unravelling the importance of diversification of carbon resources for diverse heterotrophs, such as fungi and insects.


51
For centuries, ecologists have struggled to understand and explain spatial and temporal variation in 52 biodiversity, with increasing societal attention motivated by the global biodiversity crisis. Most 53 models and theories of biodiversity refer to specific taxonomic groups and ecosystems (1-6), leaving 54 us with no general rules and ecological theories of diversity and without prioritization tools at the 55 local spatial scales of practical conservation planning and management. Several studies have 56 investigated the potential for using selected species groups to represent the wider conservation 57 interests, but based on a global meta-analysis, Westgate, Barton, Lane and Lindenmayer (7)  58 concluded that the data undermines the assumption that a taxonomic subset can represent the 59 wider biodiversity. The disappointing conclusion has been somewhat contradicted or modified lately 60 by studies showing promising potential for cross-taxon congruence in species composition and 61 compositional turnover (β-diversity) (8) and also for species richness, but only after accounting for 62 environmental variation (3,9). Along this line of reasoning, we set out to test the hypothesis that 63 terrestrial multi-taxon diversity can be predicted across contrasting environments without detailed 64 consideration of an intractable diversity of taxonomic groups, response groups or habitat types. We 65 are thus not testing the surrogacy hypothesis per se, but rather the idea that multi-taxon diversity 66 can be predicted from a low-dimensional ecological space, ignoring the possible multitude of 67 response shapes of the individual taxonomic groups or species. 68 We applied the recently proposed ecospace framework (2-4, 10, 11) for a formal and structured 69 quantification of environmental variation. Ecospace represents the total environment in space and 70 time, in which the individual (species) colonizes, grows, reproduces, and dies or goes extinct. 71 Ecospace may help reduce environmental complexity to a tractable number of dimensions and 72 measurable variables. We have proposed to subdivide ecospace into three components each 73 signifying important aspects of an area for its potential biota: 1) The abiotic environment (ecospace 74 position), 2) The accumulation and diversification of organic matter (ecospace expansion) and 3) The 75 spatio-temporal continuity (10). 76 The position of a site in n-dimensional environmental hyperspace (e.g., mean values of soil moisture, 77 pH, soil fertility, and temperature) is essential to sessile organisms like plants and soil-fungi unable 78 to move across their local environment (6), but even mobile animals respond to abiotic conditions 79 when they select their habitat (12). Expansion, i.e. the accumulation and diversification of organic 80 matter, is particularly important as an energy source for consumers (13), but may also provide 81 substrate for e.g. epiphytic plants and lichens (14). Expansion presupposes primary production and 82 subsequent differentiation into leaves, roots, stems, flowers, bark, wood, dung etc. On evolutionary 83 time scales, every differentiated pool of live or dead organic matter has provided opportunity for 84 heterotrophic niche differentiation and speciation (15). 85 The spatial continuity of habitats is expected to be particularly important for short-lived and poorly-86 dispersed species moving among patches varying in suitability over time, but less important for 87 species with long distance dispersal (16). Temporal continuity on the other hand should be 88 particularly important for organisms with limited dispersal ability and high persistence, such as some 89 plants and fungi, which are sensitive to changing habitat conditions (6). 90 The aim of this study was to investigate the hypothesis that multi-taxon α-diversity, including 91 'genetic richness' (5) from environmental DNA, can be treated as a general, predictable biotic 92 response to environmental variation represented by a low number of key factors. We assess this 93 using ecospace as a framework for guiding the study design, environmental mapping and data 94 analysis. 95

96
Individually, the nine species richness models explained between 27 % (flying insects) and 62 % 97 (decomposing fungi, herbivorous arthropods) of the variance in species richness across 98 environmental gradients (Fig. 2). When we pooled species into producers and consumers, the 99 explained variation based on cross-validated predictions was comparable to the individual 100 models (producers: 49 %, predictions producers: 50 %; consumers: 73 %, predictions consumers: 101 69 %). The model for species richness of all groups explained 54 % of the variation in biodiversity 102 across sites, considerably less than the predictions from the sum of all nine separate species 103 group models (74 %). However, the summed predictions from the models of producers and 104 consumers explained 77 % of total species richness. Based on this result, we focused further 105 reporting of results on models of producer and consumer biodiversity (Fig. 3). Producer richness 106 was primarily explained by position (Fig. 3a) S1). Presence of a shrub layer also 112 promoted producer richness. Finally, producer richness increased with temporal continuity. 113 Variance partitioning revealed that position explained most variation in producer richness with 114 minor contributions from expansion and continuity (Fig. 3a). Consumer richness increased with 115 presence of a shrub layer, higher flower abundance, and a high index of insect host plant 116 abundance (Fig. 3b). For position, consumer richness increased with air temperature and 117 decreased with incoming light. We found no effect of continuity on consumer richness. Variance 118 partitioning revealed that most variation in consumer richness could be explained by expansion 119 compared to a minor contribution from position (Fig. 3b). 120 The following expansion variables promoted species richness for the individual response groups: 121 litter mass for decomposing fungi, soil carbon content for arthropod detritivores, dung for 122 decomposing and symbiotic fungi as well as total species richness, dead wood for decomposing 123 fungi and floral abundance for total species richness, consumers, flying insects and herbivorous 124 arthropods ( Table 1). Richness of all consumer species groups increased with either plant species 125 richness or the indices of host plant availability for fungi or insects indicating a bottom-up effect 126 going from primary producer richness to consumer richness. The presence of a shrub layer had a 127 consistent and positive effect on the richness of most response groups. We found very few 128 significant effects of within-site heterogeneity on species richness indicating that these were of 129 little importance compared to the effects of between-site variability (Table 1). 130 In general, we found linear richness responses to underlying abiotic gradients, with the exception 131 of unimodal responses of predatory arthropods to pH, and bimodal responses of lichens, vascular 132 plants and producers to soil moisture. For soil fertility, we observed unimodal responses for 133 vascular plants, mosses, symbiotic fungi as well as the pooled groups of producers and total 134 species richness (Fig. 4). 135 Models for 'genetic richness' of soil and insect trap DNA were in the lower range of model 136 performance with cross-validated R 2 values of 23 % for eukaryotes, 27 % for fungi and 29 % for 137 flying insects. While fungi and flying insects confirmed the consistent positive response to plant 138 richness, none of the DNA groups responded to presence of a shrub layer (Table 1). Soil pH, 139 moisture, fertility and litter mass affected soil fungal and eukaryotic DNA richness in various ways 140 indicating the importance of the soil environment for its biota (Table 1). 141 We found no justification for including a random variable (Generalized Linear Mixed Model) to 142 account for spatial patterns in the data and only two of 15 models (the flying insects and DNA 143 flying insects models) showed significant spatial autocorrelation after modelling. Spatial signals 144 seem to be of minor importance in this study. 145

146
Across major terrestrial ecosystems within a region, as much as 77 % of the variation in multi-147 taxon richness was accounted for by our models. The remarkably high explanatory power arose 148 after grouping species into producers (autotrophic organisms) and consumers (heterotrophic 149 organisms). Studies assessing the surrogacy power of multiple environmental variables on 150 multiple taxonomic groups simultaneously, are rare. A meta-study found the surrogacy power of 151 environmental variables to be very poor compared to taxonomic surrogacy (17) while multi-152 metric site-conditions explained 48 % of variation in total species richness of a range of taxa 153 (ants, beetles, spiders, wasps, flies, butterflies, reptiles, birds, vascular plants, bryophytes and 154 lichens) in Australian woodlands (18). Our regression models for individual species groups show 155 considerable variation in the selection of variables and this could easily lead to the erroneous 156 conclusion that each taxonomic group needs individual consideration. Nevertheless, the 157 predictions derived from models dedicated to taxonomically and ecologically more specialized 158 species groups explained less variation in α-diversity than the two models based on high level 159 aggregation. This level of generalization is striking considering the contrasting life history traits 160 and modes of resource acquisition of the species aggregated as producers or consumers, and it 161 challenges commonly expressed concerns that α-diversity is necessarily contingent on taxonomy 162 and ecology (e.g., 19, 20, 21). 163 On the other hand, we found a large drop in explained variation when we tried to aggregate 164 producers and consumers and model all species in one model. This result emphasizes the 165 fundamental difference between sessile autotrophic organisms, whose resource acquisition is 166 largely controlled by abiotic conditions, and heterotrophic organisms, whose resource acquisition 167 relies on biotic resources. The important split between heterotrophic and autotrophic organisms 168 is underpinned by the contrasting role of ecospace position and expansion in the producer and 169 consumer models. Further, plant richness is an important predictor of consumer richness, but to 170 avoid circularity plant richness is excluded from the total richness model. 171 Based on island biogeography (14) and metapopulation theory (15) Table 1) corroborates a similar species pool effect for consumers (24). Despite 182 recent advances in the translation of eDNA data into diversity metrics (25, 26) the relatively poor 183 performance of models for DNA groups might point to remaining metagenomic challenges (25, 184 27). It is also possible however that the variation in below-ground and above-ground biodiversity 185 is determined by different factors and that ecospace expansion in particular may need to be 186 refined to include and differentiate between organic matter pools that support below-ground 187 biodiversity (28). 188 The most consistent predictors of consumer richness were the presence of a shrub layer and high 189 vascular plant richness. Specialized organic matter, such as dead wood, dung and flowers were 190 found to be important to fungal and insect richness (29-31), judged from their representation in 191 the detailed species group models. In this way, our study reveals a bottom up regulation of five regions across Denmark (Fig. 1). We allocated 100 sites to span the most important natural 212 gradients affecting biodiversity in Denmark i.e. gradients in soil fertility, soil moisture and 213 successional stage: from nutrient rich to nutrient poor, from dry over moist to wet and from 214 open, closed and forested vegetation. 90 of these sites were selected randomly within 18 215 predefined strata, whereas 10 sites were selected by the amateur natural historian community to 216 represent biodiversity hotspots for different species groups. The remaining 30 sites were sampled 217 randomly from six strata cultivated for production: plantations (beech, oak or spruce) and fields 218 (rotational, grass leys or set aside). Randomization was achieved by selection of sampling areas 219 and potential sites from the office desk based on geo-referenced information -for some rare 220 strata we also consulted local experts. To minimize spatial autocorrelation, the minimum 221 distance among sites was 500 m with a mean nearest distance among sites of 2291 m. Within 222 each site, we sampled vascular plants, mosses lichens, fungi, and arthropods, and we included 223 metabarcoding of DNA from soil samples and insect traps to reflect the 'genetic richness' (i.e. 224 number of operational taxonomic units -OTUs (5) of eukaryotes, fungi and flying insects. 225 Furthermore, we collected data reflecting ecospace i.e. abiotic position, biotic expansion and 226 spatio-temporal continuity (described below). For further details on site selection and data 227 collection, see Brunbjerg,et al. (11). All field work and sampling was conducted in accordance 228 with Responsible Research at Aarhus University and Danish law.  (Table S2). photographs and additional information from land mapping of woodlands, fields, grassland, 331 heathland, meadows, salt marshes and mires. The four buffer sizes were similar and highly 332 correlated. The 500 m buffer was used for analyses as most of the studied species were expected 333 to have relatively limited dispersal and small area requirements. 334 Temporal continuity: Temporal continuity was estimated by time since major land use change 335 within the 40 m × 40 m site. For each site, a temporal sequence of aerial photos and historical 336 maps was inspected starting with the most recent photos (photos from 2014, 2012, 2010, 2008, 337 2006, 2004, 2002, 1999, 1995, 1968, 1954, 1945)

Response variables: 363
We divided species into response groups according to taxonomy (plants, mosses), trophic level 364 (macrofungi) and trophic level and mobility (invertebrates). Grouping is complicated for insects 365 because the biology and mobility may depend on live stage. Hoverflies for example have larval 366 stages spanning from detritivores over predators to galling herbivores, whereas the adult flies 367 mainly feed on flowers. An entirely trophic categorization is further intractable given that many 368 resource strategies in insect larvae are unknown. Moreover, the species richness response may 369 rely on the mobility and preferences of the observed imago rather than the occupation of the 370 larvae. We therefore decided to follow a pragmatic division based on biological reasoning. 371 The response groups of our study included vascular plants, mosses, lichens, decomposing fungi, 372 symbiotic fungi, flying insects (highly mobile insects dependent on ephemeral food sources such 373 as dung, flowers, fungi and dead wood), herbivores (mobile insects dependent on live, sessile 374 plants), detritivores (invertebrates dependent on dead carbon sources), predators (arthropods 375 dependent on live animals). For details on grouping, see Table S3. In order to investigate the 376 potential for generalization across species groups we pooled species into the larger groups of 377 producers (vascular plants, mosses and lichens), consumers (fungi and arthropods) and total 378 (producers and consumers). In addition, we used richness of OTUs (operational taxonomic units 379 (5)) of eukaryotes and fungi from soil eDNA and arthropods from eDNA extracted from ethanol 380 from Malaise traps. For details on collecting species and eDNA data see Brunbjerg,et al. (11). 381 eDNA datasets: The preparation of the fungal and eukaryote eDNA datasets have been published 382 in Frøslev,et al. (44) and Fløjgaard,et al. (45) respectively. The insect DNA dataset was produced 383 by extracting DNA from the ethanol from the bulk insect Malaise traps and metabarcoding with 384 an insect specific 16S primer. 45 ml ethanol and 1.5 ml of 3M sodium acetate were added to a 50 385 ml centrifuge tube, and left in a freezer for DNA precipitation overnight, then centrifuged for 40 386 minutes. The dried pellet was extracted with the Qiagen DNeasy blood and tissue kit (Qiagen,387 Germany) with minor modifications. The extracted DNA was normalized, amplified, sequenced 388 and analyzed according to the overall procedures described for 16S insects in Fløjgaard,et al. 389 (45), but including curation of OTUs with LULU (26) and taxonomic assignment with a custom 390 script, and exclusion of non-arthropod sequences. Data from the two different collecting events 391 were handled separately and the sequences were then combined for each site. Bioinformatic 392 processing including links to all data is documented at 393 https://github.com/tobiasgf/biowide_synthesis. 394

Explanatory variables and statistical analyses: 395
We built generalized linear models (GLMs) to predict species richness of selected response 396 groups based on the best selection of ecospace variables. In addition, we built GLMs to predict 397 the summed richness of vascular plants, mosses and lichens (producers) and fungi and 398 invertebrates (consumers). Finally, we built an overall species richness model predicting the total 399 richness of all observed species. 400 For each model we made a preliminary screening and selection of relevant variables, excluding 401 variables with no hypothesized relationship to the species group or variables dependent on the 402 response variable. We further constrained the response direction and shape to ecologically 403 plausible responses (46, 47) implying an exclusion of negative effects of expansion, continuity 404 and heterogeneity variables on species richness -more resources, more diverse resources, more 405 environmental variation and increasing temporal and spatial continuity are all hypothesized to 406 increase species richness if anything. This decision is justified by the large number of variables in 407 our study increasing the risk of including spurious correlations in the models and thereby 408 covering important causal relationships (Table S4) (leave-one-out cross validation) for these models while applying the ΔAIC < 2 rule. 426 Model performance of the final models were evaluated using leave-one-out cross-validation. To 427 evaluate the effect of generalizing we compared the performance of models on aggregated 428 species groups with the performance of the corresponding models on the individual species 429 groups. in order to do this pearson correlations between the sum of predictions from specific 430 group models and the sum of observed species of producers, consumers and total (referred to as 431 the summed predictions from nine species group models and summed predictions for producers 432 and consumers), were calculated respectively. 433 Variation partitioning was calculated on final models for each component of ecospace (position, 434 expansion, continuity and co-variables) as follows: 435 Data exploration was applied following Zuur, Ieno and Elphick (47)  richness. For example, the hatched dark green bar is the correlation between the pooled 606 predictions from the plant, moss and lichen models respectively and the observed richness for all of 607 these groups. This compares to the solid dark green bar representing the correlation between the 608 predictions from a model including all producer groups and the observed values. For all groups, the 609 first hatched bar represents the correlation for the summed predictions from all nine species 610 group models, whereas the second hatched bar represents the correlation for the summed 611 predictions for the producer and the consumer models. 612 soil moisture (SMI) and pH and the species richness of plants, mosses, lichens, producers, 623 symbiotic fungi, predatory arthropods, total (producers+consumers), and eDNA eukaryotes in the 624 respective multiple regression models. 625