No single driver of biodiversity: divergent responses of multiple taxa across land use types

Abstract

Understanding the responses of biodiversity to different land use regimes is critical for managing biodiversity in the face of future land use change. However, there is still significant uncertainty around how consistent the responses of different taxonomic groups to land use change are. Here, we use a combination of high‐throughput environmental DNA sequencing and traditional field‐based survey methods to examine how patterns of richness and community composition correlate among four domains/kingdoms (bacteria, fungi, plants, and metazoans) and the four most‐abundant animal taxonomic groups (arachnids, Collembola, insects, and nematodes) across five different land use types (natural forest, planted forest, unimproved grassland, improved grassland, and vineyards). Richness for each taxonomic group varied between land use types, yet different taxa showed inconsistent responses to land use, and their richness was rarely correlated. This contrasted with community composition of taxonomic groups, for which there was relatively good discrimination of land use types and there was strong correlation between group responses. We found little evidence for consistent drivers of taxonomic richness, yet identified several significant drivers of community composition that were shared across many groups. Drivers of composition were not the same as the drivers of diversity, suggesting diversity and composition are independently controlled. While land use intensification has been viewed as having generally negative effects on biodiversity, our results provide evidence that different taxa respond divergently across different land uses. Further, our study demonstrates the power of high‐throughput sequencing of environmental DNA as a tool for addressing broad ecological patterns relating to landscape biodiversity.

Introduction

As the human population continues to grow rapidly, the pressure for development of agricultural land also increases (Meyer and Turner , Matson and Vitousek ). With ongoing declines in biodiversity globally (Butchart et al. ), understanding the responses of biodiversity to different land use regimes is critical for developing management strategies to better conserve biodiversity in the face of future land use change. Land use intensification, that is, a shift from extensive (low input) to intensive (high input) land uses, is often seen as having broadscale negative impacts on local biodiversity (Newbold et al. ). For example, significant changes in land use (e.g., replacement of primary forest with crops) and intensification within land uses have been associated with declines in plant and animal diversity (e.g., Ponge et al. , Schulze et al. , Attwood et al. ), animal abundance and richness (e.g., Donald et al. ), reduced microbial beta diversity (Guan et al. ), and reductions in the diversity of other soil organisms (Minor and Cianciolo ). Other forms of land management, such as adoption of biodiversity‐friendly management (Scherr and McNeely ), decreased intensification (Carswell et al. ), or conservation protection (Bruner et al. ), can have opposite effects on biodiversity. However, exceptions to these trends exist (e.g., dipteran flies responded positively to land use intensification in Allen et al. ), and the response of many taxa is poorly known. Current understanding of the responses of a broad range of taxonomic groups to different land use regimes is therefore limited, yet, given the widespread modification of terrestrial ecosystems by humans (Hooke et al. ), is crucial to reversing current biodiversity declines.

One key uncertainty is whether the effects of land use on the basic components of biodiversity (i.e., the richness and composition of different taxonomic groups) are correlated. Most studies have focused on the response of a few specific taxa and/or a narrow range of land use types, limiting our ability to infer regional responses of biodiversity to land use. For example, Gillison et al. () reported correlated changes in plant and termite species richness along a gradient of land use intensification, and Mueller et al. () found correlated changes in plant and fungal community composition in response to forest‐to‐pasture land use change. However, other studies have found no correlation, for example, between species richness of birds, butterflies, and ants along an intensification gradient within coffee plantations (Perfecto et al. ) or communities of flies, beetles, and mammals with conversion of land for cropping (Burel et al. ). Relatively few studies have examined biodiversity responses across multiple taxa and multiple land uses. In one of the most comprehensive studies to date, Allen et al. () found an overall decline in richness across 49 fungi, bacteria, plant, vertebrate, and invertebrate taxa within managed grasslands of increasing intensification. Using the same study system, Manning et al. () found that the richness responses of many of these taxa, but not all, were positively correlated and that the strength and number of significant correlations varied with land use intensity. However, we found no prior studies that examined composition and diversity responses of multiple high‐level taxa across a wider range of land uses. As such, it remains uncertain as to whether these results, from a single land use type (grasslands), can be generalized to landscapes.

