Evaluating restoration trajectories using DNA metabarcoding of ground-dwelling and airborne invertebrates and associated plant communities.

Invertebrates are important for restoration processes as they are key drivers of many landscape-scale ecosystem functions; including pollination, nutrient cycling and soil formation. However, invertebrates are often overlooked in restoration monitoring because they are highly diverse, poorly described, and time-consuming to survey, and require increasingly scarce taxonomic expertise to enable identification. DNA metabarcoding is a relatively new tool for rapid survey that is able to address some of these concerns, and provide information about the taxa with which invertebrates are interacting via food webs and habitat. Here, we evaluate how invertebrate communities may be used to determine ecosystem trajectories during restoration. We collected ground-dwelling and airborne invertebrates across chronosequences of mine-site restoration in three ecologically disparate locations in Western Australia and identified invertebrate and plant communities using DNA metabarcoding. Ground-dwelling invertebrates showed the clearest restoration signals, with communities becoming more similar to reference communities over time. These patterns were weaker in airborne invertebrates, which have higher dispersal abilities and therefore less local fidelity to environmental conditions. Although we detected directional changes in community composition indicative of invertebrate recovery, patterns observed were inconsistent between study locations. The inclusion of plant assays allowed identification of plant species, as well as potential food sources and habitat. We demonstrate that DNA metabarcoding of invertebrate communities can be used to evaluate restoration trajectories. Testing and incorporating new monitoring techniques such as DNA metabarcoding is critical to improving restoration outcomes.

INTRODUCTION

Fauna are often overlooked in restoration monitoring in favour of vegetation (Borges et al., 2021; Cross et al., 2019; Ruiz‐jaen & Aide, 2005), with the general assumption that animals will naturally recolonize an area with the return of plant communities (Palmer et al., 1997). However, this is not always the case (Cristescu et al., 2013), and understanding the recovery of animals is important because they play a vital role in many ecosystem functions, including pedogenesis, seed dispersal, pollination and nutrient cycling (Bronstein et al., 2006; Catterall, 2018; Hunter, 2001; Ness et al., 2004; Sekercioglu, 2006). Recently, greater attention has been paid to fauna to both assess and facilitate ecological restoration (Catterall, 2018; Cross et al., 2020; Majer, 2009).

Invertebrates are of particular interest as they have long been used as indicators of ecosystem recovery in both aquatic and terrestrial systems (Andersen et al., 2002; Andersen & Sparling, 1997; Folgarait, 1998; Majer, 2009). They are sensitive to disturbances and are essential for ecosystem function (Folgarait, 1998; Rosenberg et al., 1986), while being useful indicators because of their abundance, ease of capture and high diversity, particularly of trophic types (Gaston, 1991). Because studies tend to investigate particular groups of arthropods, responses to restoration are mixed, depending on the target taxa (Cristescu et al., 2012). Some of the variation in responses to restoration among different arthropod classes may be attributed to dispersal ability. For example, beetles with high dispersal abilities are able to recolonize more quickly than millipedes in a regenerating forest (Magura et al., 2015). However, it is unknown whether ground‐dwelling invertebrates show recovery trajectories better than airborne invertebrates (Moir et al., 2005), or if patterns are consistent across multiple locations, as most studies have been limited to a chronosequence of restoration sites in a single ecosystem (e.g., Fernandes et al., 2019; Magura et al., 2015). Several recent studies on the recovery of various taxa, including soil microbial communities (van der Heyde, Bunce, Wardell‐Johnson, et al., 2020) and vertebrates (van der Heyde, Bateman, Bunce et al., 2021) show that recovery patterns are complex and vary among locations, ecosystems and taxa. Therefore, we regard the inclusion of multiple study locations as being an important feature in study design (see also Catterall et al., 2004).

Despite being excellent indicators of ecosystem change, the high diversity within invertebrate communities makes it difficult to identify captured invertebrate specimens. Thus, many expert person‐hours from multiple taxonomists specializing in different invertebrate taxa are often required (Majer, 1983). This process is costly and time consuming, and is dependent on taxonomic expertise that is dwindling worldwide (Majer et al., 2013; Pearson et al., 2011). Additionally, many invertebrate taxa are cryptic (Smith et al., 2005) or have yet to be identified, especially in Australia with its high degree of endemism (Austin et al., 2004; Rix et al., 2015) and where as much as 75% of arthropod diversity is undescribed (Austin et al., 2004; Yeates et al., 2003). Consequently, most studies examining invertebrate responses to restoration have targeted particular taxa, either because they have been previously shown to be good bioindicators (Andersen et al., 2002), or they are threatened and therefore of regulatory and conservation value (i.e., Lepidoptera) (Majer, 2009).

Some of the difficulties associated with invertebrate monitoring can be reduced using DNA metabarcoding to provide community composition profiles. This process uses high‐throughput sequencing of small barcoding regions of the genome to determine invertebrate diversity (Beng et al., 2016; Ji et al., 2013; Yu et al., 2012). Compared to morphological identification, where each specimen has to be identified individually, DNA metabarcoding has been shown to be accurate, reliable, and faster than conventional morphological methods (Beng et al., 2016; Ji et al., 2013). As an added benefit, the sequencing data can be readily stored and analysed by a third party, such as regulators (Fernandes et al., 2018). Although abundance estimates using DNA metabarcoding are often skewed by primer bias (Elbrecht & Leese, 2015), and/or DNA extraction method (Majaneva et al., 2018), presence/absence data has been used to demonstrate arthropod responses to restoration post mining (Fernandes et al., 2019) and to land‐use change (Beng et al., 2016).

