1. Ecology

Rarity is a more reliable indicator of land-use impacts on soil invertebrate communities than other diversity metrics

  1. Andrew Dopheide  Is a corresponding author
  2. Andreas Makiola
  3. Kate H Orwin
  4. Robert J Holdaway
  5. Jamie R Wood
  6. Ian A Dickie
  1. Manaaki Whenua - Landcare Research, New Zealand
  2. Bio-Protection Research Centre, Lincoln University, New Zealand
  3. Manaaki Whenua – Landcare Research, New Zealand
  4. Bio-Protection Research Centre, School of Biological Sciences, University of Canterbury, New Zealand
https://doi.org/10.7554/eLife.52787
Share this article
Cite this article
  1. Andrew Dopheide
  2. Andreas Makiola
  3. Kate H Orwin
  4. Robert J Holdaway
  5. Jamie R Wood
  6. Ian A Dickie
(2020)
Rarity is a more reliable indicator of land-use impacts on soil invertebrate communities than other diversity metrics
eLife 9:e52787.
https://doi.org/10.7554/eLife.52787
  2. Download BibTeX
  3. Download .RIS

Abstract

The effects of land use on soil invertebrates – an important ecosystem component – are poorly understood. We investigated land-use impacts on a comprehensive range of soil invertebrates across New Zealand, measured using DNA metabarcoding and six biodiversity metrics. Rarity and phylogenetic rarity – direct measures of the number of species or the portion of a phylogeny unique to a site – showed stronger, more consistent responses across taxa to land use than widely used metrics of species richness, effective species numbers, and phylogenetic diversity. Overall, phylogenetic rarity explained the highest proportion of land use-related variance. Rarity declined from natural forest to planted forest, grassland, and perennial cropland for most soil invertebrate taxa, demonstrating pervasive impacts of agricultural land use on soil invertebrate communities. Commonly used diversity metrics may underestimate the impacts of land use on soil invertebrates, whereas rarity provides clearer and more consistent evidence of these impacts.

eLife digest

Living within the Earth’s soil are millions of insects, worms and other invertebrates, which help keep the ground healthy and fertile. There is a growing concern that changing land-use habits, such as agriculture and urban development, are causing these populations of invertebrates to decline. However, to what extent different types of land use negatively impact soil invertebrates is not clear.

Healthy habitats often have a greater variety of species. This biodiversity can be measured in a number of ways, ranging from counting the number of species, to more complex approaches that calculate a species’ role in an ecosystem or how close it is to extinction. Finding a way to sensitively measure the biodiversity of soil invertebrates could further researcher’s understanding of how different types of land use are affecting these communities.

A new method known as DNA metabarcoding has made it easier to distinguish between different species and calculate the biodiversity of entire populations. Now, Dopheide et al. have used this technique to study invertebrate communities from 75 sites across New Zealand which have been impacted by different land-use habits. This revealed that the most reliable and consistent way to uncover how land use affects soil invertebrates was to measure the rarity of species (i.e. the number of unique species present at each site).

Dopheide et al. found that agriculture negatively affected soil invertebrates and that most types of invertebrates responded in a similar way. Horticulture – such as orchards and vineyards – had the most severe impact, with the lowest variety of species compared to grassland or forest.

Other measurements of biodiversity, such as the number of different species, may underestimate the negative impact agriculture is having on invertebrate communities. The findings of Dopheide et al. highlight why developing strategies to preserve and restore these communities is so important. However, more work is needed to understand what specifically is causing biodiversity to decline and how this effect can be reversed.

Introduction

Land-use changes through deforestation, agricultural development, and urbanisation have caused worldwide impacts on the biodiversity of terrestrial communities and ecosystems (Dirzo et al., 2014; Newbold et al., 2015). Invertebrates are the most diverse and abundant component of animal biodiversity worldwide and are major contributors of terrestrial ecosystem services such as pollination, soil formation, and nutrient cycling (Lavelle et al., 2006; Wagg et al., 2014; Yang and Gratton, 2014). Long-term declines in the richness and biomass of insects and other terrestrial invertebrates are predicted to have major impacts on food webs and ecosystem functions (Eisenhauer et al., 2019; Hallmann et al., 2017; Potts et al., 2010). Despite this, most invertebrate species remain undescribed, and there is an incomplete understanding of land-use effects on invertebrate biodiversity, particularly for those that reside in soils (Cameron et al., 2018; Eisenhauer et al., 2019).

Biodiversity loss is typically measured as reductions in species richness (i.e., total number of species; e.g. Forister et al., 2010; George et al., 2019; Newbold et al., 2015). Despite widespread concern about biodiversity loss, evidence for impacts of anthropogenic land use on terrestrial invertebrate species richness is mixed, with studies often detecting richness declines for some taxa or groups but not others (Allan et al., 2014; Attwood et al., 2008; Blaum et al., 2009). Among the few studies that have examined land-use impacts on below-ground invertebrate communities, one detected negative impacts of long-term disturbance on soil invertebrate richness (Callaham et al., 2006), another detected increasing alpha diversity and homogenisation of soil invertebrates with increasing grassland intensification (Gossner et al., 2016); while others detected inconsistent richness patterns among different soil invertebrate taxa across land uses (George et al., 2019; Tsiafouli et al., 2015; Wood et al., 2017). These inconsistent patterns make it difficult to draw general conclusions about the impacts of land use on soil invertebrate biodiversity (Allan et al., 2014), and make the use of individual taxa as bioindicators problematic (Gerlach et al., 2013).

Inconsistent patterns in biodiversity measurement may reflect limitations of the diversity index used. In particular, species richness provides no indication of the distribution, taxonomy or function of species or communities (Fleishman et al., 2006; Hillebrand et al., 2018), potentially overlooking the nature and extent of land-use impacts on soil invertebrate communities. In contrast, rarity (sometimes termed ‘endemism richness’; Kier et al., 2009) measures the extent to which species are widely distributed generalists or limited to particular sites or land-use types. Rarity may thus indicate homogenising effects of land use on communities (McKinney and Lockwood, 1999; Smart et al., 2006), and the conservation value of sites (Kier and Barthlott, 2001). Furthermore, rare species can contribute disproportionately to ecosystem functioning (Dee et al., 2019; Leitão et al., 2016; Lyons et al., 2005; Mouillot et al., 2013). Rarity may therefore more accurately reflect the impacts of land use on soil invertebrate communities than species richness.

Rarity and other diversity metrics can also be placed in a phylogenetic context. Phylogenetic diversity reflects the evolutionary history and taxonomic range of communities and associated traits and functions (Faith, 1992), thus providing robust information for conservation assessment purposes (Faith, 1992; Forest et al., 2007; González-Orozco et al., 2015; Mishler et al., 2014). Phylogenetic diversity can also act as a proxy for functional diversity, albeit imperfectly (Mazel et al., 2018; Srivastava et al., 2012; Winter et al., 2013). Phylogenetic rarity, calculated as the portion of a phylogeny that is unique to a region or habitat (Mishler et al., 2014; Rosauer et al., 2009), combines elements of both rarity and phylogenetic diversity; high phylogenetic rarity implies that a community contains a taxonomically distinct assemblage of species and associated ecosystem functions. Mean pairwise distance, meanwhile, measures the phylogenetic relatedness of species within a community, which may reflect land-use driven filtering or competitive exclusion processes (Webb et al., 2002). The additional information represented by rarity and phylogenetic biodiversity metrics suggests that land-use related patterns based on these values may be clearer and more consistent among soil invertebrate taxa than those based on species richness and other non-phylogenetic diversity measures. Furthermore, rarity and phylogenetic rarity may be more sensitive indicators of land-use impacts on soil invertebrate communities than richness or phylogenetic diversity, because the former metrics reflect the distribution of species and lineages whereas the latter do not. These possibilities remain untested.

Here we present a comprehensive analysis of soil invertebrate biodiversity across different land-use types at a national spatial scale. We use modern DNA metabarcoding methods to measure invertebrate responses, as this enables the rapid and detailed identification of large numbers of invertebrate specimens from multiple taxonomic groups simultaneously (Drummond et al., 2015; George et al., 2019; Wood et al., 2017; Yu et al., 2012) and allows more efficient calculation of biodiversity metrics than previously possible. We analysed the invertebrate faunas in soil samples collected from 75 sites distributed across five different major land-use categories (natural forest, planted forest, low-producing and high-producing grassland, and perennial cropland) throughout New Zealand. Based on these data, we calculated six different biodiversity metrics: species richness, effective species numbers, rarity, phylogenetic diversity, phylogenetic rarity, and mean pairwise distance; as well as standardised effect size (SES) values for the latter phylogenetic metrics. We used these metrics to assess the impacts of land use on a comprehensive range of soil invertebrate taxa. We tested the following hypotheses: 1) all soil invertebrate taxa show the same biodiversity trends across the five land-use types; 2) patterns of soil invertebrate rarity, phylogenetic diversity, and phylogenetic rarity across the five land-use types are more consistent among taxa than species richness or non-phylogenetic diversity; 3) rarity and phylogenetic rarity of soil invertebrates are more sensitive to land use than richness, diversity, or phylogenetic diversity.

Results

Overall community composition