The mechanisms underpinning the response of different taxonomic groups to land use are likely to determine whether or not taxa have correlated responses. Two of the main factors associated with different land uses are changes in the community composition and diversity of plants. As the plant community forms the basis of both above‐ and belowground food webs, we might expect to see coordinated changes in composition and diversity of multiple taxa across land uses if taxa are primarily bottom‐up (resource) controlled, if tight linkages among species are common (Hooper et al. , Wardle ), or if higher trophic levels respond to habitat heterogeneity created by increased plant species richness (Stein et al. ). Indeed, plant species richness and community composition have been linked to the richness and community composition of higher trophic levels (e.g., Qian and Ricklefs , Scherber et al. ). Similar responses to abiotic environmental drivers associated with a given land use may also result in correlated responses (Wolters et al. ). For example, different land uses are often associated with different geomorphology (e.g., parent material and slope), as this influences how suitable land is for a given use, and with different soil properties (e.g., fertility) due to variation in both geomorphology and management practices (Fig. ). Such factors have been linked to changes in the composition and diversity of above‐ and belowground taxa, although effects are often taxon‐specific (e.g., Mueller et al. ). The relative importance of abiotic and biotic drivers for different taxa is likely to determine whether or not they respond similarly to land use (Fig. ). However, most studies examining links between multiple taxa across land use types have not assessed whether results can be explained by differences or similarities in underlying drivers.

DNA metabarcoding, from environmental or multi‐organism samples, offers significant promise for characterizing biological communities and studying biodiversity patterns (Cristescu ). Here, we use this approach, together with traditional field‐based survey methods, to examine how patterns of community composition and richness correlate between four high‐level taxa (domains and kingdoms: bacteria, fungi, Metazoa, and plants) and the four most‐abundant mid‐level animal groups (classes and orders: Collembola, arachnids, Chromadorea nematodes, and insects) across five different land use types (natural forest, planted forest, unimproved and improved grasslands, and vineyard). Using this multi‐taxonomic‐level and multiple land use approach, we test the following hypotheses: (1) that land use type significantly affects the richness and community composition of different taxonomic groups; (2) that richness and community composition of different taxonomic groups are correlated; and (3) that significant correlations among richness and community composition of different taxa will be underpinned by each group having similar environmental drivers, and the absence of such correlations will be linked to each group having divergent drivers.

Methods

Study sites and sampling

We selected 30 study sites positioned on a 4 × 4 km grid within the Wairau River Catchment, Marlborough, New Zealand (Appendix S1: Fig. S1). The sites were <1000 m altitude and included six plots from each of natural forest (e.g., dominated by ectomycorrhizal Nothofagaceae with arbuscular–mycorrhizal hardwoods and ferns), planted forest composed of Monterey pine (Pinus radiata) or Douglas–fir (Pseudotsuga menziesii), unimproved semi‐natural grassland (typically sheep and beef cattle grazing, including minor shrubs and small trees such as ectomycorrhizal Myrtaceae), improved grasslands of European origin (intensively managed, usually ryegrass or clover dominated, with irrigation, cultivation, and intensive fertilization), and vineyards (Appendix S1: Table S1). Further details of site selection were provided by Orwin et al. () and in Appendix S1. At each site, we established a 20 × 20 m plot and performed a survey of all vascular plants following standard protocols (Hurst and Allen ). This included recording the identity and foliar cover of all vascular species within each plot using an ordinal scale in five height classes (Appendix S1). As a measure of plant species abundance, we calculated a single site‐level percentage cover score for each species by converting the ordinal scores to percentage cover using the geometric mean for each cover class, averaged across all height classes and plots within a site. We subdivided each plot into 16 subplots, each 5 × 5 m, and took 15 cm deep soil cores using a 4.75 cm diameter corer (AMS, Idaho, USA) at regular intervals along the edges of each subplot (24 cores per plot; Appendix S1: Fig. S2). We removed the surface litter layer prior to coring to maximize comparability of samples where litter was common (e.g., natural forest) with those where litter was generally sparse (e.g., improved grassland). We pooled the soil cores and stored at 4°C until analysis.

Subsampling and DNA extraction