One of the advantages of DNA metabarcoding over morphology based approaches is its ability to detect invertebrate diversity and composition and also provide information on plant species that they are using as forage and habitat (Jurado‐Rivera et al., 2009; Pornon et al., 2016). In the case of arthropods, previous studies suggest that DNA from arthropod samples should be able to identify which plant species the pollinators have visited (Pornon et al., 2016) and which plant species they have consumed (Jurado‐Rivera et al., 2009). However, these studies have hitherto not been undertaken in a restoration context, so the utility of such approaches for restoration monitoring is unknown. Assessing these communities can demonstrate interactions between invertebrates and plants during restoration programs. However, since the invertebrates may carry plant DNA from outside the restoration area (van der Heyde, Bunce, Dixon, et al., 2020), they may not necessarily have high fidelity to local conditions.

Our earlier work has explored the use of DNA metabarcoding of ground‐dwelling invertebrates to monitor mine site restoration (Fernandes et al., 2019). However, that study used a single reference site per mine and the results were spatially correlated, in that older sites were also closest to the reference sites. Here we used two spatially separated reference sites per mine to reduce such a correlation bias, and two trap types that capture ground‐dwelling and airborne invertebrates. We also use study sites at three locations in different climates and ecosystems. This study evaluates the applicability of DNA metabarcoding of invertebrates to evaluate restoration trajectories (convergence to reference communities) in restored sites. We have three hypotheses: (i) Ground‐dwelling invertebrates will show recovery trajectories more effectively than airborne invertebrates because they more clearly reflect local environmental conditions. (ii) Trajectories of recovery vary with location and environment. (iii) Metabarcoding plant sequences from bulk invertebrate samples provides plant species occurrence and habitat information.

MATERIALS AND METHODS

Study sites

Restoration and reference sites were sampled from three locations up to 1000 km apart in Western Australia, namely: Swan Coastal Plain (SCP); Jarrah Forest (JF); and Pilbara (PB). There was consistency in restoration approaches, soil type, climate and site aspect of the sites within each location. At each location, sites of different restoration age were sampled along with two spatially separated reference sites (Figure 1, see Figure S1 for maps). At all three locations, we sampled at least two sites less than nine years old (Young), and at least two sites older than nine years (Older). These sites are previously described in van der Heyde, Bunce, Wardell‐Johnson, et al. (2020), and briefly below. At all locations, two reference sites were selected on the basis of the following criteria: similarity to ecosystems that are the target of restoration efforts, proximity to restoration sites, topographical similarity, and spatially separate from each other to account for variation in reference communities. All restored sites were established to address requirements for site rehabilitation post‐mining rather than for our study objectives. As a result, there is a lack of site replication. Despite this, we suggest that our conclusions provide meaningful insight into the return of invertebrate communities following restoration, and represent a case study on the application of DNA metabarcoding to restoration monitoring.

FIGURE 1

Chronosequences of mining restoration where invertebrate samples were collected. Restoration sites shown with the number of years restoration from 1 to 22 years. Reference sites shown below. JF, Jarrah Forest; PB, Pilbara; SCP, Swan Coastal Plain

The coastal plain site (SCP) has a warm‐summer Mediterranean climate with mild cool wet winters; mean minimum temperature 12.8°C, mean maximum 24.7°C, and with 757 mm mean annual rainfall (Australian Bureau of Meteorology). This location is part of the broader region of south‐western Australia, a globally recognized biodiversity hotspot (Myers et al., 2007). The mine is located on the silicaceous Bassendean dunes, with high acidity and low water‐holding capacity (Dodd & Heddle, 1989; McArthur, 1991). The ecosystem is referred to as Banksia woodland after the dominant tree species, Protaceae Banksia attenuata and B. menziesii. Other trees include less dominant Myrtaceae Eucalyptus todtiana and Loranthaceae Nuytsia floribunda. The understory consists of woody species of Myrtaceae, Ericaceae, Proteaceae, and nonwoody species in Asparagaceae, Stylidiaceae, Cyperaceae, and Haemodoraceae (Trudgen, 1977). In October 2018, we sampled eight sites at a Hanson Construction Materials sand quarry in Lexia (31.76°S, 115.95°E), with two reference sites and restoration sites 1, 3, 7, 11,14, and 22 years old. The sites have been restored with the aim of returning mined areas to the surrounding native Banksia woodlands. All restoration was done by Hanson and previous mine owners and reflect best practice in mine restoration through direct transfer of fresh topsoil ripping, and seeding with native plant species. Plant species richness and density tended to be higher in restoration than reference sites, and percent cover has increased with restoration age and is highest in reference sites (Benigno et al., 2013).

The second location in the Jarrah (Eucalyptus marginata) forest (JF) is also part of the Southwest Australia Global Biodiversity hotspot (Myers et al., 2007) and has a similar hot‐summer Mediterranean climate; mean minimum temperature of 8.6°C, mean maximum of 23.7°C, and 668.9 mm annual mean rainfall (Australian Bureau of Meteorology). The lateritic soils are nutrient poor and high in gravel, with surfaces rich in iron and aluminum (McArthur, 1991). The vegetation is dominated by E. marginata, with E. patens, and E. wandoo also being common. The understory consists of sclerophyllous shrubs from several families including Asparagaceae, Fabaceae, Asteraceae, Proteaceae, Dasypogonaceae, and Myrtaceae (Havel, 1975). In October 2018, we sampled six sites from the bauxite mine which is now run by South32 (32.96°S, 116.48°E); two reference sites and restoration sites 2, 6, 11, and 20 years old. All restoration was undertaken by South32 or the previous mine owners. After mining the landscape was shaped using waste material and gravel. Fresh topsoil was directly transferred from newly mined areas to the restoration area and supplemented with stockpiled topsoil as needed. The sites were then ripped, seeded with over 100 native species, recalcitrant plants (mostly grasses) were planted, and a one‐time treatment of superphosphate was applied (Data from South32). Reference and restoration sites are dominated by Myrtaceae and Fabaceae species. Total cover has increased with age of restoration to similar cover values of reference sites (Data from South32).