We detected a total of 11,284 operational taxonomic units (OTUs), of which 4549 (40.3%) were identified as terrestrial invertebrates. The remainder were identified as protists (37.6%), fungi (14.9%), non-terrestrial metazoans (5%), bacteria (1.7%), and plants (0.5%). The terrestrial invertebrate OTUs mostly belonged to the phylum Arthropoda (2,626 OTUs, among which insects were most common), followed by Rotifera (772 OTUs), Nematoda (656 OTUs), Mollusca (219 OTUs), Annelida (204 OTUs), Platyhelminthes (44 OTUs), Tardigrada (22 OTUs), Gastrotricha (four OTUs), and Onychophora (two OTUs) (Appendix 1—figures 1 and 2).

Non-metric MDS ordinations showed clear differences between overall invertebrate community composition in samples from different land-use categories (Figure 1). Natural forest samples formed a distinct cluster with no overlap with any other land-use categories. Samples from the other four land-use categories overlapped, with planted forest communities most similar to those from low-producing grassland followed by high-producing grassland communities, and least similar to those from perennial cropland. Similar trends were observed when only Arthropoda, Mollusca, Nematoda, or Rotifera OTUs were included, whereas Annelida OTUs showed less distinction between land-use categories. PERMANOVA tests for composition differences among different land-use categories detected a significant difference based on the overall invertebrate community (F4,61 = 1.804, p≤0.001), and based on each of the main phyla detected (Annelida, Arthropoda, Mollusca, Nematoda and Rotifera; F4,44-61 = 1.447–2.288, p≤0.001; Figure 1—source data 1A).

Figure 1 with 3 supplements see all
Download asset Open asset
Soil invertebrate community composition differs between land-use categories.

Non-metric MDS ordinations showing differences in the composition of soil invertebrate communities detected by DNA metabarcoding in five land-use categories, for overall communities, and for individual phyla with ≥ 100 OTUs. Ordinations are based on binary Jaccard distances.

Figure 1—source data 1

Results of PERMANOVA tests for differing soil invertebrate community composition, and ANOVA tests for differing multivariate homogeneity of sample dispersions, beta diversity, and phylogenetic beta diversity, between five land-use categories.

https://cdn.elifesciences.org/articles/52787/elife-52787-fig1-data1-v1.docx

To test for homogenisation effects of land use on soil invertebrate communities we compared multivariate heterogeneity/homogeneity of sample dispersions, mean pairwise beta diversity, and mean pairwise phylogenetic beta diversity, between land-use categories. For overall invertebrate communities, each of these measures differed significantly among land uses (F4, 61-442 = 3.59–14.99, p≤0.011), being highest in natural forest sites and lowest in grassland and/or cropland sites (Figure 1—source data 1B; Figure 1—figure supplements 1–3). Similar trends were observed for Arthropoda and Nematoda communities based on all three measures, and for Annelida and Mollusca communities based on phylogenetic beta diversity and multivariate heterogeneity of sample dispersions, whereas Rotifera communities showed different patterns.

A heatmap based on the 1000 most relatively abundant terrestrial invertebrate OTUs detected suggested that low-producing grassland, high-producing grassland, and perennial cropland samples each had relatively consistent assemblages of abundant OTUs, both within and between each land-use category, whereas planted forest samples, and especially natural forest samples, each had more distinctive assemblages of abundant OTUs (Figure 2 and Figure 2—figure supplement 1). In particular, most of the natural forest samples had a subset of abundant OTUs that were not detected in any other sample.

Figure 2 with 1 supplement see all
Download asset Open asset
Distribution of the 1000 most abundant soil invertebrate OTUs across samples and land-use categories.

The proportional abundance and distribution among samples and five land-use categories of the 1000 most proportionally abundant soil invertebrate OTUs detected by DNA metabarcoding, showing that natural forest sites have more heterogeneous assemblages of soil invertebrate OTUs than agricultural sites. Samples are ordered on the x-axis by land-use category and increasing latitude.

Overall invertebrate biodiversity differences among land-use categories

All biodiversity metrics (except for mean pairwise distance) showed a general trend of declining overall invertebrate biodiversity (i.e. the biodiversity of the entire invertebrate community) from forested and/or low-producing grassland sites to high-producing grassland and/or perennial cropland sites (Figures 3 and 4). Rarity and phylogenetic rarity metrics showed the largest and most consistent land-use-related biodiversity declines, with the highest mean values in natural forest sites followed by planted forest sites and low-producing grassland sites, and high-producing grassland sites, and lowest values in perennial cropland sites. Removing species found in only a single site did not substantially change these trends (Appendix 1—figures 3–5). Significant differences between mean biodiversity of overall invertebrate communities in different land-use categories were detected according to richness, rarity, phylogenetic diversity, phylogenetic rarity, and phylogenetic diversity and rarity SES metrics (F4,64 = 3.56 to 17.986, p = 0.012 to <0.001), but not effective species numbers, mean pairwise distance, or mean pairwise distance SES metrics (Figure 3—source data 1A; Figure 4—source data 1A). ANOVA tests of derived land-use rank trends provided similar results, with significant trends identified for all metrics except for mean pairwise distance and mean pairwise distance SES (F1,67 = 4.66–31.94, p = 0.034 to <0.001; Appendix 1—table 1).

Figure 3 with 2 supplements see all
Download asset Open asset
Biodiversity estimates for overall soil invertebrate communities detected in different land-use categories.

The biodiversity of soil invertebrate communities detected by DNA metabarcoding declines from forested to agricultural sites according to most metrics, with the clearest declines shown by rarity metrics. Diamonds and whiskers represent mean values ± standard errors, with individual data points represented by circles. ANOVA test statistics and trend splines are shown for cases with statistically significant biodiversity differences among land-use categories, with letters indicating differences between land-use categories detected by post-hoc Tukey HSD tests.

Figure 3—source data 1

Results of ANOVA tests for differing soil invertebrate biodiversity between different land-use categories, according to six biodiversity metrics.

https://cdn.elifesciences.org/articles/52787/elife-52787-fig3-data1-v1.docx

The mean rarity of overall invertebrate communities was significantly lower in all four other land uses compared with natural forest (t23-27 = −31.6 to −62.4, P.adj = 0.03 to <0.001). Similarly, the mean phylogenetic rarity of overall invertebrate communities was significantly lower in all four other land-use categories compared with natural forest (t23-27 = −3.34 to −6.90, P.adj = 0.043 to <0.001), and in perennial cropland compared with planted forest (t24 = −3.55, P.adj = 0.046). In contrast, the mean richness and phylogenetic diversity of overall invertebrate communities were similar in natural forest, planted forest, and low-producing grassland samples, and significantly lower in perennial cropland compared with natural forest (t23 = −78.3, P.adj = 0.023, and t23 = −14.6, P.adj = 0.008, respectively) and compared with low-producing grassland (t23 = −84.2, P.adj = 0.012, and t23 = −13.3, P.adj = 0.019, respectively; Figure 3). Mean phylogenetic diversity SES was significantly lower in low-producing grassland compared with natural forest (t23-27 = −2.20, P.adj = 0.048), but did not otherwise differ between land-use categories, while phylogenetic rarity SES differences between land-use categories matched those based on non-SES phylogenetic rarity (t23-27 = −3.68 to −8.61, P.adj = 0.031 to <0.001; Figure 4).

Figure 4 with 2 supplements see all
Download asset Open asset
Phylogenetic biodiversity SES estimates for overall soil invertebrate communities detected in different land-use categories.

Phylogenetic biodiversity SES estimates for soil invertebrate communities detected by DNA metabarcoding tend to decline from natural forest to agricultural sites, with the clearest decline shown by phylogenetic rarity SES. Diamonds and whiskers represent mean values ± standard errors, with individual data points represented by circles. ANOVA test statistics and trend splines are shown for cases with statistically significant biodiversity differences among land-use categories, with letters indicating differences between land-use categories detected by post-hoc Tukey HSD tests.

Figure 4—source data 1

Results of ANOVA tests for differing soil invertebrate biodiversity between different land-use categories, according to three phylogenetic biodiversity SES metrics.

https://cdn.elifesciences.org/articles/52787/elife-52787-fig4-data1-v1.docx

A mixed-model ANOVA test for effects of derived land-use rank, land-use category, and taxonomic group effects showed that derived land-use rank and taxonomic group (and interactions) were the most consistently significant predictors of the diversity metrics (F1-16 = 7.74 to 32.14, p = 0.007 to <0.001; Appendix 1—table 2). The further addition of land-use category to models already containing derived land-use rank did not explain additional variation for effective species, rarity, phylogenetic rarity and mean pairwise distance, but did for richness and phylogenetic diversity (in the form of significant interactions between land-use category and taxonomic group; F48 = 1.41 and 1.82, p = 0.037 and <0.001).