We isolated genomic DNA from eight different fractions of the pooled soil cores from each plot (two bulk soil samples, roots, and five invertebrate concentrations). We used two different approaches to isolate DNA from bulk soil. First, we isolated DNA from up to 5 g of bulk soil using the PowerMax DNA Isolation Kit (MO BIO Laboratories, Carlsbad, California, USA). Second, we isolated DNA from 15 mL of phosphate buffer (Taberlet et al. ), in which we had saturated and shaken 250 g of soil for 25–30 min, using the PowerMax DNA Isolation Kit from step 7 (solution C2) in the manufacturer's protocol. We isolated DNA from a subsample of milled roots using a modified protocol for the PowerSoil DNA Isolation Kit (MO BIO Laboratories). Due to the expected low population densities of some invertebrate groups, we concentrated five different invertebrate fractions from the total soil volume: macroinvertebrates, nematodes (gravity and centrifuge concentration methods), and mites (heat lamp and heptane flotation concentration methods; details of methods are provided in Appendix S1). We isolated DNA from each of these using specific lysis buffers and modified protocols for the DNeasy Blood & Tissue Kit (Qiagen, Valencia, California, USA). Full details of the DNA extraction methods are provided in Appendix S1.

PCR and DNA sequencing

Details of the fusion primers (incorporating 454 Lib‐L type adapters) that we used to amplify DNA are provided in Appendix S1: Table S2. Full details of the combinations of DNA extracts and primers used, and the thermocycling conditions for each primer pair are provided in Appendix S1. PCRs were performed in final volumes of 25 μL and included 2 μg/mL RSA, 1 × PCR buffer, 1 mmol/L MgCl2, 80 μmol/L each dNTP, 0.4 μmol/L each primer, 1.25 U FastStart Taq (Roche Applied Science, Indianapolis, Indiana, USA), and 2 μL of DNA template.

We visualized amplicons on a 3.5% (for products <300 bp) or 2% (for products >300 bp) agarose gel and pooled those produced using the same primer (i.e., from different DNA extracts) for each site. Where primer artifacts were absent, we purified the pooled amplicons using a mix of exonuclease I (Exo) and shrimp alkaline phosphatase (SAP) at a ratio of 0.8U Exo: 1U SAP: 5–10 μL PCR product. This was incubated at 37°C for 45 min, 80°C for 15 min, and 15°C for 3 min. Where non‐target amplification or primer artifacts gave multiple bands, we purified the pooled amplicons using SPRIselect (Beckman Coulter, Indianapolis, Indiana, USA). For amplicon sizes <250 bp, we followed the right‐side selection protocol, and for amplicon sizes >250 bp, we followed the left‐side size selection protocol. Final elutions were performed with TE buffer. We quantified the purified PCR products using a Qubit 2.0 fluorometer (Life Technologies, Grand Island, New York, USA).

We performed a final pooling step (of all amplicons for each plot), accounting for the preferential bias toward short DNA fragments in emulsion PCR (see Appendix S1 for details). We used 3 μL of each of the 30 pooled amplicons for quality control purposes, including quantification of DNA and analysis using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, California, USA) and a GX DNA High Sensitivity LabChip (Caliper Life Sciences, Hopkinton, Massachusetts, USA). The pooled amplicons were submitted to Macrogen (Korea) where each was sequenced on 1/8 of a 454 GS‐FLX titanium plate.

Abiotic drivers

To assess the potential drivers of biotic communities, we used geomorphology/climate and soil property variables along with abundance of different plant species (Fig. ; Appendix S1: Table S1). Geomorphology/climate variables included parent material, mean annual rainfall, and mean annual temperature data and were obtained from Land Environments of New Zealand (Leathwick et al. ). Slope was measured in the field with a clinometer, and slope shape and elevation were taken from topographic maps. We measured several soil properties as potential drivers of changes in community composition and richness. Methods for total soil C and N, total P, Olsen P, and pH are given in Orwin et al. (). We also measured cation exchange capacity; base saturation; and Ca, Mg, K, and Na contents (using the ammonium acetate leaching method; Blakemore et al. ).

Bioinformatics and statistics

DNA sequences were clustered and linked to taxonomic identity using a bioinformatic pipeline based on Uparse (Edgar ). Sequences were trimmed to 250 bp length (except for nematode and invertebrates which were trimmed to 300 bp based on examination of quality scores), filtered to a maximum expected error of 1.0, dereplicated, and clustered into Operational taxonomic unit (OTUs) using Usearch, which removes putative chimeric sequences. Singletons were not used to form clusters, but could be matched to a cluster. To identify sequences, all OTUs were matched against a database of identified sequences (see Appendix S1) using BlastN.