The third location, the Pilbara (PB), is in north‐western Australia. The Pilbara has a hot, arid climate, with most rainfall occurring in summer, and associated with irregular cyclonic activity (McKenzie et al., 2009) causing unpredictable flooding. Temperatures have a mean minimum of 15°C and mean maximum of 30.6°C, with 263.8 mm mean rainfall (Australian Bureau of Meteorology). Soils are acidic stony loams with low fertility, which support open woodlands of snappy gum (E. leucophloia) over hummock grasses (Poaceae Triodia wiseana, T. basedowii, T. lanigera) and low Fabaceae Acacia shrubs (McKenzie et al., 2009). The harsh climate, large variation in yearly rainfall, and low soil fertility, result in low productivity when compared to the other study sites. The Pilbara is a significant mining region and accounts for 39% of global iron ore production (Government of Western Australia, 2019). We sampled six sites at a BHP iron ore mine (22.84°S, 118.95°E) in September 2018, with two reference sites and restoration sites 4, 7, 11, and 15 years old. Restoration was conducted by the mine owners; landscapes were reformed and stockpiled topsoil (average age 10 years) was applied and then ripped. Restoration areas tended to have higher coverage of woody shrubs (Acacia), while reference sites and older restoration areas have more hummock grasses (Triodia) and a sparse shrub stratum. Restoration areas also had invasive species such as buffel grass (Poaceae Cenchrus ciliaris) and kapok bush (Amaranthaceae Aerva javanica), which were absent in reference sites (Data from BHP).

Sample collection

At each site we collected 10 invertebrate samples, (five from vane traps and five from pitfall traps), there were eight SCP sites and six sites from PB and JF, for a total of 200 samples. Each vane trap sample included the contents of a yellow and blue vane trap containing 150 ml of ethylene glycol, with traps remaining on site for seven days. Vane traps of different colours tend to capture different invertebrates (Hall, 2018), and for this study the vane traps were pooled to provide a more complete view of the airborne invertebrate community. Each pitfall trap sample included the contents of four pitfall traps (4 cm diameter, 12 cm deep with ethylene glycol as a capture fluid), and was also left in the field for 7 days. For each sample point, pitfall traps were spaced 10 m apart in a square around the vane traps in the centre.

Sample processing

For DNA extraction, we first rinsed off the ethylene glycol with deionized water using 20‐µm sieves that were sterilized in bleach and under UV light between every sample and visible plant material was removed. We used two legs of all specimens larger than a bee (~15 mm) and the whole body for smaller specimens to minimize the effect of body mass on the sequence abundance (Elbrecht et al., 2017; Ji et al., 2013). Samples were then homogenized using a TissueLyser (Qiagen) for 2 min in 30 s increments at 30/s in 50 ml falcon tubes with four steel balls (4 mm diameter). Then, 400 μl of the homogenate was digested overnight and the DNA extracted using the DNeasy Blood and Tissue kit (Qiagen) on the QiaCube Connect automated platform (Qiagen). The final elution volume was 200 μl, and extraction controls (blanks) were carried out for every set of extractions. Quantitative PCR (qPCR) was done on neat extracts and a 1/10 dilution to see if samples exhibited inhibition, and to determine optimal DNA input for PCR for each sample to maximize input relative to any inhibitors (Murray et al., 2015). Two assays were used in this study to target invertebrate and plant diversity. The invertebrate assay used the primers fwhF2/fwhR2n (Vamos et al., 2017) to amplify a 205 bp section of the cytochrome c oxidase I (COI) region. For plants we used the trnlc/h primers (Taberlet et al., 2007) which targets the chloroplast trnL (UAA) intron.

The qPCRs were run on a StepOne Plus (Applied BioSystems) real‐time qPCR instrument with the following conditions: 5 min at 95°C, 40 cycles of 95°C for 30 s, 30 s at the annealing temperature (50°C for invertebrates, 52°C for plants) and 45 s at 72°C, a melt curve stage of 15 s at 95°C 1 min at 60°C and 15 s at 95°C, ending with 10 min elongation at 72°C. The PCR mix for quantitation had a 25 μl volume and contained: 2.5 mM MgCl2 (Applied Biosystems), 1× PCR Gold buffer (Applied Biosystems), 0.25 mM dNTPs (Astral Scientific), 0.4 mg/ml bovine serum albumin (Fisher Biotec), 0.4 μmol/l forward and reverse primer, 1 U AmpliTaq Gold DNA polymerase (Applied Biosystems) and 0.6 μl of a 1:10,000 solution of SYBR Green dye (Life Technologies) and 2 μl of DNA template. Extraction control and nontemplate controls were included in qPCR assays.