Most environmental variables showed clear land use-related trends of increasing or decreasing values in the order of natural forest, planted forest, low-producing grassland, high-producing grassland, and perennial cropland (Appendix 1—figure 6). An ANOVA test of spatial attributes (latitude and altitude) plus land-use category showed latitude had no effect on overall soil invertebrate biodiversity according to any metric, whereas altitude had significant effects on biodiversity of all metrics except for mean pairwise distance (F1 = 9.41 to 22.33, p = 0.003 to <0.001). In addition to altitude, land-use category had a significant effect only on rarity and phylogenetic rarity metrics (F1 = 4.40 and 4.60, p = 0.003 and 002; Appendix 1—table 3). The first three components of a PCA incorporating latitude, altitude, and soil chemistry variables explained 70.25% of variance. According to an ANOVA test of these three PCA components plus land-use category, the first component had significant effects on the rarity, phylogenetic diversity and phylogenetic rarity of the overall soil invertebrate biodiversity (F1 = 4.79 to 15.25, p = 0.032 to <0.001), and the second component on the former three metrics plus richness (F1 = 7.00 to 10.24, p = 0.010 to 0.002). The third component did not have a significant effect on any of the metrics. The addition of land-use category to these models explained further variation for richness, rarity, and phylogenetic rarity metrics only (F4 = 2.71 to 4.72, p = 0.038 to 0.006; Appendix 1—table 4), indicating that there was some confounding between the environmental PCAs and land-use category.

Biodiversity differences among invertebrate taxa

Biodiversity metrics for the main insect orders (Coleoptera, Diptera, Hymenoptera, Lepidoptera, Hemiptera, and all other insects), other arthropod taxa (Collembola, mites, non-mite Arachnida, Malacostraca, myriapods), and non-arthropod phyla (Annelida, Mollusca, Nematoda, Platyhelminthes, Rotifera, and Tardigrada) that were detected showed a general trend of declining biodiversity from forested to agricultural sites. Rarity, phylogenetic diversity, and phylogenetic rarity patterns were most consistent among different taxonomic groups (Appendix 1—figures 7–12), while land-use trends were most clear and consistent across taxonomic groups according to rarity and phylogenetic rarity (Figure 3—figure supplements 1 and 2). ANOVA tests detected significant differences among land-use categories for ten of the 17 taxonomic groups based on rarity (all insect groups, non-mites, Annelida, Nematoda, and Platyhelminthes; F4 = 2.60 to 13.26, p = 0.048 to <0.001); nine groups based both on phylogenetic rarity (all insect groups except Hemiptera, mites and non-mites, Annelida, and Platyhelminthes; F4 = 2.74 to 11.07, p = 0.036 to <0.001) and phylogenetic diversity (all insect groups, Annelida, Mollusca, and Nematoda; F4 = 3.14 to 6.41, p = 0.047 to <0.001); eight groups based on richness (all insect groups, Nematoda, and Platyhelminthes; F4 = 2.55 to 6.32, p = 0.048 to <0.001); five groups based on effective species numbers (Diptera, Hymenoptera, Lepidoptera, mites, and Annelida; F4 = 2.73 to 4.36, p = 0.037 to 0.004); and three groups based on mean pairwise distance differences (Hymenoptera, mites, and Rotifera; F4 = 3.53 to 6.24, p = 0.012 to <0.001; Figure 3—source data 1B). Tests of derived land-use rank trends for each metric and taxonomic group provided concordant results, with the same groups (with few exceptions) showing significant trends for each metric (Appendix 1—table 5).

Post-hoc Tukey HSD tests showed that biodiversity was most commonly significantly higher in natural forest compared with perennial cropland (Figure 3—figure supplements 1 and 2). This was observed for nine taxonomic groups based on rarity (t14-23 = 1.92 to 7.31, P.adj = 0.040 to <0.001), eight groups based on phylogenetic rarity (t20-28 = 0.054 to 1.19, P.adj = 0.024 to <0.001), five groups based on phylogenetic diversity (t20-23 = 1.16 to 2.63, P.adj = 0.032 to <0.001), four groups based on richness (t14-23 = 3.69 to 9.47, P.adj = 0.026 to <0.001), three groups based on mean pairwise distance (t22-23 = 0.03 to 0.35, P.adj = 0.014 to 0.003), and just one group based on effective species numbers (t25 = 3.00, P.adj = 0.012). Biodiversity was also significantly higher in natural forest compared with high-producing grassland (for two to six groups according to each of five metrics; t21-27 = 0.02 to 6.86, P.adj = 0.029 to <0.001), low-producing grassland (one to five groups, four metrics; t20-26 = 0.04 to 4.71, P.adj = 0.040 to <0.001), and planted forest (one to three groups, three metrics; t24-27 = 0.64 to 4.61, P.adj = 0.041 to 0.007); in planted forest, low-producing grassland, or high-producing grassland compared with perennial cropland (one to two groups, two to five metrics; t12-24 = 0.38 to 16.92, P.adj = 0.045 to 0.001); and in planted forest or low-producing grassland compared with high-producing grassland (one or two groups, two metrics; t23-30 = 2.14 to 3.33, P.adj = 0.036 to 0.023). All of the pairwise differences together implied a land-use category rank order of natural forest > planted forest > low producing grassland > high producing grassland > perennial cropland.

Non-parametric bootstrapping of ANOVA sum of squares values for the (non-SES) biodiversity metrics and taxonomic groups for which significant land-use differences were detected showed that phylogenetic rarity followed by (non-phylogenetic) rarity explained the largest proportions of land-use category variance across the 17 taxonomic groups, while mean pairwise distance and richness explained the least variance (Figure 5). A Kruskal-Wallis test detected significant differences among the biodiversity metrics (Chi square = 4782.6, df = 5, p<0.001), with post-hoc tests indicating that the distributions of all metrics differed significantly from each other (p<0.05).

Figure 5
Download asset Open asset
Proportions of sample variance explained by land use according to different biodiversity metrics.

The proportions of sample variation (sum of squares) explained by land use were estimated for different biodiversity metrics by non-parametric bootstrapping, based on the combinations of biodiversity metric and soil invertebrate taxonomic group for which significant land-use differences were detected by ANOVA tests. Observed mean values and 95% confidence interval limits are indicated by orange and blue vertical bars, respectively.

Phylogenetic biodiversity metric SES differences among taxa

Patterns of phylogenetic rarity SES values among land-use categories were more consistent across taxonomic groups, and with their corresponding non-SES metric patterns, than patterns of phylogenetic diversity SES and mean pairwise distance SES values (Figure 4—figure supplements 1 and 2). ANOVA tests detected significant differences among land-use categories for 11 of the 17 taxonomic groups based on phylogenetic rarity SES (Collembola, Coleoptera, Diptera, Lepidoptera, other insects, mites and non-mites, Annelida, Mollusca, Nematoda, and Rotifera; F4 = 3.10 to 8.91, p = 0.022 to <0.001), six groups based on phylogenetic diversity SES (Hymenoptera, Lepidoptera, mites, Malacostraca, Nematoda, and Rotifera; F4 = 2.76 to 7.39, p = 0.035 to <0.001); and four groups based on mean pairwise distance SES (Lepidoptera, mites, Malacostraca, and Rotifera; F4 = 4.40 to 11.28, p = 0.016 to <0.001; Figure 4—source data 1B). All of the 11 taxonomic groups with significant phylogenetic rarity SES differences showed a consistent pattern of declining rarity from natural forest to planted forest to agricultural land-use categories. Post-hoc Tukey HSD tests detected significantly higher phylogenetic rarity SES values in natural forest (for 11 groups) and in planted forest (for four groups) compared with at least two of the agricultural land-use categories in each case (t22-28 = −1.73 to −3.77, P.adj = 0.047 to <0.001). In contrast, only two groups (mites and Rotifera) showed this pattern based on either phylogenetic diversity SES (t22-28 = −1.52 to −3.15, P.adj = 0.031 to <0.001) or mean pairwise distance SES values (t22-28 = −1.83 to −2.89, P.adj = 0.047 to <0.001). Otherwise, Lepidoptera phylogenetic diversity SES values were significantly lower in both planted forest and high-producing grassland compared with both natural forest and perennial cropland (t21-27 = −1.07 to −1.44, P.adj = 0.035 to 0.004), whereas Hymenoptera, Malacostraca and Nematoda phylogenetic diversity SES values were higher in one or more of the anthropogenic land use categories compared with natural forest (t3-28 = 1.46 to 2.87, P.adj = 0.020 to 0.005). Patterns of mean pairwise distance SES values across land use categories and taxonomic groups closely matched those observed for phylogenetic diversity SES values (except significant differences among land-use categories were not detected for Hymenoptera or Nematoda).

Discussion

This research provides clear evidence of adverse impacts of agricultural land use upon soil invertebrate communities. Effects of land use on biological communities are usually measured as shifts in species richness. However, rarity metrics were much more sensitive to land use and more consistent among taxa than richness or effective species numbers in our study, suggesting that the latter metrics may underestimate land-use impacts on biodiversity. Rarity is a function of the number of species with limited distributions or narrow habitat specificity. These rare species can have important roles in ecosystem processes (Dee et al., 2019; Leitão et al., 2016; Lyons et al., 2005), and are inherently more vulnerable to extinction. Overlooking species rarity, as richness does, therefore obscures the effects of different land uses on communities, with potential detrimental consequences for the function and resilience of ecosystems. Our results suggest that efficient DNA-based measurement of plot-level rarity improves our understanding of rare species occurrence and provides an effective basis for incorporating soil invertebrates into conservation planning.