Testing of hypothesis 1 required estimating OTU richness and determining community composition. OTU richness in next‐generation sequencing is highly sensitive to differences in numbers of sequences obtained across samples (Dickie ). We therefore used rarefaction analysis to the smallest sequencing depth within each primer to obtain a robust comparison of relative richness between samples using function “rarefy” in vegan (Oksanen et al. ). This function determines the number of species detected in random subsamples of identical sequence counts. In a few cases, plots with very low sequence counts were dropped from the analysis to avoid having very low rarefaction results for all plots. Plant richness was based on actual count of species observed. Land use effects on richness (and all subsequent analyses) were tested first at the level of domain/kingdom (bacteria, fungi, metazoans, and plants) and second for the four most dominant metazoan classes/orders (Collembola, arachnids, nematodes, and insects). These groups were present across all land use types. We focused on animals as a highly diverse group that has been demonstrated to lose biodiversity with land use intensification. The effects of land use on richness were tested for each taxonomic group using ANOVA.

Community ordinations were conducted using the metaMDS function of the vegan package in R (Oksanen et al. ) using Bray distances. As a measure of plant species abundance, we calculated a site‐level percentage cover score for each species by converting the ordinal scores to percentage cover using the geometric mean for each cover class, averaged across all height classes and plots within a site. For other taxa, we used the proportion of sequence reads as a measure of abundance, which although supported by some empirical evidence can also be biased by experimental factors (such as amplicon length and primers used; Amend et al. , Engelbrektson et al. ). Nonetheless, treating next‐generation sequencing data as simply presence–absence would over‐inflate the importance of rare sequences, which are more prone to include errors (Dickie , Lindahl et al. ). Separation of biotic communities by land use was tested using permutational multivariate analysis of variance as function “adonis” in the vegan package of R (Oksanen et al. ). To test for correlations of OTU richness across samples (hypothesis 2), we used linear regression of rarefied richness (as above) or, for plants, observed richness across the four domains/kingdoms and across dominant metazoan classes/orders. Similarly, correlations of community composition were tested using procrustes rotations and the function “protest” in the vegan package of R (Oksanen et al. ). Procrustes rotations find the maximum similarity between two configurations, and protest tests the significance of the resulting correlation.

To test hypothesis 3, we determined which abiotic drivers were best correlated with both species richness and species composition of all taxonomic groups. For richness, general linearized models (GLMs) were fit within each driver group (geomorphology, soil properties, and dominant plant species) using GLM, and then, a stepwise model selection was conducted using stepAIC to find the most parsimonious model. Dominant plant species were considered to be those with average cover of >2.5% across all plots. The penalty for degrees of freedom (k) was set to 4, insuring that most retained variables were significant. The significance of retained predictors was tested using the function drop1, with an F test for significance at P < 0.05. After finding the best individual predictors within each driver group, we then combined these predictors to test whether they were retained in a full model. This two‐step procedure was necessary, given the large number of potential predictors. For species composition, we used the function bioenv in vegan to find the best subset of environmental variables to maximize rank correlations with community dissimilarities. The best subset was found within geomorphology, soil properties, and dominant plants, followed by testing whether these predictors were retained in a full model. Results were visualized using the function envfit.

Results

DNA sequencing and taxonomic assignment

A total of 1,573,388 reads were obtained (mean of 52,446 per plot, range of 29,172–86,313 per plot). The number of reads that passed quality control and clustered to OTUs was 134,049 for the bacterial 16S primers, 294,313 for the fungal ITS primers, 248,586 for the nematode 18S primers, 37,664 for the metazoan 18S primers, and 124,897 for the insect COI primers. The somewhat lower sequence count for metazoan 18S primers reflected many sequences with short read lengths, rather than low overall yield. Although we did not include a positive control in the sequencing, the results suggest cross‐contamination was minimal. For example, two fungal families known to be host‐specific to pines: Suillaceae and Rhizopogonaceae, formed 2.25% and 7.3% of sequences in pine forest, respectively, but formed 0 and <0.002% of sequences under other land uses, respectively.

Richness and composition across land use types (hypothesis 1)

Richness for each taxonomic group varied with land use type (Fig. ). Significant effects of land use type on richness were observed for bacteria (ANOVA, F = 4.46, df = 4,25, P = 0.0073), Metazoa (ANOVA, F = 2.01, df = 4, 24; P = 0.00047), plants (ANOVA, F = 5.34, df = 4, 25; P = 0.0030), arachnids (ANOVA, F = 3.48, df = 4, 25; P = 0.022), and nematodes (ANOVA, F = 4.23, df = 4, 25; P = 0.0094). However, there was no consistent pattern between taxa in regard to the trends across different land uses. For example, bacterial richness was lowest in natural forests, where metazoan richness was highest (Fig. ).