After optimal DNA input was determined by qPCR, each sample was assigned a unique combination of multiplex identifier (MID) tags for each primer assay. These MID tags were incorporated into fusion tagged primers, and none of the primer‐MID tag combinations had been used previously in the laboratory to prevent cross contamination. Fusion PCRs were done in duplicate and to minimize PCR stochasticity, the mixes were prepared in a dedicated clean room before DNA was added. The PCRs were done with the same conditions as the standard qPCRs described above. Samples were then pooled into approximately equimolar concentrations to produce a PCR amplicon library that was size‐selected to remove any primer‐dimer that may have accumulated during fusion PCR. Size selection was performed (150–450 bp) using a PippinPrep 2% ethidium bromide cassette (Sage Science). Libraries were cleaned using a QIAquick PCR Purification Kit (Qiagen) and quantified using Qubit Fluorometric Quantitation (Thermo Fisher Scientific). Single‐end sequencing was performed on the Illumina MiSeq platform using the 300 cycle V2 as per manufacturer's instructions.

DNA sequence analysis

Sequences were demultiplexed, removing the primers and MID tags, using a demultiplex function in the insect package (Wilkinson et al., 2018) on the R 3.5.3 platform (R Core Team, 2018). Further sequence processing was performed in R using the DADA2 package (Callahan et al., 2016) where sequences were quality filtered with a minimum length of 100 bp, maxEE = 2, maxN = 0, and phiX removed. The error rates were estimated for each sequencing library separately using the learnErrors function in DADA2, an algorithm which uses a parametric error model learnt from the sequence data by alternating inference of sample composition and estimations of error rates until they converge on a solution (Callahan et al., 2016). This includes both within read (chimera) and between read errors. The error rates were then used with the core sample inference algorithm to remove sequences likely to be errors and leave amplicon sequence variants (ASV) used to construct a sequence table. These ASVs are equivalent to denoised zero radius operational taxonomic units (ZOTUs) in usearch (Edgar, 2016) and the sequence table produced is essentially a higher‐resolution version of the OTU table produced by other methods. The sequence tables for each library were then merged and chimeras removed. Finally, we used LULU to remove spurious ASVs based on sequence similarity and co‐occurrence patterns (Frøslev et al., 2017) creating a curated ASV table for further analyses. Taxonomy was determined using the Basic Local Alignment Search Tool (blastn) on a high‐performance cluster computer (Pawsey Supercomputing Centre) to search against the online reference nucleotide database GenBank (https://www.ncbi.nlm.nih.gov/genbank/) with a minimum percent ID of 85%, minimum query coverage of 90%, and a maximum 10 hits per ASV. Invertebrate sequences were also searched against and arthropod COI reference sequences extracted from the Barcode of Life Database (BOLD: https://www.barcodeoflife.org), because there are reference sequences that are found uniquely on one of the two databases. MEGAN (Huson et al., 2007) was used to assign taxonomy to each sequence by applying a bit‐score threshold of 205 (min‐score), retaining only the hits within 10% of the best hit (top‐percent) to assign sequences to the lowest common ancestor of the matched species.

Statistical analysis

All statistics were run using R 3.5.3 (R Core Team, 2018). Samples with low sequencing depth were removed and ASVs that were present in the extraction controls were removed from the data set (Figure S2). The invertebrate extraction controls did not amplify, so no ASVs were removed. The plant extraction controls only contained 111 reads, which removed two ASVs from the data set. We selected ASVs in the phylum “Arthropoda” for the invertebrate assay and “Plantae” for the plant assay according to the classification on MEGAN. Copy numbers in each sample were filtered to a minimum of 0.5% within sample abundance. We also used a more conservative 1% threshold to avoid ASV inflation and overestimation of ASV richness, a common issue in COI metabarcoding (Andújar et al., 2021), but this made no difference to observed patterns (Figures S3 and S4) and we have therefore presented only the results for the more relaxed filtering parameter. We verified there was no correlation between sequencing depth and ASV richness using a Pearson correlation test, before continuing. Read counts were transformed to presence/absence to reduce the effects of biases (Elbrecht & Leese, 2015; Majaneva et al., 2018). Spatial autocorrelation was tested using the Mantel test in the ade4 package in R (Mantel, 1967). Where there is significant spatial autocorrelation, this would indicate that distance between samples is an important factor explaining the variation in communities, limiting inferences from our other variables (i.e., restoration age).

Three criteria were examined to determine if communities showed a trajectory of recovery or convergence to the reference community. First, community composition should be different between younger restoration, older restoration, and reference sites. This was visualized using nonmetric multidimensional scaling (NMDS), based on a presence/absence ASV table, and with Jaccard similarity because the data were presence/absence. The ordiellipse function from the vegan R package was used to draw ellipses showing the 95% confidence interval of the group (Oksanen et al., 2018). Differences between restoration and reference sites were tested using permutational multivariate analysis of variance (PERMANOVA). Second, establishing a restoration trajectory requires directional change with restored communities expected to become more similar to a reference community. Replicates at each site were pooled and the similarity between each site and both of the reference sites was calculated. This relationship was tested using linear models separately for each assay and location. For this analysis, we used Bray‐Curtis similarity because the pooled sites included the number of replicates with positive detection as a proxy for abundance, so the data was not presence/absence. Keeping the site similarity for both of the references separate implies potential pseudoreplication, but since a site could be more similar to one reference site than another, we felt it relevant to keep them separate. Doing so also allowed us to separate the effect of distance from restoration age using the variance partitioning function (varpart) of the R package vegan (Oksanen et al., 2018) to quantify the variance in community similarity to the reference sites that was explained by each factor (geographic distance vs. restoration age). Third, we expect that the proportion of “reference” ASVs, that is, ASVs that were found in reference sites, would increase over time. This relationship was tested using a simple linear model. For all three, we tested the SCP data with and without the extra two sites (7 and 11 years) to ensure that any comparisons of trajectory between the locations were fair, since the other locations only had four restoration ages while the SCP had six. We recognize that trajectories may not necessarily converge with reference communities and may settle on a new stable state, but we focus on recovery trajectories in this study as it is the restoration target. Finally, to understand the taxa associated with restoration and reference sites, we ran a multipattern analysis for each site using the R package indicspecies, with the sample as the unit of analysis (De Caceres & Legendre, 2009).

RESULTS

In total, 14,780,759 quality‐filtered invertebrate sequences were generated from 196 samples with a mean sequencing depth of 82,934 (±8411 SE) and a minimum of 3000 reads/sample. Out of 5862 initial ASVs, 2635 belonged to the phylum Arthropoda and 951 ASVs remained after abundance filtering. The remaining ASVs were either unidentified or fungi, and only made up 23.7% of the read count. In the plant assay, we generated 13,441,527 filtered plant sequences from 197 samples with a mean sequencing depth of 63,754 (±3870 SE) and a minimum of 5600 sequences/sample. From the initial 511 plant ASVs, 205 remained post filtering and these accounted for 82.4% of the sequences. Overall, there were fewer invertebrate ASVs in the Pilbara (323 ASVs) compared to the Coastal Plain (377 ASVs) or Jarrah (344 ASVs), especially in the pitfall traps where the Pilbara had 17%–28% fewer invertebrate ASVs (Table S1).

Community composition

Invertebrate diversity in the vane traps was dominated by Hymenoptera, Coleoptera, Diptera, Hemiptera, and Lepidoptera. Some of these (Hymenoptera, Coleoptera, and Hemiptera) also made up most of the diversity in the pitfall traps, along with Collembola and Araneae. Collembola were largely absent from the Pilbara, which had more Orthoptera ASVs. The majority (67%) of invertebrate ASVs could not be identified beyond order level. However, 99% of plant ASVs could be identified to family level. Plant diversity in the SCP and JF sites were dominated by Myrtaceae, Fabaceae, Dilleniaceae, and Proteaceae, while in the PB sites, the richest families were Fabaceae, Poaceae and Malvaceae (Figure 2). Because of the poor taxonomic assignments, we confined our considerations to ASVs for our subsequent analyses.

FIGURE 2

Composition of invertebrates (above) and plant (below) communities detected from pitfall and vane traps. Shows the number of ASVs in each order (invertebrates) or family (plants) at all restoration and reference sites. JF, Jarrah Forest; PB, Pilbara; SCP, Swan Coastal Plain

There were significant differences in community composition between younger restoration, older restoration and reference sites in all locations for both trap types and assays (Figure 3, PERMANOVA, p < .05). The Mantel tests showed no significant spatial autocorrelation in the invertebrate communities from pitfall and vane traps (Table 1). For site similarity estimates based on plant sequence data, the spatial autocorrelation between samples was significant only for the coastal plain vane traps (Table 1).

FIGURE 3

NMDS ordinations of invertebrate and plant communities in restoration and reference sites. Ellipses were drawn using Ordiellipse in the vegan R package and indicate 95% confidence interval of the group. Stress values and significance of PERMANOVA tests indicated in the bottom of each facet. JF, Jarrah Forest; PB, Pilbara; SCP, Swan Coastal Plain

TABLE 1

Results of the Mantel test showing the correlation between spatial distances and community dissimilarity. Results for the samples separately, and pooled (sites) are shown

Trap	Assay	Location	r	p‐value
Pitfall	Invertebrate	JF	.236	.093
		PB	–.074	.606
		SCP	–.084	.558
	Plant	JF	.074	.051
		PB	.094	.170
		SCP	–.064	.745
Vane	Invertebrate	JF	.237	.101
		PB	–.074	.553
		SCP	–.084	.517
	Plant	JF	.059	.117
		PB	.042	.324
		SCP	.280	.001

Abbreviations: JF, Jarrah Forest; PB, Pilbara; SCP, Swan Coastal Plain.

Similarity to reference communities

The invertebrate communities showed clear directional changes in the pitfall traps from the coastal plain (p = .001) and the forest (p = .003) (Figure 4). This trajectory was present but weaker in the SCP vane traps (p = .029) and entirely absent in the vane traps of the Jarrah forest. There were no observed directional changes in invertebrate community composition in the Pilbara. The results from the plant communities were different. In the Coastal Plain, there was a significant relationship between similarity to reference communities and age of restoration in the vane traps (p = .026), but not the pitfall traps (p = –.376). The directional changes in forest plant communities were similar to the invertebrate communities, with an increase in similarity over time in the pitfall traps (p = .031) and no relationship in the vane traps (p = .711). Similarly, the plant communities in the Pilbara showed no relationship between restoration age and similarity to reference communities in either the pitfall (p = .659) or the vane traps (0.693) (Figure 4). These results were also reinforced by the variance partitioning, which showed that restoration age explained more of the variance in community similarity except in the Pilbara, where distance to the reference site had greater explanatory power (Table 2).

FIGURE 4

Similarity (Bray‐Curtis) of restoration sites of different ages (years) to communities in both reference sites. (a) Shows the invertebrate communities and (b) the results of the plant assay. Lines indicate linear models with 95% confidence interval shown with shading. p‐values for the linear models shown for each plot. Removing the two extra sites in the SCP (711 years) did not change the relationships or the significance of the models. JF, Jarrah Forest; PB, Pilbara; SCP, Swan Coastal Plain

TABLE 2

Partitioning the variance in community similarity to reference site that can be explained by restoration age and distance to reference site

Assay	Substrate	Location	Restoration	Both	Distance	Residual
Inverterbate	Pitfall	JF	0.839	<0	<0	0.295
		PB	<0	<0	0.177	0.955
		SCP	0.713	<0	0.226	0.146
	Vane	JF	0	<0	0.363	0.698
		PB	<0	<0	0.005	1.052
		SCP	0.217	0.115	<0	0.721
Plant	Pitfall	JF	0.606	<0	<0	0.559
		PB	<0	0.0326	<0	1.183
		SCP	<0	0.0129	<0	1.129
	Vane	JF	<0	<0	0.661	0.478
		PB	<0	<0	0.536	0.598
		SCP	0.854	<0	0.412	0.242

Abbreviations: JF, Jarrah Forest; PB, Pilbara; SCP, Swan Coastal Plain.

Proportion of “reference” associated ASVs

Only the invertebrate communities from the pitfall trap samples from the coastal plain and the forest showed significant increases in the proportion of “reference” ASVs over time. For plant sequences, only vane traps in the coastal plain showed increasing “reference” ASVs over time (Figure 5). Overall, the vane traps had a higher proportion of ASVs that were shared with reference samples than pitfall traps. This was true for both the invertebrate assay (49.7% vs. 22.9% “reference” ASVs) and the plant assay (51.8% vs. 37.0% “reference” ASVs). Between the two reference sites, there was variation in the number of ASVs shared with each other. The pitfall traps in the Pilbara only had 8% invertebrate ASVs shared between the two reference sites. The amount of shared invertebrate ASVs was higher between the coastal plain and forest pitfall traps (28% and 21%, respectively).

FIGURE 5

Proportion of reference associated ASVs in the different sites, separated into (a) invertebrate and (b) plant communities. p‐values indicate the significance of the relationship between the proportion of reference ASVs and age of restoration (years). JF, Jarrah Forest; PB, Pilbara; SCP, Swan Coastal Plain

Multipattern analysis

Across the three locations, there were 66 invertebrate ASVs with significant association (p < .05) with younger restoration (<9 years), older (>9 years), reference sites, or a combination (Table S2). Of these, 34 were assigned to family, 14 to genus, and only three to species level. This includes the ant Iridomyrmex sanguineus, which was associated with younger restoration in the Pilbara and the ant Monomorium rothsteini, associated with reference sites in the Pilbara. Most Coleoptera (8/12) were associated with older restoration or reference sites and 11 of those were from vane trap samples. For the plant assay, there were 35 ASVs with significant association (Table S3), 31 of which were assigned to family, nine to genus, and three to species level. Among these were the family Fabaceae, associated with younger restoration in the Jarrah forest and Pilbara, and the genus Anigozanthus (Haemodoraceae) associated with younger restoration in the coastal plain.

DISCUSSION

Terrestrial invertebrate fauna are key indicators of ecosystem change (Andersen et al., 2002; Majer, 2009; Majer et al., 2007), and in this study, we show that even with limited taxonomic information, DNA metabarcoding of invertebrate samples can be used to rapidly assess complex biological interactions and establish restoration trajectories. These trajectories of community recovery were more evident in older restored sites, and in ground‐dwelling invertebrates with lower dispersal ability than airborne invertebrates. Plant species identified from bulk invertebrates also showed indications of directional changes in community composition.

Different signal strengths from ground‐dwelling and airborne invertebrates

Vane traps did not show the same local fidelity as pitfall traps and, as expected, tend to have weaker indications of community recovery (Figures 3 and 4). Vane traps capture airborne invertebrates, often pollinators (Hall, 2018), and can trap organisms that may come from more than 1.8 km away (Jha & Dick, 2010) while species caught by pitfall traps have more limited catchment areas (Majer, 1980; Ness et al., 2004; Ward et al., 2001). This would also explain the greater proportion of shared taxa in the vane traps compared to the pitfall traps (Figure 5). Beyond the differences in attraction distance of the traps, our results also suggest quicker recolonization of airborne invertebrates as evidenced by the number of “reference” associated taxa is similar to reference sites within a few years (Figure S4, SCP, PB). Variation in dispersal abilities is important, as those with more mobility are able to recolonize areas more quickly (Magura et al., 2015) and from greater distance (Knop et al.,2011). Fortunately, there is no sign of thermophilic or other barriers (Cranmer et al., 2012; Tomlinson et al., 2018) preventing invertebrates from accessing and using restoration sites. Because of their more sedentary nature, ground‐dwelling invertebrates are good indicators of organisms that are probably reproducing in situ, while airborne invertebrates can indicate the forage support and attractiveness of a site.

Our findings indicate that invertebrate communities are demonstrating an ability to recover without intervention subject to suitable source populations being available. This conforms with the “Field of Dreams” hypothesis which states that if suitable habitat can be re‐established, species will colonize it, leading to the restoration of function (Palmer et al., 1997). Again, this is dependent on the presence of source populations with migration ability. In this study, all sites were near remnant vegetation that could act as a taxa pool; in cases of isolated restoration sites, it may be more difficult to evaluate restoration trajectories using invertebrate communities.

Patterns vary among ecosystems

In older restored sites on the coastal plain and forest we recorded significant increases in the proportion of “reference” taxa, which shows a directional change in community composition toward that of the reference community. In contrast, the Pilbara location did not show a similar trajectory of invertebrate community recovery. These trends match those of previous studies using morphological identification, at bauxite mines close to the one reported in this paper (JF) (e.g., Majer, Heterick, et al., 2013), sand mines to the south (e.g., Davieson & Majer, 1983) and north (e.g., Bisevac & Majer, 1999) of our mine (SCP), and a range of Pilbara sites (e.g., Dunlop et al., 1985; Fletcher, 1990).

One possible explanation is variation in the climate and productivity among the three study sites. The Pilbara Region can be classified as a “harsh” environment due its high temperatures, arid climate, poor soils and unpredictable flooding in monsoonal rains (Charles et al., 2013; Sudmeyer, 2016), which limits productivity, and results in open, unvegetated patches, with overall lower percent plant cover than found in coastal woodlands or forests (McKenzie et al., 2009). Dunlop et al. (1985) and Fletcher (1990), observed that ant richness rapidly recovered in young Pilbara rehabilitation, but, similar to our results, the species composition remained different between natural and restored sites. In the Pilbara, the main factors driving compositional turnover in terrestrial fauna are regolith/soil and landform/hydrogeologic, as well as climate (Gibson et al., 2015). All were factors shared between Pilbara restored and reference sites. Here, the structure of the revegetation rapidly came to resemble the structure of the original predominantly grassland habitat (see Figure 1), which is in marked contrast to the other two locations. In that regard, the reference areas may provide conditions that are as unpredictable, unfavourable and unproductive as the areas under restoration; and compared with the other two regions, they are also less rich in species. Thus, recolonization of Pilbara sites may be more stochastic and less influenced by selection pressures than in the coastal plain and Jarrah forest. However, this hypotheses must be treated with caution, as we did not have replication within the different ecosystems, and declines seen at the Pilbara location may be due to a variety of reasons specific to the site (e.g., management practices), or specific to that point in time (e.g., weather events).

Information on plant species

Plant assay metabarcoding of bulk invertebrate samples provided information on local plant species occurrences with directional changes in plant community composition identified (Figures 3 and 4). Changes in plant community composition detected by eDNA metabarcoding are similar to successional changes known to occur at the three study sites. For example, a higher richness of Fabaceae ASVs, many of which are coloniser species, were found in younger restoration sites. A further example, Anigozanthos was significantly associated with younger restoration sites in the coastal plain and was observed in great abundance in the SCP restoration sites. This group is fast growing and rapidly establishes post restoration (Table 3). While some plant DNA may originate from debris falling into traps (this in itself is useful as its still provides information on local plant occurrences), there is an indication that at least some of the plant species detected were likely to have been ingested or otherwise visited by invertebrates. For example, plants in the family Goodeniaceae require insect pollination (Jabaily et al., 2012; Keighery, 1980) and were flowering at two study sites during sample collection (PB and JF). While there are virtually no Goodeniaceae ASVs in the pitfall traps, they are present in most sites in vane traps (PB and JF, Figure 2) suggesting that flying invertebrates visited the flowers of nearby Goodeniaceae species.

TABLE 3

Taxa of interest, based on general observations of the data and indicator species analysis

Taxa of interest	Name	Reason
Melophorus	Australian genus of ant	Associated with younger restoration sites in both the SCP and PB. Species in this genus are known as “sun‐loving” (Andersen et al., 2002) and are often found in restoration sites with large areas of bare ground (Andersen et al., 2003)
Iridomyrmex sanguineus	Northern meat ant	Associated with younger restoration sites in PB. Iridomyrmex species are among the first to colonize revegetated sites (Andersen, 1993). Much like Melophorus, Iridomyrmex species are attracted to areas of bare ground in newly revegetated sites (Andersen et al., 2003)
Hemiptera	Order of sucking insects	Higher richness in younger restoration in the JF pitfall traps. Hemipteran species composition is linked to the presence of host plants, vegetation structure, and soil pH (Orabi et al., 2010). In newly restored sites, there is generally a higher abundance of generalist Hemiptera species, with a slower recolonization of specialist Hemiptera species dependant on vegetation structure reassembly (Moir et al., 2005)
Apidae	Family of bees	Found primarily in the younger restoration sites in JF. Newly restored sites with less ground cover offer optimal nesting area for ground nesting bees (Seitz et al., 2019). Apidae have also been used as important indicators of pollution and stress in ecosystems (Rabea et al., 2010; Schindler et al., 2013)
Julida ‐ Ommatoiulus	Portugese millipede	Invasive detritivore species found in great abundance in the SCP, particularly in older restoration and reference sites. Feeds on litter, which is more available in those sites. The Portuguese millipede is known to be widespread throughout southern Australia (Baker et al., 2013)
Fabaceae	Legume family	ASVs in this family are strongly associated with younger restoration in JF and PB. Acacia shrubs tend to establish rapidly at restored sites in these locations (data from BHP, Data from South32)
Goodenia microptera	Narrow‐winged Goodenia	An insect pollinated species found predominantly in the vane traps of PB reference sites. Common Pilbara arid plant species, known to flower between February–October (Barrett & Barrett, 2014)
Anigozanthos	Kangaroo paw	Associated with younger restoration in SCP pitfall traps. These grow quickly (within a year) in SCP restoration. Anigozanthos are known to be predominantly pollinated by nectar‐feeding birds (Ayre et al., 2020) and could be encouraging of faunal recolonization of restored sites

Abbreviations: JF, Jarrah Forest; PB, Pilbara; SCP, Swan Coastal Plain.

Unfortunately, we cannot identify which invertebrates are interacting with which plant species. This would require isolating invertebrates and extracting DNA from each specimen separately (Bell et al., 2017; Pornon et al., 2016). Alternatively, eDNA from vegetative surfaces could be used to detect the associated invertebrates, for example, using flowers to identify possible pollinators (Thomsen & Sigsgaard, 2019). However, these studies require species‐specific sampling and therefore far more samples and greater costs. Overall, this study demonstrates that using bulk arthropod samples is a cost, time and resource efficient method that allows researchers to gain an informative snapshot of the invertebrate community and the plants they utilise.

Limitations

We were able to demonstrate how DNA metabarcoding can reveal restoration trajectories of invertebrate communities and provide useful information on their associated plant communities. However, there were some limitations. For example, as this study was conducted in the period of optimum plant growth and flowering, we cannot confirm whether the same patterns exist throughout the year. Seasonality affects invertebrate communities (Santorufo et al., 2014; Shimazaki & Miyashita, 2005), plant communities, and especially the interaction between the two (CaraDonna et al., 2017; Rico‐Gray et al., 1998). A previous study conducted during autumn (April) in the coastal plain sites using pitfall traps also detected directional changes in invertebrate communities (Fernandes et al., 2019), but no differences in plant communities from reference or restoration sites (Fernandes et al., 2019). This study offers preliminary testing of consistency in restoration patterns across space, but not within or between years and seasons.

While this and other studies (e.g., Beng et al., 2016) have demonstrated the utility of taxonomic independent analyses to investigate changing community profiles, a lack of taxonomic information limits the utility of these data. For example, more complete reference libraries allow greater resolution of taxonomic assignment, enabling species rather than family level identifications (Dormontt et al., 2018), and reducing the number of unassigned ASVs that may be removed from subsequent analyses (Schenekar et al., 2020; Stoeckle et al., 2020). Populating barcode reference libraries is a solution to this issue but it is especially challenging for invertebrates because of their high diversity and difficulty finding and funding the taxonomic expertise necessary for identifications (Austin et al., 2004). Arguably, demonstrating the utility of DNA metabarcoding for such monitoring projects, as in our study, provides enhanced incentive and need to further support barcoding efforts, systematics and taxonomy to improve representation of insect diversity in publicly available reference databases.

Greater taxonomic identification would also allow analyses to be conducted based on species identification rather than sequence based units such as ASVs or OTUs. Similarly, ASV or OTU analyses provide the means to measure sequence diversity. The fact that there may be multiple ASVs in a single species (Callahan et al., 2016), especially in a highly variable region such as COI, may result in an overestimation of diversity and influence the results of richness. For example, the same species may occur in both restoration and reference sites, but with sequence variants (i.e., haplotypes) amongst the gene region studied, in which case some ASVs of that species would be present in the reference sites but not the restoration sites. However, the striking similarity of trends detected in this study and previous morphology based projects (Davieson & Majer, 1983; Dunlop et al., 1985; Majer, Gunawardene, et al., 2013) gives us a high degree of confidence that the trends reported here are real and meaningful.

Finally, this study emphasises the need for restoration projects to be designed to test questions critical to restoration ecology. Unfortunately, ad hoc study systems as used in this study are a necessary approach (due to the slow maturity of these shrub dominated ecosystems) rather than deliberately designed experimental study systems (Prober et al., 2018). However, this study does provide confidence in the potential benefits and limitations of using DNA metabarcoding to monitor invertebrate recovery. This includes showing where this method may demonstrate recovery trajectories (Mediterranean woodlands and forests), as well as where it may fail to do so (hot, arid deserts). Further development will be required to control variation between different seasons, and to provide comparative equivalency between different ecosystems.

Conclusion

We have demonstrated the use of high throughput sequencing of invertebrate samples to establish restoration trajectories. Defining the likely trajectory of a restored site is important as it enables the definition of success criteria, and the required time scales for restoration monitoring. We show that trajectories towards reference ecosystems were more evident in ground dwelling invertebrates in older restored sites. Despite the lack of abundance data, metabarcoding can detect recovery of ecosystem function by showing whether invertebrates are interacting with the plant community. Understanding restoration trajectories using DNA metabarcoding will require additional validation research to determine the effects of seasonal variation, and consistency of patterns across multiple years and different ecosystems. Further, because ecosystems are dynamic, determining whether sites have been fully restored depends heavily on the selection of appropriate reference sites to capture natural variation in the reference ecosystem. The Bonn Challenge goal to restore 350 million km2 of degraded terrestrial ecosystems by 2030 (Suding et al., 2015) and ambitions of the UN Decade of Ecosystem Restoration means that effective tools such as metabarcoding are necessary to audit, manage and to inform interventions when trajectories are failing while protecting the considerable investments needed to meet these ambitious global restoration targets. Refining the emerging toolkit of rapid monitoring techniques such as DNA metabarcoding and evaluating where they are beneficial is critical to incorporation in restoration projects, to ultimately improve restoration outcomes.

CONFLICT OF INTERESTS

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

AUTHOR CONTRIBUTIONS

Mieke van der Heyde conducted the study and wrote the manuscript. Mieke van der Heyde, Paul Nevill, Michael Bunce, Nicole E. White, and Grant Wardell‐Johnson were involved in the experimental design. Samples were collected and processed by Mieke van der Heyde; molecular and bioinformatic work was performed by Mieke van der Heyde; all data was analysed and processed by Mieke van der Heyde; statistical analysis was done by Mieke van der Heyde; the manuscript was edited by all authors.

Supplementary Material

Click here for additional data file.