Rare species include not only habitat specialists, but also transient and conditionally rare taxa. It is possible that OTUs that were rare in this study may be more common in locations not sampled. Nonetheless, our observation that patterns of rarity among land-use categories were the most consistent among different taxa suggests that rarity is an ecologically meaningful measure of ecosystem biodiversity. This is supported by prior studies suggesting that rare species are particularly sensitive to ecosystem change. For example, rare plant and fungal species appear to be particularly sensitive to changes in environment (Avis et al., 2008; Dickie et al., 2009; Dickie and Reich, 2005; McIntyre and Lavorel, 1994), and pollinating insect losses are concentrated among rare species (Powney et al., 2019). Furthermore, most of the terrestrial invertebrate species currently considered to be at risk or threat of extinction in New Zealand are naturally uncommon (Stringer and Hitchmough, 2012).

Phylogenetic diversity – and especially phylogenetic rarity – explained larger proportions of land-use variance across taxa than their non-phylogenetic counterparts, and phylogenetic rarity was overall the most sensitive metric to land-use differences. Phylogenetic metrics incorporate evolutionary and functional aspects of biodiversity (Faith, 1992; Faith, 2015; Mazel et al., 2018). New Zealand has a long history of geographic isolation and glaciation, reflected by the presence of many deeply divergent invertebrate lineages (Buckley et al., 2015; Trewick et al., 2011). The high levels of invertebrate phylogenetic rarity in natural forest sites likely reflects assemblages of long-present soil invertebrates that are highly adapted to these habitats, but ill-suited to the modified land-use types included in the study. These trends might differ in regions with greater connectivity, longer-term agriculture, and different geological history. Phylogenetic diversity SES and mean pairwise distance SES values showed different evidence of land-use effects compared with their non-SES counterparts, suggesting, for example, that Lepidoptera, mite and Rotifera communities are less dispersed, suggesting loss of lineages, in agricultural sites compared with forest habitats. In contrast, Malacostraca communities appear to be under-dispersed in natural forest sites, and to gain lineages due to anthropogenic land use. Phylogenetic rarity SES values further support the finding of consistently reduced rarity in agricultural sites, independent of species richness effects. Together, these observations indicate that phylogenetic information provides additional insights into soil invertebrate biodiversity patterns, as has been observed for other groups (González-Orozco et al., 2015; Mishler et al., 2014).

Land-use impacts

The low beta diversity, heterogeneity, and rarity values detected in agricultural sites, and the overlap of samples from these sites in MDS ordinations, together strongly imply that these habitats tend to have relatively similar assemblages of species across locations. Agricultural practices have effects at a wide range of scales, from local-scale use of chemical fertilisers and pesticides to landscape-scale habitat simplification (Tscharntke et al., 2005). Together these factors lead to homogenisation of communities and functions among sites, in which specialists in diverse natural communities are replaced by a smaller number of generalists that thrive in anthropogenic habitats (Börschig et al., 2013; Clavel et al., 2011; Gámez-Virués et al., 2015; McKinney and Lockwood, 1999; Smart et al., 2006).

In contrast to the agricultural sites, the high diversity and rarity observed in natural forest sites indicates that these habitats tend to have richer and more unique assemblages of species. Forested sites tend to have greater physical habitat complexity and heterogeneity, providing more varied resources and niches for diverse communities including various specialists (Jonsson et al., 2009; Stein et al., 2014). Furthermore, natural forest habitats tend to be more disconnected, and located in more rugged and less accessible areas than agricultural sites, with more physical barriers to limit the dispersal of invertebrate fauna. Consequently, the distinct assemblages detected in natural forest sites are likely to reflect natural historical biogeographic distribution and evolutionary processes (Buckley et al., 2015; Trewick et al., 2011).

Despite their varying sensitivity, most metrics of rarity and diversity (not mean pairwise distance, phylogenetic diversity SES, or mean pairwise distance SES) showed a consistent trend of lower biodiversity in agricultural land-use categories than in forested land-use categories. Further, while not all taxa showed significant evidence of declining biodiversity in relation to agricultural land use, no taxa responded positively. Many taxa not showing significant biodiversity declines had few species (e.g. myriapods, Malacostraca and tardigrades), suggesting there was insufficient data to infer land-use differences. Among the most species-rich groups that did not show significant declines (collembola, mites and rotifers), many of the diversity metrics were nonetheless lowest in grassland or perennial cropland sites, suggesting that while these groups may be more resilient to impacts of agricultural land use than others, the general trend was similar. These biodiversity declines are in contrast to previous research that suggested soil fauna are resilient to grassland intensification (Gossner et al., 2016), likely because our study encompasses a broader range of land-use types. While it is likely that spatial and environmental factors associated with particular land uses contribute to these patterns, the fact that land use explained additional variation of richness and rarity metrics after these factors were statistically accounted for strongly indicates an independent role of land management practices.

While rarity and phylogenetic rarity metrics showed the most consistent responses across land-use categories, the rank order of land-use categories implied by these (and other) metrics were not easily predicted prior to measurement. Planted forests, which were predominantly Pinus radiata plantations, are sometimes perceived as being biologically depauperate, while low-producing grasslands are frequently perceived as semi-natural in New Zealand (Hobbs et al., 2006). Despite this, we found rarity and diversity in planted forest sites to be similar to those in low-producing grassland sites and higher than those in high-producing grassland or perennial cropland sites, consistent with suggestions that plantations can play an important role in insect biodiversity conservation (Pawson et al., 2009; Pawson et al., 2010). Similarly, high-productivity grasslands are often perceived as a more severe land use than perennial cropland due to high homogeneity of vegetation cover, low habitat complexity, and high fertiliser use. Nonetheless, our data suggest perennial cropland supports the lowest levels of invertebrate diversity and rarity of any of the measured land-use categories. This may reflect high chemical input in and intensive management of fruit production systems (Manktelow et al., 2005).

Overall, our results suggest pervasive impacts of agricultural land use upon soil invertebrate communities, with likely adverse consequences for ecosystem services. This adds to widespread evidence of declines in invertebrate biomass and diversity in response to anthropogenic land-use change and habitat loss (Attwood et al., 2008; Hallmann et al., 2017; Hendrickx et al., 2007; Powney et al., 2019), and suggests that efforts to conserve and restore soil invertebrate communities may be needed.

Conservation implications

Invertebrates tend to be neglected by conservation initiatives, due to the challenges of determining their identities, functions, and distributions (Leandro et al., 2017). Indirect preservation of communities via flagship or umbrella species protection schemes tends to be ineffective (Andelman and Fagan, 2000; Oberprieler et al., 2019; Schuldt and Assmann, 2010), and similarly, biomonitoring based on individual species is problematic. By allowing the efficient assessment of invertebrate community composition and distribution across large spatial scales, DNA metabarcoding methods may enable more informative biomonitoring and improved targeting of conservation initiatives based on multiple invertebrate taxa, if not entire invertebrate communities. While rarity and phylogenetic rarity were the most informative metrics of community change in this case, it is likely that consideration of these alongside richness and phylogenetic measures of diversity would provide the most comprehensive information for purposes such as biomonitoring and conservation planning (Fleishman et al., 2006). Our results suggest that conserving a network of sites with high invertebrate diversity and rarity would preserve a diverse assemblage of species, communities, and functional traits, thus providing resilience of communities and ecosystem processes to environmental changes (Balvanera et al., 2006; Yachi and Loreau, 1999). While diversity and rarity was typically highest in our natural forest sites (of which many are protected), certain grassland and cropland sites with unusually high rarity values (outliers on Figure 3) might be logical targets for further investigation and potential incorporation into conservation initiatives.

In conclusion, our analysis of soil invertebrate biodiversity across land-use categories at a national scale shows that most soil invertebrate taxa have consistent rarity responses to land use, and that agricultural land use tends to cause the homogenisation and loss of soil invertebrate biodiversity. This research adds to evidence of widespread impacts of anthropogenic land use on invertebrate biodiversity, but also implies that these impacts may have been underestimated due to a widespread emphasis on species richness. DNA metabarcoding methods offer an efficient basis for measuring the diversity and rarity of invertebrate communities at large scales. Incorporating this information into conservation schemes would enable the protection of a broader range of biodiversity and enhance the preservation of terrestrial ecosystems.

Materials and methods

Key resources table
Reagent type
(species) or resource		DesignationSource or referenceIdentifiersAdditional
information
Sequence-based reagentmICOIintFDOI:10.1186/1742-9994-10-34GGWACWGGWTGAACWGTWTAYCCYCC
Sequence-based reagentHCO2198PMID:7881515TAAACTTCAGGGTGACCAAAAAATCA
Commercial assay or kitNucleoSpin Tissue kitMacherey-Nagel740741.4
Software, algorithmcutadapthttps://github.com/marcelm/cutadaptv 1.11
Software, algorithmUSEARCHhttps://www.drive5.com/usearch/v 9.0.2132_i86linux32
Software, algorithmVSEARCHhttps://github.com/torognes/vsearchv 2.4.0
Software, algorithmRhttps://www.r-project.org/v 3.52
Software, algorithmphylo.endemismhttps://davidnipperess.blogspot.com/2012/07/phyloendemism-r-function-for.html

Sample collection

Request a detailed protocol