Community composition multidimensional scaling (MDS) plots showed relatively good discrimination of land use types across taxa (Fig. ). All taxa were significantly affected by land use (permutation multivariate analysis of variance results: bacteria: F = 4.2_4,25, R² = 0.402, P = 0.001; fungi: F = 2.58_4,25, R² = 0.292, P = 0.001; metazoans: F = 1.67_4,25, R² = 0.211, P = 0.001; plants: F = 4.88_4,25, R² = 0.439, P = 0.001; Collembola: F = 2.09_4,25, R² = 0.251, P = 0.003; arachnids: F = 1.39_4,25, R² = 0.182, P = 0.044; nematodes: F = 2.97_4,25, R² = 0.322, P = 0.001; insects: F = 1.65_4,24, R² = 0.216, P = 0.001). Overlap between communities of different land use types tended to be between vineyard and improved grassland (e.g., fungi), and unimproved grassland and planted forests (e.g., plants and nematodes). For fungi, animals, plants, nematodes, and insects, the natural forest polygons were completely separated from the polygons of other land use types (Fig. ). Strong patterns within individual taxa supported these results (Appendices S2 and S3). In bacterial communities, for example, there were a relatively high abundance of Solibacterales and Chromatiales in natural forest and a low abundance of Gemmatimonadales and MC47, which were proportionally more abundant in all other land use types. Similarly, the fungal community of natural forest was dominated by Cortinariaceae, Hyaloscyphaceae, and Russulaceae which were largely absent from other land uses. Nematode communities in the grasslands were dominated by the plant pathogenic Hoplolaimidae and Tylenchulidae (the latter also being present in vineyards), while the soil bacteria‐feeding Cephalobidae were relatively common in each land use type. Insects stood out as showing strong separation of natural forest from all other land uses, driven by multiple taxa found almost only in this land use (e.g., Tephritidae, Sarcophagidae, and Staphylinidae).

Correlation of richness and composition between taxa (hypothesis 2)

There was poor correlation of richness between pairs of taxonomic groups. Fungal richness correlated with animal and plant richness (Fig. ), yet there were no significant correlations in species richness for class‐ and order‐level metazoan groups (Fig. ). In contrast, there was strong correlation of community composition between pairs of taxonomic groups. For domain‐ and kingdom‐level taxa, pairwise correlation coefficients ranged from 0.49 (Metazoa vs. bacteria) to 0.70 (Metazoa vs. plants) and all pairwise correlations were significant (P = 0.001; Fig. ). For class‐ and order‐level animal taxa, pairwise correlation values coefficients ranged from 0.16 (insects vs. Collembola) to 0.42 (nematodes vs. Collembola) and, again, all pairwise correlations were significant (Fig. ).

Drivers of taxonomic richness and composition (hypothesis 3)

There was little evidence of consistent drivers of taxonomic richness. Considered independently of land use, there were significant geomorphological, soil chemistry, and/or dominant plant species factors correlated with species richness in all taxonomic groups, but no single factor was included in the lowest Akaike information criterion (AIC) model across all taxonomic groups, and most were correlated with the richness of only one to three taxonomic groups (Table ). For example, soil parent material, mean annual rainfall, organic C, total N, and several plant taxa were significantly correlated with the richness of only one group each. Cation exchange capacity was the only factor correlated with three taxa, but was positively correlated with Collembola richness and negatively correlated with plant and bacterial richness.

Standardized linear coefficients of best environmental correlates of rarefied richness within each taxonomic group based on stepAIC within each category (geomorphology, soil chemistry, dominant plants) and then for the combination of the best predictors across the three categories