Soil invertebrate communities were sampled from a total of 75 sites distributed across five different major land-use categories throughout New Zealand (Figure 6), during dry weather between November 2014 and March 2015. The five land-use categories (natural forest, planted forest, low-producing grassland, high-producing grassland, and perennial cropland) represent differing states of anthropogenic modification (Figure 6—source data 1). The site locations were selected from a nationwide 8 km grid used for regular monitoring of native species and pests. For each land-use category, 15 replicate sites were randomly selected from the nationwide monitoring grid, excluding any that were >1000 m altitude and ensuring they were distributed across the length of New Zealand (Makiola et al., 2019). At each site, a 20 m × 20 m plot was established according to a standardised protocol (Hurst and Allen, 2007). Twenty-four soil cores were collected within each plot on a regular grid (min 3.54 m distance between cores) to a depth of 15 cm using a sterile corer (5.08 cm diameter), following Wood et al. (2017). Surface litter was removed prior to coring. The 24 soil cores were pooled together, homogenised, and stored at 4°C until laboratory processing. Invertebrates were extracted from a one-litre subsample of homogenised soil material from each site using Berlese-Tullgren funnels and stored in ethanol until DNA extraction.

Figure 6
Download asset Open asset
Location and land-use category of 75 sample sites.

Site locations were randomly selected from a nationwide 8 km grid used for regular monitoring of native species and pests, excluding any that were >1000 m altitude and ensuring they were distributed throughout New Zealand. X- and y-axes represent longitude and latitude, respectively.

Figure 6—source data 1

Defining attributes of land-use categories.

https://cdn.elifesciences.org/articles/52787/elife-52787-fig6-data1-v1.docx

The altitude and latitude of plots were determined from topographic maps. Soil chemistry variables (pH, C, N, C:N ratio, Olsen P, Total P, Ca, Mg, K, Na, cation exchange capacity, base saturation) were determined for each plot according to Orwin et al. (2016) and Wood et al. (2017).

Molecular laboratory procedures

Request a detailed protocol

Bulk invertebrate concentrates were centrifuged for three minutes at 2,500 rpm (1258 rcf), after which ethanol was removed until <5 ml remained. The concentrates were then transferred into 5 ml tubes and homogenised with eight steel balls in a bead mill operated at 15 Hz for six intervals of 20 s each. A 1.5 ml aliquot of homogenised invertebrate concentrate from each sample was removed into a 1.5 ml microtube and centrifuged for one minute at 13,000 rpm (11,337 rcf), after which any ethanol was removed. The pelleted material was resuspended in purified water, re-centrifuged as before, then resuspended in 200 µl digestion buffer (10 mM Tris buffer, 10 mM NaCl, 5 mM CaCl2, 2.5 mM EDTA, 2% SDS, 0.04 M dithiothreitol, and 0.1 M proteinase K) with vortexing, and incubated overnight at 56 °C with shaking at 450 rpm (Campos and Gilbert, 2012). DNA was extracted from the digested samples using a Macherey-Nagel NucleoSpin Tissue kit (MACHEREY-NAGEL GmbH and Co. KG, Düren, Germany), omitting sample lysis steps but otherwise according to the manufacturer’s directions, with a JANUS workstation laboratory robot (PerkinElmer, Waltham, MA, USA). The DNA concentration was quantified in each extract using an Invitrogen Quant-iT PicoGreen dsDNA quantitation assay kit (Thermo Fisher Scientific, Waltham, MA USA), and standardised across samples to 3 ng/µl.

COI barcodes were amplified by PCR from each sample using metazoan-targeted primers mICOIintF (5'-GGWACWGGWTGAACWGTWTAYCCYCC-3') (Leray et al., 2013) and HCO2198 (5'-TAAACTTCAGGGTGACCAAAAAATCA-3') (Folmer et al., 1994), which were respectively modified at their 5' ends with the linker sequences 5'-TCGTCGGCAGCGTC-3' and 5'-GTCTCGTGGGCTCGG-3'. PCRs were carried out in 20 µl volumes, containing 200 nM of the forward and reverse COI primers, 0.2 mM of each dNTP, 1.5 mM MgCl2, 2 µg rabbit serum albumin, 0.5 U KAPA Plant 3G enzyme (Kapa Biosystems, Wilmington, MA, USA), and 2 µl (6 ng) DNA template. The PCR amplification protocol was 95 °C for 3 min; 35 cycles of 95 °C for 20 s, 52 °C for 15 s, and 72 °C for 30 s; and 1 min at 72 °C. Illumina sequencing adapters and sample-specific barcodes were added to the COI amplicons in a second round of PCR, carried out in 25 µl volumes containing the same reagents and concentrations as the first PCR, except for Illumina-tagged sequencing adaptors instead of COI primers, and 2 µl of the first PCR amplicon as template. The second-round PCR amplification protocol was 95 °C for 3 min; five cycles of 95 °C for 20 s, 54 °C for 15 s, and 72 °C for 30 s; and 1 min at 72 °C. The resulting libraries were purified and size-selected using a Pippin Prep system (Sage Science, Beverly, MA, USA), to remove primer dimers and high molecular weight DNA, quantified, pooled, and sequenced on an Illumina MiSeq system with a 2 × 250 sequencing kit at the Australian Genome Research Facility Ltd.

Bioinformatic processing

Request a detailed protocol

Demultiplexed forward and reverse DNA reads were merged and relabelled by sample using USEARCH (Edgar, 2013). Linker sequences and primers were trimmed from the merged sequences using cutadapt (Martin, 2011). The trimmed sequences were quality filtered to remove any with >1 maximum expected errors and dereplicated using VSEARCH (Rognes et al., 2016). Non-singleton sequences (i.e. those represented by at least two identical sequences) were clustered into OTUs at a sequence identity threshold of 97% and simultaneously filtered for chimeras using the UPARSE algorithm in USEARCH (Edgar, 2013). OTU abundance was inferred by mapping the trimmed sequences back to the OTU centroid sequences at a sequence identity threshold of 97%. The OTUs were assigned a taxonomic identity using the RDP Naïve Bayesian classifier (Wang et al., 2007) in combination with an RDP-formatted animal mitochondrial COI sequence database (Porter and Hajibabaei, 2018), which includes bacterial, fungal, and protist COI sequences to enable the detection of non-metazoan OTUs. We excluded any OTUs that were not identified as belonging to an expected terrestrial invertebrate phylum.

Biodiversity analyses and statistics

Request a detailed protocol

Data analyses were carried out using R version 3.5.1 (R Development Core Team, 2016) and RStudio (RStudio team, 2015). Extraction blanks, negative and positive controls were examined for contamination. Tag jumping (Schnell et al., 2015) was accounted for by using a regression of contaminant abundances versus the maximum of total abundances in all other samples, after which the coefficient estimate for the 90th quantile regression was used to subtract that many sequences from the abundances of all OTUs (Makiola et al., 2019).

Comparisons of multivariate community composition and homogeneity between land-use categories were carried out for the overall terrestrial invertebrate dataset and the main terrestrial invertebrate phyla detected using the R package vegan v2.4–3 (Oksanen et al., 2017). Non-metric MDS ordinations and PERMANOVA tests for community composition differences among land uses were based on the Jaccard distance metric and presence/absence data. Any samples with unusually low sequence abundance (defined as less than 5% of the mean sequence abundance per sample for a given phylum) were excluded from MDS ordinations. For the Mollusca-based MDS ordination, one further sample that resulted in an uninterpretable plot was excluded. To test for homogenisation effects of land use, multivariate homogeneity of sample dispersions was determined for each land-use category and compared between categories using the function betadisper in the R package vegan. Similarly, mean pairwise beta diversity and phylogenetic beta diversity (UniFrac distances; Lozupone and Knight, 2005) were calculated for each land-use category, and compared between land-use categories using ANOVA and post-hoc Tukey HSD tests. Heatmaps of relative OTU abundance and distribution among sites were generated using phyloSeq (McMurdie and Holmes, 2013), for the 1000 terrestrial invertebrate OTUs with the highest proportional abundances across sites.

Biodiversity estimates were calculated for each sample based on the overall terrestrial invertebrate communities, and for each of the main invertebrate groups detected, in such a way that all terrestrial invertebrate OTUs were represented: (1) the dominant insect orders detected (Coleoptera, Diptera, Hymenoptera, Lepidoptera, and Hemiptera, each represented by >150 OTUs); a further 18 insect orders represented by 1 to 36 OTUs were considered as a single pooled group (‘other insects’); (2) non-insect arthropod groups (non-mite Arachnida, mites, Collembola, Malacostraca, myriapods); and (3) non-arthropod phyla (Annelida, Mollusca, Nematoda, Platyhelminthes, Rotifera, and Tardigrada). Because many OTUs were only found in a single site, biodiversity estimates were also calculated with these OTUs excluded, to check whether this affected the results. Species richness and effective species numbers (exponential of Shannon entropy; Jost 2006), were calculated for each invertebrate group using the R packages vegan v2.4–3 (Oksanen et al., 2017) and vegetarian v1.2 (Charney, 2012) respectively. To calculate rarity, a weighting factor (w) was determined for each OTU as the reciprocal of its occurrence across all samples (regardless of land use), so that w = 1 for OTUs that occur in only in a single sample, and w approaches zero for OTUs that occur in many samples. For each sample, values of w were then summed for all OTUs occurring in that sample. In other words, rarity represents the number of OTUs per sample adjusted for their occurrence across all samples (Kier et al., 2009; Kier and Barthlott, 2001).

To calculate phylogenetic diversity, phylogenetic rarity, and mean pairwise distance, OTU sequences were aligned using MAFFT v7 (Katoh and Standley, 2013), and phylogenetic trees constructed. Initially, phylogenetic trees were constructed separately for each phylum using both FastTree 2 (Price et al., 2010) and RAxML v8 (Stamatakis, 2014), and for the overall invertebrates using FastTree 2 (construction of the overall invertebrates tree using RAxML failed). As phylum-level trees based on each method and the overall tree pruned to each phylum yielded similar results, the overall tree was used for estimation of phylogenetic biodiversity metrics per sample and taxonomic group. Phylogenetic diversity, in the form of total branch length per sample, and mean pairwise distance were calculated for each taxonomic group using the R package Picante (Kembel et al., 2010). Phylogenetic rarity, in the form of the branch length unique to each sample (based on occurrences across all samples), was calculated for each taxonomic group and sample according to Rosauer et al. (2009) using the R function phylo.endemism (Niperess, 2010). In addition, standardised effect size values were calculated for each of the phylogenetic metrics, by comparing observed values per site to a null distribution generated by 999 randomisations of the data using a regional null model (Kembel et al., 2010; Miller et al., 2017).

ANOVA was used to test for significant differences among mean biodiversity values between land-uses, for overall invertebrate communities and for each of the taxonomic groups, based on each of the biodiversity metrics. We considered land use as an unordered categorical factor in these tests, because we had no a priori expectation about the relative intensity or impact of all five land uses. Any statistically significant ANOVA tests were followed with post hoc two-sided Tukey HSD tests to identify significant pairwise differences among land-use categories. Subsequently, based on our observed rank order of land uses, we derived a numeric rank of 1 to 5 in the order natural forest > planted forest > low producing grassland > high producing grassland > perennial cropland. We refer to this numeric rank as derived land-use rank (DLUR in tables), to make clear that it is derived from our observed results, rather than on any a priori hypothesis as to which land uses might be considered more intense than others. We tested whether this provided the same conclusions as treating land use as a categorical factor for each metric and taxonomic group. We also included DLUR in a further ANOVA test for biodiversity and taxonomic group differences, to test whether different taxonomic groups showed the same patterns.

We also investigated whether environmental covariates might explain biodiversity trends of overall soil invertebrate communities. To do so, we carried out ANOVA tests for effects of spatial variables (latitude and altitude) plus land-use category effects on overall biodiversity estimates for each metric. In addition, we generated a PCA based on spatial and soil chemistry variables. We then tested whether the most important PCA components, plus land-use category, had significant effects on overall biodiversity estimates for each metric.

To investigate whether the biodiversity metrics differed in their sensitivity to land use, non-parametric bootstrapping stratified by taxonomic group with 999 replicates was used to estimate the proportion of variance attributable to land-use effects with 95% confidence intervals, across the set of taxonomic groups and metrics for which significant land-use differences were detected by ANOVA. These results were plotted as a histogram and compared between metrics using a non-parametric Kruskal-Wallis test.

Appendix 1

Appendix 1—figure 1
Download asset Open asset
Taxonomic composition of invertebrate OTUs and sequences.

Phylum and class-level taxonomic composition of terrestrial invertebrate OTUs detected in soil samples from 75 sites distributed across five land-use categories.

Appendix 1—figure 2
Download asset Open asset
A phylogeny of terrestrial invertebrate COI OTU sequences detected in soil samples from 75 sites distributed across five land-use categories.
Appendix 1—figure 3
Download asset Open asset
Biodiversity estimates for overall soil invertebrate communities detected in different land-use categories, with species detected in a single site excluded.

Diamonds and whiskers represent mean values ± standard errors, with individual data points represented by circles. ANOVA test statistics and trend splines are shown for cases with statistically significant biodiversity differences among land-use categories, with letters indicating differences between land-use categories detected by post-hoc Tukey HSD tests.

Appendix 1—figure 4
Download asset Open asset
Biodiversity estimates for soil arthropod groups in different land-use categories, with species detected in a single site excluded.

‘Other insects’ consists of all insect orders other than Coleoptera, Diptera, Hemiptera, Hymenoptera, and Lepidoptera. ‘Non-mites’ consist of Araneae, Opiliones, and Pseudoscorpiones. Diamonds and whiskers represent mean values ± standard errors, with individual data points represented by circles. ANOVA test statistics and trend splines are shown for cases with statistically significant biodiversity differences among land-use categories, with letters indicating differences between land-use categories detected by post-hoc Tukey HSD tests.

Appendix 1—figure 5
Download asset Open asset
Biodiversity estimates for non-arthropod soil invertebrate phyla in different land-use categories, with species detected in a single site excluded.

Diamonds and whiskers represent mean values ± standard errors, with individual data points represented by circles. ANOVA test statistics and trend splines are shown for cases with statistically significant biodiversity differences among land-use categories, with letters indicating differences between land-use categories detected by post-hoc Tukey HSD tests.

Appendix 1—figure 6
Download asset Open asset
Patterns of environmental variables across five land-use categories.
Appendix 1—figure 7
Download asset Open asset
Richness correlations between different taxonomic groups.

Numbers indicate Pearson correlation coefficients. Ellipse shape and colour represent the magnitude of correlations with p-values≤0.05.

Appendix 1—figure 8
Download asset Open asset
Effective species number correlations between different taxonomic groups.

Numbers indicate Pearson correlation coefficients. Ellipse shape and colour represent the magnitude of correlations with p-values≤0.05.

Appendix 1—figure 9
Download asset Open asset
Rarity correlations between different taxonomic groups.

Numbers indicate Pearson correlation coefficients. Ellipse shape and colour represent the magnitude of correlations with p-values≤0.05.

Appendix 1—figure 10
Download asset Open asset
Phylogenetic diversity correlations between different taxonomic groups.

Numbers indicate Pearson correlation coefficients. Ellipse shape and colour represent the magnitude of correlations with p-values≤0.05.

Appendix 1—figure 11
Download asset Open asset
Phylogenetic rarity correlations between different taxonomic groups.

Numbers indicate Pearson correlation coefficients. Ellipse shape and colour represent the magnitude of correlations with p-values≤0.05.

Appendix 1—figure 12
Download asset Open asset
Mean pairwise distance correlations between different taxonomic groups.

Numbers indicate Pearson correlation coefficients. Ellipse shape and colour represent the magnitude of correlations with p-values≤0.05.

Appendix 1—table 1
Results of ANOVA tests for significant derived land-use rank (DLUR) trends for overall soil invertebrate communities and each biodiversity metric.
MetricTermDfSum Sq.Mean Sq.F stat.R2P
RichnessDLUR133428.0733428.078.180.110.006
Residuals67273678.564084.750.89
Effective SpeciesDLUR11163.401163.404.660.070.034
Residuals6716728.67249.680.93
RarityDLUR125771.7425771.7431.940.32<0.001
Residuals6754061.04806.880.68
Phylogenetic DiversityDLUR11234.541234.5411.170.140.001
Residuals677404.70110.520.86
Phylogenetic RarityDLUR1311.22311.2231.710.32<0.001
Residuals67657.649.820.68
Mean Pairwise DistanceDLUR10.000.000.990.010.324
Residuals670.110.000.99
Phylogenetic Diversity SESDLUR125.6325.635.620.080.021
Residuals67305.464.560.92
Phylogenetic Rarity SESDLUR1549.79549.7948.950.42<0.001
Residuals67752.5611.230.58
Mean Pairwise Distance SESDLUR111.1611.163.280.050.075
Residuals67228.003.400.95
Appendix 1—table 2
Results of mixed-model ANOVA tests for derived land-use rank (DLUR), land-use category (LCAT), and taxonomic group differences and interactions for each biodiversity metric.
MetricTermDfSum sq.Mean sq.F stat.P
RichnessDLUR1277.45277.457.740.007
LCAT3205.3568.451.910.137
Group1615737.52983.6027.43<0.001
DLUR:Group161014.9163.431.770.031
LCAT:Group482425.3750.531.410.037
Effective SpeciesDLUR166.4566.459.280.003
LCAT336.9312.311.720.173
Group163000.93187.5626.19<0.001
DLUR:Group16155.719.731.360.155
LCAT:Group48293.706.120.850.749
RarityDLUR1222.86222.8624.71<0.001
LCAT317.235.740.640.594
Group163082.62192.6621.36<0.001
DLUR:Group16421.7926.362.92<0.001
LCAT:Group48318.826.640.740.908
Phylogenetic DiversityDLUR112.9712.9712.830.001
LCAT34.731.581.560.208
Group16520.0332.5032.14<0.001
DLUR:Group1662.753.923.88<0.001
LCAT:Group4888.351.841.82<0.001
Phylogenetic RarityDLUR14.144.1431.77<0.001
LCAT30.280.090.720.543
Group1638.742.4218.56<0.001
DLUR:Group1610.120.634.85<0.001
LCAT:Group486.940.141.110.288
Mean Pairwise DistanceDLUR10.210.212.870.096
LCAT30.210.070.950.421
Group1619.111.1916.40<0.001
DLUR:Group162.930.182.510.001
LCAT:Group484.200.091.200.169
Appendix 1—table 3
Results of ANOVA tests for effects of spatial attributes (latitude and altitude) and land-use category on overall invertebrate community biodiversity metrics.
MetricTermDfSum Sq.Mean Sq.F stat.R2P
RichnessLatitude179.4979.490.020.0000.882
Altitude172529.2872529.2820.190.236<0.001
Land use411761.722940.430.820.0380.518
Residuals62222736.153592.520.725
Effective SpeciesLatitude1283.47283.471.190.0160.280
Altitude12241.172241.179.410.1250.003
Land use4597.76149.440.630.0330.645
Residuals6214769.66238.220.825
RarityLatitude1465.12465.120.600.0060.443
Altitude117387.7417387.7422.330.218<0.001
Land use413699.463424.874.400.1720.003
Residuals6248280.46778.720.605
Phylogenetic DiversityLatitude10.750.750.010.0000.933
Altitude11740.881740.8816.780.202<0.001
Land use4464.28116.071.120.0540.356
Residuals626433.33103.760.745
Phylogenetic RarityLatitude12.962.960.300.0030.586
Altitude1164.59164.5916.680.170<0.001
Land use4189.5747.394.800.1960.002
Residuals62611.749.870.631
Mean Pairwise DistanceLatitude10.000.000.390.0060.536
Altitude10.000.000.680.0100.411
Land use40.010.001.410.0820.241
Residuals620.100.000.902
Appendix 1—table 4
Results of ANOVA tests for effects of the first three components of a PCA on environmental covariates, plus land-use category, on overall invertebrate community biodiversity metrics.