Characteristic	Kingdom/domain	Metazoan order/class
Bacteria	Fungi	Metazoans	Plants	Collembola	Nematodes	Arachnida	Insects
Geomorphology
Slope	–	–	0.4662***	–	–	0.486***	–	–
Mean annual rainfall	–	0.6382***	–	–	–	–	–	–
Soil (schist)	−1.398***	–	–	–	–	–	–	–
Soil (alluvial)	−0.8034***	–	–	–	–	–	–	–
Soil (Greywacke)	0.6743***	–	–	–	–	–	–	–
Soil chemistry
Organic carbon	–	–	–	–	−1.763*	–	–	–
Total nitrogen	–	–	–	–	1.292*	–	–	–
Olsen phosphorus	–	–	–	−1.073***	0.6389*	–	–	–
Total phosphorus	–	–	–	0.7696**	−1.289**	–	–	–
Calcium	–	–	−0.3744**	–	−1.542*	–	–	–
Magnesium	–	–	–	0.4088*	−0.5757	–	–	–
Potassium	–	–	–	0.827***	–	–	–	–
Cation exchange capacity	−0.6798***	–	–	−0.4032*	1.656	–	–	–
Base saturation	–	–	–	−0.6529*	1.821*	–	–	–
Dominant plant cover
Agrostis capillaris	–	–	–	–	0.2741	–	–	–
Kunzea ericoides	–	–	–	–	–	–	–	−0.3994*
Lolium perenne	–	–	–	–	–	–	−0.4754**	–
Muehlenbeckia complexa	–	0.4221**	0.4454***	–	–	–	–	–
Nothofagus cliffortioides	−0.5117***	–	–	–	0.4732*	–	–	–
Nothofagus fusca	–	–	–	–	0.4128	–	–	–
Pseudotsuga menziesii	–	−0.3028*	–	–	–	0.3948**	–	–
Trifolium repens	–	0.3731*	–	–	–	–	–	–
Trifolium subterraneum	–	–	–	−0.5127***	–	–	–	–
Vitis vinifera	−0.3214**	–	–	–	–	–	–	–
Weinmannia racemosa	–	–	–	–	–	0.3298*	–	–

Notes

Dashes indicate not included. Significance for geomorphology is for the parent material factor, not each level.

*P < 0.05, **P < 0.01, and ***P < 0.001 in F tests of models with and without that variable (function drop1).

Unlike richness, there were significant drivers of community composition shared across many groups. Soil carbon‐to‐nitrogen ratio and plot slope, for example, were significantly correlated with the composition of all four domain/kingdom‐level groups, as well as Collembola and Arachnids (Fig. ). There were also drivers that appeared more specific to individual taxa. The abundance of Pinus radiata, for example, was retained in the best models for fungi, metazoans, and Collembola, but not other groups, while Trifolium subterraneum was retained in the best model only for bacteria. Curiously, nematodes responded only to plant drivers, while arachnids were the only group to not respond to any dominant plant.

Discussion

More than 50% of the global land surface is estimated to be human‐modified (Hooke et al. ), with potentially substantial effects on local biodiversity (Newbold et al. , Gossner et al. ). While land use intensification has been viewed as having generally negative effects on biodiversity, our results suggest that different taxa respond divergently to land use. These inconsistencies were reflected in our driver analysis, with few environmental correlates of species richness in common across taxa. Species composition, in contrast, was strongly correlated across taxa, with a number of the same environmental drivers being significant for most. Drivers of composition were not the same as the drivers of diversity, suggesting diversity and composition are independently controlled.

Responses of richness and community composition to land use

As would be expected, we found that all taxonomic groups responded to land use; however, the direction of some responses was unexpected. It has often been demonstrated that natural ecosystems or less intensively managed ecosystems have higher species richness than highly modified ecosystems (Allen et al. , Newbold et al. ), but we found that this was only true for some groups. Others, for example, bacteria, had the lowest richness in natural forest and the highest richness in more intensively managed systems. Even for plants, the highest diversity was in managed unimproved grasslands rather than natural forest, and the intensively managed vineyards showed a similar diversity to natural forest.

Why were our results so distinct from prior suggestions of biodiversity loss with land use intensification (Allen et al. , Newbold et al. , Gossner et al. )? In part, this may be because our measurements are across multiple land uses, whereas some previous studies that have shown consistent negative effects of intensification on the diversity of multiple taxa examine single land uses and, in particular, grasslands (Allen et al. , Gossner et al. ). Indeed, if we compare grasslands within our data, we also see a reduction in the diversity of all taxa with intensification. Interestingly, this pattern is not evident between the natural and managed forested systems (Fig. ). This shows that patterns across multiple land uses may give a strikingly different pattern than looking only within grasslands.