A PCA was carried out on spatial (latitude and altitude) and soil chemistry variables (pH, C, N, C:N ratio, Olsen P, Total P, Ca, Mg, K, Na, cation exchange capacity, base saturation), of which the first three components explained 70.25% of variation.

MetricTermDfSum Sq.Mean Sq.F stat.R2P
RichnessPC1112142.5712142.573.250.0400.076
PC2126135.3626135.367.000.0850.010
PC31414.99414.990.110.0010.740
Land use440528.7010132.182.710.1320.038
Residuals61227885.013735.820.742
Effective SpeciesPC11497.15497.152.020.0280.161
PC21618.84618.842.510.0350.118
PC3111.1911.190.050.0010.832
Land use41725.81431.451.750.0960.151
Residuals6115039.07246.540.841
RarityPC1111905.8811905.8815.250.149<0.001
PC217487.417487.419.590.0940.003
PC31233.36233.360.300.0030.587
Land use412569.113142.284.020.1570.006
Residuals6147637.01780.930.597
Phylogenetic DiversityPC11503.99503.994.790.0580.032
PC21812.16812.167.720.0940.007
PC3110.9810.980.100.0010.748
Land use4897.84224.462.130.1040.087
Residuals616414.27105.150.742
Phylogenetic RarityPC11147.45147.4515.350.152<0.001
PC2198.3898.3810.240.1020.002
PC3110.3410.341.080.0110.304
Land use4126.5431.633.290.1310.017
Residuals61586.159.610.605
Mean Pairwise DistancePC110.0020.0021.4320.0210.236
PC210.0010.0010.4780.0070.492
PC310.0020.0020.9960.0150.322
Land use40.0060.0010.8630.0510.491
Residuals610.1050.0020.906
Appendix 1—table 5
Results of ANOVA tests for significant derived land-use rank (DLUR) trends for each taxonomic group and biodiversity metric.
MetricGroupTermDfSum Sq.Mean Sq.F stat.R2P
RichnessCollembolaDLUR128.93528.9352.6490.0390.108
Residuals65710.11010.9250.961
ColeopteraDLUR1335.849335.8499.9610.1310.002
Residuals662225.21033.7150.869
DipteraDLUR1615.926615.92618.0120.214<0.001
Residuals662256.83934.1950.786
HymenopteraDLUR1293.041293.04119.0240.229<0.001
Residuals64985.82315.4030.771
LepidopteraDLUR1476.785476.78519.3280.227<0.001
Residuals661628.08324.6680.773
HemipteraDLUR126.78826.7882.9200.0440.092
Residuals64587.1669.1740.956
other insectsDLUR174.53174.5318.9680.1230.004
Residuals64531.9098.3110.877
non-mitesDLUR145.41245.4122.9240.0430.092
Residuals651009.45415.5300.957
mitesDLUR10.2810.2810.0080.0000.930
Residuals662359.71935.7531.000
MalacostracaDLUR10.0440.0440.0460.0010.831
Residuals3432.7060.9620.999
myriapodsDLUR11.5181.5180.5550.0170.462
Residuals3287.4532.7330.983
AnnelidaDLUR167.63767.6374.4560.0650.039
Residuals64971.39315.1780.935
MolluscaDLUR10.2400.2400.0140.0000.906
Residuals641082.24516.9101.000
NematodaDLUR1110.491110.4910.5590.0080.457
Residuals6713243.798197.6690.992
PlatyhelminthesDLUR10.1750.1750.1950.0040.661
Residuals4943.9820.8980.996
RotiferaDLUR1201.555201.5550.6350.0100.428
Residuals6620954.136317.4870.990
TardigradaDLUR10.1770.1770.4690.0150.499
Residuals3011.3230.3770.985
Effective SpeciesCollembolaDLUR10.5300.5300.2730.0040.603
Residuals65126.2001.9420.996
ColeopteraDLUR113.59513.5953.8660.0550.053
Residuals66232.0673.5160.945
DipteraDLUR144.60244.60210.9670.1420.002
Residuals66268.4114.0670.858
HymenopteraDLUR121.82621.8264.6140.0670.036
Residuals64302.7684.7310.933
LepidopteraDLUR162.71962.71912.3840.1580.001
Residuals66334.2575.0640.842
HemipteraDLUR10.6660.6660.2810.0040.598
Residuals64151.7752.3710.996
other insectsDLUR11.1861.1860.8490.0130.360
Residuals6489.4621.3980.987
non-mitesDLUR16.5396.5391.9970.0300.162
Residuals65212.8053.2740.970
mitesDLUR11.0641.0640.2560.0040.615
Residuals66274.7164.1620.996
MalacostracaDLUR10.0070.0070.0360.0010.852
Residuals346.7540.1990.999
myriapodsDLUR10.6370.6370.6100.0190.441
Residuals3233.4141.0440.981
AnnelidaDLUR115.53015.5309.8020.1330.003
Residuals64101.4001.5840.867
MolluscaDLUR14.3294.3291.0210.0160.316
Residuals64271.4454.2410.984
NematodaDLUR116.68216.6820.5710.0080.452
Residuals671955.70129.1900.992
PlatyhelminthesDLUR10.0940.0940.3400.0070.563
Residuals4913.6220.2780.993
RotiferaDLUR1131.612131.6122.6150.0380.111
Residuals663321.55050.3270.962
TardigradaDLUR10.1740.1740.6300.0210.434
Residuals308.2910.2760.979
RarityCollembolaDLUR12.8912.8911.3310.0200.253
Residuals65141.1732.1720.980
ColeopteraDLUR1263.306263.30625.7860.281<0.001
Residuals66673.93010.2110.719
DipteraDLUR1422.250422.25051.6910.439<0.001
Residuals66539.1398.1690.561
HymenopteraDLUR1102.088102.08819.3990.233<0.001
Residuals64336.8095.2630.767
LepidopteraDLUR1256.685256.68536.2000.354<0.001
Residuals66467.9837.0910.646
HemipteraDLUR116.81616.8167.0480.0990.010
Residuals64152.6952.3860.901
other insectsDLUR141.82841.82814.9030.189<0.001
Residuals64179.6312.8070.811
non-mitesDLUR165.61465.61413.7570.175<0.001
Residuals65310.0204.7700.825
mitesDLUR137.89537.8955.6750.0790.020
Residuals66440.6986.6770.921
MalacostracaDLUR10.1260.1260.2410.0070.627
Residuals3417.7450.5220.993
myriapodsDLUR14.4264.4263.5110.0990.070
Residuals3240.3471.2610.901
AnnelidaDLUR147.96647.96617.0570.210<0.001
Residuals64179.9722.8120.790
MolluscaDLUR115.08115.0812.9920.0450.088
Residuals64322.5355.0400.955
NematodaDLUR1197.691197.6916.0080.0820.017
Residuals672204.66432.9050.918
PlatyhelminthesDLUR10.3330.3330.8540.0170.360
Residuals4919.1120.3900.983
RotiferaDLUR1157.017157.0172.0330.0300.159
Residuals665097.29277.2320.970
TardigradaDLUR10.3040.3041.1970.0380.283
Residuals307.6130.2540.962
Phylogenetic DiversityCollembolaDLUR10.7110.7111.5050.0230.224
Residuals6530.7180.4730.977
ColeopteraDLUR121.95221.9529.6880.1280.003
Residuals66149.5442.2660.872
DipteraDLUR140.10640.10619.3840.227<0.001
Residuals66136.5572.0690.773
HymenopteraDLUR121.84821.84812.7310.1660.001
Residuals64109.8301.7160.834
LepidopteraDLUR142.88042.88020.4420.236<0.001
Residuals66138.4462.0980.764
HemipteraDLUR15.7735.7733.8830.0570.053
Residuals6495.1471.4870.943
other insectsDLUR118.45618.45613.6310.176<0.001
Residuals6486.6521.3540.824
non-mitesDLUR113.61213.6127.0640.0980.010
Residuals65125.2591.9270.902
mitesDLUR17.9757.9752.9260.0420.092
Residuals66179.9242.7260.958
MalacostracaDLUR10.4680.4680.7550.0220.391
Residuals3421.0690.6200.978
myriapodsDLUR10.0020.0020.0050.0000.944
Residuals3215.6130.4881.000
AnnelidaDLUR110.57110.57113.5780.175<0.001
Residuals6449.8250.7790.825
MolluscaDLUR11.0721.0720.3790.0060.540
Residuals64181.1562.8310.994
NematodaDLUR10.0020.0020.0000.0000.985
Residuals67305.6604.5621.000
PlatyhelminthesDLUR11.3141.3142.0570.0400.158
Residuals4931.3050.6390.960
RotiferaDLUR15.4445.4442.6650.0390.107
Residuals66134.8092.0430.961
TardigradaDLUR10.1740.1741.5940.0500.217
Residuals303.2790.1090.950
Phylogenetic RarityCollembolaDLUR10.1650.1654.0700.0590.048
Residuals652.6400.0410.941
ColeopteraDLUR18.1518.15126.0350.283<0.001
Residuals6620.6630.3130.717
DipteraDLUR110.22110.22140.0130.377<0.001
Residuals6616.8590.2550.623
HymenopteraDLUR13.5363.53619.3190.232<0.001
Residuals6411.7130.1830.768
LepidopteraDLUR16.9316.93136.8770.358<0.001
Residuals6612.4050.1880.642
HemipteraDLUR10.4670.4674.0950.0600.047
Residuals647.3060.1140.940
other insectsDLUR12.7522.75215.7550.198<0.001
Residuals6411.1800.1750.802
non-mitesDLUR14.0734.07316.1590.199<0.001
Residuals6516.3840.2520.801
mitesDLUR12.3032.30311.0520.1430.001
Residuals6613.7520.2080.857
MalacostracaDLUR10.0210.0210.3140.0090.579
Residuals342.2360.0660.991
myriapodsDLUR10.2000.2004.7080.1280.038
Residuals321.3590.0420.872
AnnelidaDLUR12.0242.02425.9660.289<0.001
Residuals644.9890.0780.711
MolluscaDLUR11.4181.4183.5700.0530.063
Residuals6425.4240.3970.947
NematodaDLUR11.5741.5744.9430.0690.030
Residuals6721.3280.3180.931
PlatyhelminthesDLUR10.0000.0000.0020.0000.966
Residuals492.8730.0591.000
RotiferaDLUR10.7010.7013.0200.0440.087
Residuals6615.3210.2320.956
TardigradaDLUR10.0950.0953.5430.1060.070
Residuals300.8010.0270.894
Mean Pairwise DistanceCollembolaDLUR10.0280.0281.1930.0180.279
Residuals651.5020.0230.982
ColeopteraDLUR10.0020.0020.0970.0010.757
Residuals661.3380.0200.999
DipteraDLUR10.0050.0050.1270.0020.722
Residuals662.3750.0360.998
HymenopteraDLUR10.4210.4217.0010.0990.010
Residuals643.8500.0600.901
LepidopteraDLUR10.0020.0020.0640.0010.801
Residuals662.3250.0350.999
HemipteraDLUR10.1130.1131.2110.0190.275
Residuals645.9680.0930.981
other insectsDLUR10.0390.0390.6310.0100.430
Residuals643.9860.0620.990
non-mitesDLUR10.2110.2112.8250.0420.098
Residuals654.8470.0750.958
MitesDLUR10.4670.46713.8850.174<0.001
Residuals662.2220.0340.826
MalacostracaDLUR10.6890.6891.9570.0540.171
Residuals3411.9680.3520.946
myriapodsDLUR10.1370.1370.5980.0180.445
Residuals327.3180.2290.982
AnnelidaDLUR10.4360.4367.8050.1090.007
Residuals643.5740.0560.891
MolluscaDLUR10.2280.2282.5770.0390.113
Residuals645.6510.0880.961
NematodaDLUR10.0010.0010.1680.0030.683
Residuals670.5510.0080.997
PlatyhelminthesDLUR10.5400.5401.5280.0300.222
Residuals4917.3130.3530.970
RotiferaDLUR10.0070.00714.1780.177<0.001
Residuals660.0330.0010.823
TardigradaDLUR10.1340.1341.3770.0440.250
Residuals302.9100.0970.956

Data availability

Sequence data, metadata, processed data, bioinformatic processing and analysis code used to generate the results in the manuscript (with one exception, detailed below) are deposited in the Manaaki Whenua-Landcare Research DataStore, accessible at: https://doi.org/10.7931/w3j3-5v40. Our sample sites include many Māori and/or privately-owned locations. We have have removed site location details from our metadata out of respect for concerns of Māori and other landowners. The removal of site location details precludes recreation of the map of sample site locations (Figure 6) and analyses of latitude effects, but otherwise has no impact on our results.

The following data sets were generated

References

    1. Folmer O
    2. Black M
    3. Hoeh W
    4. Lutz R
    5. Vrijenhoek R
    (1994)
    DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates
    Molecular Marine Biology and Biotechnology 3:294–299.
  36. Book
    1. Hurst JM
    2. Allen RB
    (2007)
    A Permanent Plot Method for Monitoring Indigenous Forests: Field Protocols
    Manaaki Whenua-Landcare Research.
  63. Book
    1. Pawson SM
    2. Brockerhoff EG
    3. Meenken ED
    4. Didham RK
    (2009) Non-native plantation forests as alternative habitat for native forest beetles in a heavily modified landscape
    In: Brockerhoff E. G, Jactel H, Parrotta J. A, Quine C. P, Sayer J, Hawksworth D. L, editors. Plantation Forests and Biodiversity: Oxymoron or Opportunity?. Springer. pp. 203–224.
    https://doi.org/10.1007/s10531-008-9380-x
    1. Pawson SM
    2. Ecroyd CE
    3. Seaton R
    4. Shaw WB
    5. Brockerhoff EG
    (2010)
    New Zealand’s exotic plantation forests as habitats for threatened indigenous species
    New Zealand Journal of Ecology 34:342–355.
  69. Software
    1. R Development Core Team
    (2016) R: a language and environment for statistical computing
    R Foundation for Statistical Computing, Vienna, Austria.
    https://www.r-project.org/

Article and author information

Author details

  1. Andrew Dopheide

    Manaaki Whenua - Landcare Research, Auckland, New Zealand
    Contribution
    Conceptualization, Formal analysis, Software, Bioinformatics, Data curation, Visualization, Writing - original draft, Writing - reviewing and editing
    For correspondence
    dopheidea@landcareresearch.co.nz
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1554-9832

  2. Andreas Makiola

    Bio-Protection Research Centre, Lincoln University, Lincoln, New Zealand
    Contribution
    Investigation, Fieldwork, Writing - reviewing and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9611-9238

  3. Kate H Orwin

    Manaaki Whenua – Landcare Research, Lincoln, New Zealand
    Contribution
    Conceptualization, Funding acquisition, Writing - reviewing and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0906-728X

  4. Robert J Holdaway

    Manaaki Whenua – Landcare Research, Lincoln, New Zealand
    Contribution
    Conceptualization, Funding acquisition, Investigation, Fieldwork, Writing - reviewing and editing
    Competing interests
    No competing interests declared

  5. Jamie R Wood

    Manaaki Whenua – Landcare Research, Lincoln, New Zealand
    Contribution
    Funding acquisition, Investigation, Molecular analysis, Writing - reviewing and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8008-6083

  6. Ian A Dickie

    Bio-Protection Research Centre, School of Biological Sciences, University of Canterbury, Christchurch, New Zealand
    Contribution
    Conceptualization, Funding acquisition, Investigation, Fieldwork, Formal analysis, Software, Bioinformatics, Writing - reviewing and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2740-2128

Funding

Ministry of Business, Innovation and Employment (MBIE Contract No. C09X1411)

  • Kate H Orwin
  • Robert J Holdaway
  • Jamie R Wood
  • Ian A Dickie

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank landowners, and the New Zealand Department of Conservation, for allowing access to sampling sites, Manaaki Whenua – Landcare Research staff, especially Chris Morse, Larry Burrows, Thomas Easdale, Karen Boot, Nicola Bolstridge, Hamish Maul, and Carina Davis, for their extensive support in the field and laboratory, and peer reviewers for their helpful insights. This research was funded by the New Zealand Ministry of Business, Innovation and Employment (MBIE Contract No. C09X1411) via Manaaki Whenua – Landcare Research.

Copyright

© 2020, Dopheide et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,219
    views
  • 356
    downloads
  • 24
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Citations by DOI

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Andrew Dopheide
  2. Andreas Makiola
  3. Kate H Orwin
  4. Robert J Holdaway
  5. Jamie R Wood
  6. Ian A Dickie
(2020)
Rarity is a more reliable indicator of land-use impacts on soil invertebrate communities than other diversity metrics
eLife 9:e52787.
https://doi.org/10.7554/eLife.52787

Share this article

https://doi.org/10.7554/eLife.52787

Categories and tags

Further reading

  1. Listen to Andrew Dopheide discuss the diversity of insects in soil

    Podcast

    Fewer species of insects live in the soil of agricultural land than in natural environments.