Another reason why our results may differ is that previous studies may have focused on different taxonomic groups and potentially were biased toward those likely to show responses (e.g., plants; Kleijn et al. ), or may have omitted or down‐weighted the importance of diverse groups such as bacteria (e.g., Manning et al. ). Even in our study, bacterial diversity was addressed with a single primer pair, and bacterial OTUs defined based on this are likely higher in taxonomic rank than OTUs resolved for other taxa. Compared to other studies, our results may also be more strongly influenced by belowground biota, whose diversity response to land use intensification appears to be less consistently negative than that of aboveground biota (Gossner et al. ). Notwithstanding this, our results suggest that the impacts of land use intensification on biodiversity are taxon‐, land use‐, and scale‐dependent.

In contrast to the variable responses of richness to land use among taxa, the significant correlations in Procrustes rotations (Figs. , ) reflect the fact that all land uses had similar rank order in ordinations of the composition of different taxa (Fig. ). In all eight taxonomic groups, this followed the order: natural forest, planted forest, unimproved grasslands and then (in various orders) improved grasslands, and vineyards. The consistent responses of taxa with land use seem to support an almost Clementsian view of communities responding as coherent units. This may reflect the large scale of change between the land uses measured. The observed rank order also follows a rough gradient in disturbance frequency (a successional gradient), with high‐producing grasslands and vineyards both receiving high anthropogenic inputs and frequent soil disturbance, and progressively lower inputs and disturbance frequency in unimproved grasslands, planted forest, and natural forest. These results are consistent with other studies on successional gradients such as pasture to forest (Mueller et al. ) and cultivated cropland to forest (Lauber et al. ) suggesting that different land uses may be viewed as different successional states, defined by their disturbance regimes and functional composition. Our understanding of biodiversity and ecosystem responses to land use may therefore be enhanced by considering land use through the lens of ecological succession (Odum ).

Some of the observed changes of composition with land use were expected, for example, members of the Acidobacteria phylum (e.g., Solibacterales, Acidobacteriales) and many ectomycorrhizal fungi (e.g., Cortinariaceae and Russulaceae in native forest and Rhizopogonaceae in plantation forest) along with some decomposer fungi (Hyaloscyphaceae) being most abundant in forest ecosystems (Lauber et al. ); the increase in particular plant‐feeding nematodes in more managed ecosystems and the dominance of Tephritidae, Sarcophagidae, and Ditomyiidae flies in natural forest where there may be a greater amount of decomposing organic matter to feed on (Sousa et al. ). Other changes were more unexpected. For example, Burkholderiales, the dominant order in the beta‐Proteobacteria across plots, are thought to be copiotrophic (having a relatively fast growth rate, high‐nutrient requirement, and preference for using labile, rather than recalcitrant, carbon; Fierer et al. ), yet had a higher relative abundance in natural forest than in the improved grassland and vineyard sites.

Drivers of change in richness and composition

Our results indicated that the species richness of the taxonomic groups studied here was influenced by geomorphology, soil properties, and by dominant plant species identity. However, in accordance with the divergent responses of each taxonomic group to land use, the drivers of each group rarely overlapped (Table ). This is consistent with other studies that have also found drivers of richness to be taxon‐specific (Philpott et al. , Mueller et al. ). Environmental heterogeneity, which has been termed a universal driver of species richness (Stein et al. ), did not appear to be a strong driver of the richness of our taxonomic groups, as plant species richness, which is one of the main causes of environmental heterogeneity, was not correlated with the richness of any other group. This, combined with the lack of correlation among the richness of most other groups suggests that bottom‐up control was weak and that species interactions were general rather than specialized for our taxa (Hutchinson , Scherber et al. ). Overall, our results suggest that accurately predicting the responses of diversity of different taxonomic groups to land use requires a detailed knowledge of the complex interactions between the taxon of interest, its environment, and the other taxonomic groups interacting with it. In light of this, striving for a universal rule to predict the effects of land use change and intensification on biodiversity may ultimately be an unachievable aspiration. Further, the wide variation in responses suggests that proposed cross‐taxon measures of diversity (e.g., multidiversity; Allen et al. ) may obscure important changes in diversity within individual taxonomic groups.

In contrast to the results for richness, the composition of all taxonomic groups was correlated (Figs. , ), and soil C:N ratio was a significant driver for many taxa (Fig. ), along with the abundance of dominant early‐successional woody shrubs (Kunzea ericoides) and a suite of base cation‐related drivers (pH, Ca, K, cation exchange capacity, base saturation). These drivers are consistent with the idea that disturbance (e.g., soil disturbance and anthropogenic inputs) is the underlying mechanism differentiating slow‐cycling, high C:N ecosystems, which favor immobilization, bacteria, and fungi adapted to low‐nutrient environments, and limits plant available nutrients, from faster cycling low C:N systems with greater bacterial dominance and diversity (Odum , Vitousek et al. , Leff et al. ). The change in nutrient cycling and bacterial and fungal dominance, in turn, has potential implications for animal communities which feed on bacteria and fungi (e.g., nematodes and Collembola) and subsequent trophic levels. The abundance of K. ericoides, which is notable as an indicator of secondary succession as well as being ectomycorrhizal and having prolific floral resources and phytochemicals (Wardle ), was also important to fungi and several animal groups.

Of the geomorphological drivers, slope was a particularly key factor. Slope is an important determinant of land use (Mohammad ; Fig. ), perhaps reflecting the fact that while soil chemistry or moisture can be modified, slope is generally intractable (short of terracing). As a driver, slope was somewhat orthogonal to C:N ratio and the main land use gradient (Fig. ), suggesting that it may reflect within land use responses. This may, again, reflect disturbance and soil nutrient status, as slope is a major determinant of disturbance frequency and intensity (e.g., steep slope failures during earthquakes or storms), modifies organic matter accumulation, and influences soil fertility (Swanson et al. ). Other drivers of composition were more taxon‐specific. Bacterial composition was driven by pH, in accordance with previous studies (e.g., Lauber et al. , Zhalnina et al. ). Elsewhere, ecological interactions appeared to drive composition of some taxonomic groups. For example, bacterial composition was correlated with the nodulating Trifolium subterraneum, reflecting a plant–symbiont relationship, while ectomycorrhizal Pinaceae (Pseudotsuga menziesii and Pinus radiata) were important drivers for fungal communities (Dickie et al. ).

High‐throughput sequencing of multiple taxonomic groups

As one of the first studies to use high‐throughput sequencing for this kind of assessment, the results of our study demonstrate the power of the approach for understanding patterns of richness and community composition in multiple taxonomic groups across different land use types. The potential benefits of using environmental DNA to study such ecological patterns are numerous (Taberlet et al. ). While the fieldwork required is cost‐comparable to that used in traditional surveys of biodiversity, the ability to resolve diverse taxa without the specialist expertise and time required for morphological identification of each group is a clear benefit (Yu et al. , Drummond et al. ). Moreover, the greater sequencing depth offered by high‐throughput sequencing over other molecular methods allows greater resolution of ecological communities and patterns.

Implications for biodiversity management

By characterizing diversity within a broad range of taxonomic groups across multiple land use types, we have demonstrated inconsistencies in the responses of different taxa to land use. There was not always lower diversity in the more intensively managed land use types. Therefore, to the extent that maintaining regional biodiversity is considered desirable, our results support the idea that conserving a broad range of environments or habitats (i.e., an ecosystem approach) within regions, rather than particular suite of species, is likely to preserve the greatest overall diversity including both cryptic species (bacteria, fungi, soil invertebrates) and larger organisms such as plants. This is in line with suggestions that different taxonomic groups require different management strategies, and conserving hotspots for species richness (i.e., biodiversity) in one taxon would not necessarily imply benefits for maintaining high diversity within other taxa (Lambeck , Andelman and Fagan , Roberge and Angelstam ). However, consideration also needs to be given to whether particular taxa comprising the diversity are desirable or undesirable for either conservation or productive land uses, that is, whether the biodiversity measured includes exotic species, pests, and potential pathogens.

Acknowledgments

This research was funded by the New Zealand Ministry of Business, Innovation and Employment (PROP‐29430‐ESI‐LCR). K. Drew, E. Hayman, S. Kruis, H. Maule, and G. Walls provided field assistance, and D. Park assisted with DNA analyses. We thank the landowners of the study sites. P. Bellingham provided helpful comments on the manuscript. JRW, RJH, and IAD designed the study; RJH, CM, KIB, and IAD undertook fieldwork and sample collection; KHO, KIB, CD, and NB subsampled soil samples and performed function measurements; JRW, CD, and NB performed molecular analyses; RJH and IAD performed bioinformatic and statistical analyses; JRW, RJH, KHO, and IAD wrote the manuscript.

Data Accessibility

Pipeline and OTU tables summarizing DNA results are provided online at the Landcare Research Datastore (https://datastore.landcareresearch.co.nz/dataset/ngba-phase1).