Contrasting seasonal patterns in diet and dung‐associated invertebrates of feral cattle and horses in a rewilding area

Trophic rewilding is increasingly applied in restoration efforts, with the aim of reintroducing the ecological functions provided by large‐bodied mammals and thereby promote self‐regulating, biodiverse ecosystems. However, empirical evidence for the effects of megafauna introductions on the abundance and richness of other organisms such as plants and invertebrates, and the mechanisms involved still need strengthening. In this study, we use environmental DNA (eDNA) metabarcoding of dung from co‐existing feral cattle and horses to assess the seasonal variation in plant diet and dung‐associated arthropods and nematodes. We found consistently high diet richness of horses, with low seasonal variability, while the generally lower dietary diversity of cattle increased substantially during summer. Intriguingly, season‐specific diets differed, with a greater proportion of trees in the horses' diet during winter, where cattle relied more on shrubs. Graminoids were predominantly found in the diet of horses, but were generally underrepresented compared to previous studies, possibly due to the high prevalence of forbs in the study area. Dung‐associated arthropod richness was higher for cattle, largely due to a high richness of flies during summer. Several species of dung‐associated arthropods were found primarily in dung from one of the two herbivores, and our data confirmed known patterns of seasonal activity. Nematode richness was constantly higher for horses, and nematode communities were markedly different between the two species. Our results demonstrate complementary effects of cattle and horses through diet differences and dung‐associated invertebrate communities, enhancing our understanding of large herbivore effects on vegetation and associated biodiversity. These results are directly applicable for decision‐making in rewilding projects, suggesting biodiversity‐benefits by inclusion of functionally different herbivores.


| INTRODUC TI ON
Recent anthropogenically driven changes in the Earth's biosphere and climate system have imposed major threats on biodiversity and ecosystems on a global scale (IPBES, 2019;IPCC, 2021). There is an emerging consensus that conservation and restoration of natural processes are urgent for the safeguarding of natural habitats and biodiversity (Svenning, 2020). Accordingly, protecting and restoring self-regulating ecosystems with natural dynamics are increasingly in focus (Carver et al., 2021;Johns, 2019;Svenning et al., 2016).
Large-bodied animals ("megafauna") are important agents in these schemes. In terrestrial ecosystems, they shape plant communities through physical disturbance, herbivory, and seed dispersal, and are expected to promote heterogeneous, biodiverse landscapes (Bakker et al., 2016;Vera, 2000;Wigley et al., 2014). Within megafauna, functional diversity of the large herbivores is expected to be important for these effects (e.g., Bakker et al., 2016). As an important dimension to large-herbivore function, the presence, diversity, and relative abundance of grazers and browsers are expected to be important due to divergent ecological effects of their different feeding behaviours (Bakker et al., 2016;Hidding et al., 2013;Wigley et al., 2014).
Megafauna has been present in most of the world's terrestrial ecosystems since the mid Cenozoic and is expected to have shaped the evolution of many extant species (Galetti et al., 2018;Malhi et al., 2016;Smith et al., 2018;Titcomb et al., 2022). Since the Late Pleistocene, there has been a high extinction-selectivity for megafauna, with the result that most contemporary terrestrial ecosystems are deprived of megafauna and the functions they provide (Malhi et al., 2016;Smith et al., 2018). Thus, introductions of formerly present megafauna, or extant counterparts -known as "trophic rewilding" ) -have been suggested as a tool to achieve self-regulating systems and have already been applied in many ecosystems across the globe (e.g., Dvorský et al., 2022;Pedersen et al., 2020;Seddon et al., 2014).
In Europe, the species included in trophic rewilding efforts include species such as European bison (Bison bonasus), Eurasian elk (Alces alces), wild boar (Sus scrofa), cattle (Bos primigenius taurus), horse (Equus ferus caballus), water buffalo (Bubalus bubalus), and Asiatic wild ass (Equus hemionus). Feral cattle and horses are the most popular choices due to their availability, familiarity and because their wild ancestral forms (aurochs (Bos primigenius) and wild horses (E. ferus)) were widely distributed across Europe as late as the mid-Holocene (Sommer, 2020).
Multiple studies have shown positive effects of year-round grazing with these species on plant diversity and vegetation structure (e.g., Gilhaus et al., 2014;Köhler et al., 2016). In the Mols Laboratory in Denmark, where our study was conducted, exclosure experiments have shown higher plant species richness, indications of higher floristic uniqueness, and higher forbs-to-graminoids ratios in plots with year-round or winter grazing compared to mowing or summer grazing (Bonavent et al., 2022). These studies highlight the importance of natural grazing regimes for plant species diversity.
The complementarity of different herbivore species has been investigated by Cromsigt et al. (2018), who found that European bison and cattle included a significant amount of woody material in their diet, whereas horses almost exclusively grazed. This suggests that these large herbivores occupy different natural feeding niches, and consequently have differing influences on vegetation structure. However, other studies have shown high potential for competition and niche overlap between large herbivores (Menard et al., 2002;Scasta et al., 2016). These contrasting findings suggest that the interactions between large herbivores and their ecological impacts may be highly population-and context-specific and emphasize the limitations in our current understanding of how these species interact with each other and the ecosystems in which they occur. Thus, more empirical studies on the ecology of these species are needed, such that appropriate species can be selected for rewilding projects in different contexts.
A large part of Earth's biodiversity is composed by insects and other invertebrates. Trophic rewilding is expected to have positive effects on invertebrate diversity through promoting heterogeneity in vegetation and soils and the generation of biotic microhabitats such as dung and carcasses (e.g., Svenning et al., 2019). However, more studies testing these effects of different assemblages of large herbivores on various groups of invertebrates are needed (Nickell et al., 2018;van Klink et al., 2015). Importantly, there may be cases of real-world implication that cause less positive outcomes (as discussed in van Klink et al., 2015). It has been shown that continuous, extensive grazing positively affects dung-associated faunal communities (Buse et al., 2015;Jay-Robert et al., 2008). Numerous invertebrate species are directly or indirectly associated with the dung from large herbivores (Byk & Piętka, 2018;Skidmore, 1991), and land use change has been identified as a major driver of the decline of dung beetles (Carpaneto et al., 2007;Lobo, 2001), with veterinary treatment of animals also playing an important role (Tonelli et al., 2020;Verdú et al., 2018). Furthermore, it has been shown that downsizing of dung beetles during the late Quaternary is associated with megafauna losses, highlighting the link between these two species groups . Hence, the importance of management practices, including continuous presence of large herbivores suggests that trophic rewilding might be a beneficial management strategy for dung-fauna conservation. Nevertheless, our understanding of outcomes of trophic rewilding for dung-associated fauna would benefit from broader taxonomic coverage and further studies of different large-herbivore assemblages in a variety of natural contexts (Nickell et al., 2018;van Klink et al., 2015).
DNA metabarcoding of environmental samples (hereafter eDNA metabarcoding: Taberlet et al., 2012), has been established as an efficient and nondestructive method for retrieving community-level information from various substrates (Taberlet et al., 2018;Thomsen & Willerslev, 2015). In a trophic rewilding context, this method has recently been used to characterize the diets of bison, moose, and horses (Hagstrup et al., 2020;Hartvig et al., 2021;Kowalczyk et al., 2011Kowalczyk et al., , 2019, and thereby evaluate their influence on vegetation structure. In addition, Sigsgaard et al. (2021) showed that eDNA metabarcoding of cowpats can further be used to study dung-associated invertebrate assemblages. Thus, eDNA metabarcoding has great potential for investigating the effects of trophic rewilding on natural habitats and biodiversity.
Here, we investigate seasonal change in diet composition and dung-associated invertebrate communities for coexisting feral cattle and horses in a rewilding area with year-round grazing and without supplementary feeding using eDNA metabarcoding of dung samples.

| Study site
The study was performed in Mols Bjerge National Park, Denmark, in the area known as "The Mols Laboratory" (56 o 13′36"N, 10 o 34′33″ E; Figure 1), which is owned by the Natural History Museum of Aarhus.
From 2016, the area has been managed using a trophic rewilding approach , with feral Galloway cattle and Exmoor ponies ( Figure 1) roaming a 120-hectare enclosure. No supplementary feeding and no medical treatment of the animals (e.g., treatment against worms or antibiotic treatments) took place during the study period (Supporting Information S1). Approximately 50% of the area consists of unmanaged forest, including some parts dominated by common ash (Fraxinus excelsior) and alder (Alnus glutinosa), and generally high coverage of European beech (Fagus sylvatica) and common oak (Quercus robur). The rest of the area contains a mixture of pastures, moors, meadows, and fens (Naturhistorisk Museum, 2021b). The pastures host a variety of forbs as well as shrubs, such as raspberry (Rubus idaeus), blackberry (Rubus fruticosus), roses (Rosa spp.), and common broom (Cytisus scoparius). The moors have a high abundance of common heather (Calluna vulgaris) and wavy hairgrass (Avenella flexuosa). The meadows and fens are characterized by a high richness of plant species, especially forbs (Naturhistorisk Museum, 2021b). From extensive field observations and knowledge about the area obtained from highly regular monitoring throughout the last decade, we estimated extent of occurrence (EOO) and area of occupancy (AOO) for all registered plant genera. We subsequently assigned an abundance category ("highly abundant", "common", "rare" or "absent since 2018"), based on a combination of these estimates. The "highly abundant" category included genera with EOO estimates >75% of research area and AOO estimates >75% of suitable habitat, "common" included genera with EOO estimates from 25-75% and AOO estimates >25%, and "rare" included genera with estimated values below 25%. Some genera with a limited distribution, but with a vivid growth and estimated >1000 kg of living biomass available to the grazers (Lotus, Mentha, Juncus, Prunus, and Nardus) were upgraded from rare to common.
Changes in vegetation structures from dominance by graminoids to forbs and increased abundances of invertebrates have been observed in the area after the megafauna introductions.

| Sampling
Each month in 2020, 10 cowpats and 10 horse droppings were sampled in the study area (see Supporting Information S2 and Figure S9 for examples of sampled dung pats). Dung pats of approximately 3-4 days of age were selected based on visual appearance (wetness, crust formation) to detect the most colonizing species, while minimizing the degradation of DNA, and to increase the probability of targeting different herbivore individuals in each round of sampling.
For each of the two species, we collected five dung samples from open habitats and five from forest-, or thicket-covered habitats.
Furthermore, in every sampling month (except January and March) a field control was collected by holding a plastic knife out in the air for 30 s, then placing it inside a falcon tube for 30 s, holding it outside the tube for 30 s again, and finally sealing it in the tube. Face masks and nitrile gloves were worn during sampling and gloves were replaced between every sample. Samples were kept in a cooler bag with ice packs until returning from the field, after which they were placed in a freezer and stored at -20°C until DNA extraction.

| DNA extraction
DNA from the dung samples was extracted using the QIAamp Fast  Table S1 for amounts added). Afterwards, ~220 mg of dung from each sample was transferred to an Eppendorf tube using a metal spatula. Spatulas were cleaned between handling of different samples by placing them in a 3% chlorine bath for at least 10 min and rinsing with 70% EtOH.
The manufacturer's DNA extraction protocol was then followed, except for 2 h of incubation at room temperature after the addition of InhibitEX buffer, followed by 5 min of centrifugation at 14,000 rpm (18,407 g) . Final elution was done in two rounds, each time adding 60 μL ATE buffer, incubating for 5 min at room temperature and centrifuging for 1 min. The samples were extracted in 20 batches, and an extraction blank (CNE1-20, Table S1) was included in each extraction batch. The final eluates were stored at -20°C until PCR amplification.

| PCR amplification and sequencing
The 240 samples, 10 field controls, and 20 extraction blanks were divided across five "PCR sequencing units" (hereafter PSUs), each containing 54 samples (including field controls and extraction blanks) and four PCR blanks. DNA was amplified using two primer sets: the ITS2-S2F/ITS4 primer set (ITS) targeting plants (Fahner et al., 2016;Supporting Information S4) and the BF1/BR1 primer set (COI) targeting arthropods (Elbrecht & Leese, 2017). The ITS primer set was evaluated by in silico amplification (Ficetola et al., 2010; see Supporting Information S5 and Figure S10). Based on preliminary PCR tests, samples were diluted 1:10 for the ITS amplification, while for the COI amplification, the undiluted DNA eluates were used.
Primers were uniquely tagged before PCR amplification, with tags of six nucleotides designed with the OligoTag program (Coissac, 2012).
Two or three random bases were added before each tag to increase sequence complexity (De Barba et al., 2014). To minimize the problem of tag jumps, identical tags were used on the forward and reverse primer (Schnell et al., 2015). Annealing temperature for the ITS2-S2F/ITS4 primer set was set following Fahner et al. (2016), and for the BF1/BR1 primer set, the reaction volumes and thermal settings followed Sigsgaard et al. (2021) with only minor adjustments (Table S2). Four replicate PCR runs were performed for each PSU, resulting in four PCR replicates of each sample. Tags were used only once within each replicate but were reused between replicates. The PCR products from each replicate were pooled, resulting in 40 pools of PCR products (20 from each primer set). Each pool was thoroughly mixed and a subsample of 100 μL was then taken and purified using the MinElute PCR purification kit from Qiagen. The manufacturer's protocol was followed, except for the final elution step, which was done in two rounds, each time adding 20 μL elution buffer (EB), incubating at 37°C for 10 min, and centrifuging for 1 min. Library building and sequencing on an Illumina NovaSeq 6000 was performed by Novogene using 250 PE sequencing for ITS and 150 PE sequencing for COI. For each library, an output of 10 Gb was requested.

| Data processing
Sequencing data was processed through the MetaBarFlow pipeline (Sigsgaard et al., 2022, Supporting Information S6), which contains demultiplexing based on Cutadapt (Martin, 2011), quality trimming using the sickle tool (Joshi & Fass, 2011) and error filtering with DADA2 (Callahan et al., 2016). After initial filtering and merging of paired reads, the amplicon sequence variants (ASVs) of the ITS data set were searched against a locally downloaded version of the GenBank nt database (Sayers et al., 2020) released on 19 November 2020, while the ASVs in the COI data set were searched against a custom COI database, built using the MARES database pipeline (Arranz et al., 2020), containing a combination of eukaryotic COI sequences from the GenBank nt database and the BOLD database (Ratnasingham & Hebert, 2007; www.bolds ystems.org). Details about search terms can be found in Klepke et al. (2022). Searches were done using the BLASTn algorithm (Altschul et al., 1990), requesting a maximum of 500 hits, minimum 90% query coverage per high-scoring segment pair, and minimum 80% sequence similarity.
ASVs were then taxonomically classified using the R package taxizedb (Chamberlain & Arendsee, 2021). Each ASV was assigned to the most recent common ancestor of the taxon (or taxa) yielding the hit with highest sequence similarity and any taxa showing an overlap in sequence similarities with those of the best-matching taxon . Only blast hits with ≥90% identity and 100% query coverage (≥97% query coverage for the COI data set, to increase the amount of reference sequences considered) were included in the taxonomic assignment. Species level identification was only made if sequence similarity was at least 98%.
Analysis of the ITS data set was made at the genus level (included 94.23% of all ASVs), as species-resolution identification was not possible for many ASVs (77.64%), and plant functional groups could be inferred from genus-level identifications. For the COI data, analyses were made at the species level (36.9% of all ASVs identified to species level), as species-level resolution was important for evaluation of biodiversity and ecological function. For the COI data set, all ASVs were checked for errors in the taxonomic assignment due to erroneous sequences in the database, and all ASVs identified to a higher taxonomic level than the level used for further analysis were aggregated into putative species (COI) or genera (ITS) (e.g., Aphodius sp.1), following Jensen et al. (2021). For the ITS data set, only the 500 ASVs with the highest read counts, which represented 86% of total reads, and some additional ASVs with clearly erroneous identifications, were checked manually and aggregated, and for the rest of the ASVs, the ID assigned by the automated pipeline was used.
After these manual corrections, sequences present in controls (field blanks, extraction blanks, and PCR blanks) were removed from all samples if the maximum read count in a single control sample was higher than the maximum read count in any dung sample (Table S3).
For the ITS data set, sequences with less than 100 reads in any sample were removed to decrease the probability of false positive results (DNA from other sources than the dung itself) as suggested Accumulation curves were plotted for all samples using the function specaccum from vegan (version 2.5-7).

| Statistical analysis
For the ITS data set, statistical analysis was only made on sequences belonging to the group Streptophyta (91% of all reads). The year was divided into three seasons; "early" (January-April), "mid" (May-August), and "late" (September-December) to represent the differences in available plant biomass over the year. As the data were not normally distributed, the nonparametric Kruskal-Wallis test was used to test for the influence of herbivore species (cattle/horse), and season (early/mid/late) on genera richness. Subsequently, Wilcoxon rank-sum tests were used to identify pairwise differences, if Kruskal-Wallis tests indicated that a certain predictor had an influence.
Bray-Curtis distances were calculated using the function vegdist on the rarefied and square root transformed read counts.
Differences in plant diet composition were then evaluated through an nMDS ordination based on the Bray-Curtis distances, using the function metaMDS, specifying three nMDS axes (lowest number of axes which yielded stress<0.2). Differences between samples were evaluated by permutational analysis of variance (PERMANOVA, number of permutations = 999) using the function adonis2 on the Bray-Curtis distances. Again, herbivore species and season were used as predictors. The assumption of multivariate homogeneity of dispersion in PERMANOVA was evaluated by permutation tests with the functions betadisper and permutest, and none of the assumptions of the PERMANOVA test were violated.
Differences in the relative proportions of these functional groups of plants between herbivore species (cattle/horse) and seasons (early/ mid/late) were evaluated by the nonparametric Kruskal-Wallis test, with subsequent pairwise evaluations by Wilcoxon rank-sum tests.
For the COI data set, all dung-associated arthropods (Skidmore, 1991) and nematodes were included in the analysis, as other species were expected to represent contamination, or random visits. For the nematodes, only ASVs with at least 95% similarity to a database sequence were included, due to the large database deficiency for this group (Kvist, 2013). First, a generalized linear latent variable model (gllvm; Niku et al., 2019), was fitted to the data, with species richness within each order represented in the COI data set (19 orders of dung-associated arthropods and nematodes) as the input matrix. Orders which were present in fewer than 20 samples were removed (10 orders), to minimize their disproportionate influence on the model fit due to underrepresentation. Herbivore species (cattle/horse), season (early/mid/late), and habitat (open/forest) were used as predictors in the model, and the best model (evaluated by AIC) was a model using a poisson distribution, with all three predictors, without their interactions, and including a random roweffect. The inclusion of a random row-effect was done to handle the risk of nonindependent sampling due to our sampling design.
Subsequently, we focused the analysis on the most relevant groups which represented >90% of all Metazoan reads ( Figure 1).
Thus, the richness of (i) all dung-associated arthropods (Arthropoda), (ii) dung-associated beetles (Coleoptera), (iii) dung-associated flies and allies (Diptera), and (iv) nematodes (Nematoda) were evaluated separately, to investigate differences between the two herbivore species in different seasons. As with the ITS data set, nonparametric tests were used, as data was not normally distributed. The richness of nematodes in the dung samples was described by a generalized linear model (GLM, poisson distribution) with month number (1-12, representing January-December), herbivore species and habitat as predictors, as preliminary inspection of the data showed a continuous increase in richness over the year.
Bray-Curtis distances were calculated for each group, and distance between samples were evaluated by nMDS ordination (two nMDS axes) and PERMANOVA (number of permutations = 999), as with the plant data. Here, multivariate dispersion was not homogenous, so results should be interpreted with caution.
Finally, species that were present in five or more samples, and for which >90% of their total read count was present in either cattle or horse samples or in either open or forest samples, respectively, were defined as being predominant in dung from that herbivore species or habitat. All p-values from the nonparametric tests were corrected using the Holm-Bonferroni method, to minimize the risk of type 1 errors because of multiple testing (Holm, 1979). All functions referred to in this section were from the R-package vegan version 2.5-7 (Oksanen et al., 2020), and all analyses were performed in R version 3.6.2 (R Core Team, 2019).

| Sequencing output
For the ITS data set, a total of 823,972,154 raw reads were ob-  Tables 2 and S5). After inspection of reads matching the two herbivore species (cattle and horse), the assumed herbivore species for seven samples was corrected to match the presence of either horse or cattle reads (Table S1). The distribution of the obtained reads between different taxonomic groups are shown for both data sets in Figure 1. Rarefaction and accumulation analyses did not show any serious concerns (see Supporting Information S8). In silico analyses showed somewhat higher mismatches with the primers for sequences from grasses (Poales) and conifers (Pinales) compared to dicots (Magnoliopsida) ( Figure S10).

| Plant diet
A total of 133 plant genera were found in the dung samples (Table S4)  Note: NA, not present in the diet of the respective species; total reads, total number of nonrarefied reads obtained from each order across all samples.

TA B L E 2
List of all phyla, classes and orders of dung-associated arthropods and nematodes identified in this study. Number of families, genera and species are shown, as well as the number of species within each order that were predominantly detected (more than 90% of total reads) in samples from either cattle or horses, or from either open or forest habitat.   and Bellis, the coarse grasses Phalaris, Nardus, and Arrhenatherum, and the tree genus Fraxinus. We found 55% of all registered plant genera in the area in the dung samples, 94% of the highly abundant (only Achillea was not found), 49% of the rare, and none of those not registered since 2018. Relative proportion of reads in each month are shown for the 10 most abundant genera, and a few selected genera of interest in Figure S6.

| Diet richness and composition
The mean number of consumed plant genera was significantly higher for horses compared to cattle in early and late season (Wilcoxon test: p < .001), but not in mid-season (Wilcoxon test: p > .05, Table S6, Figure 2). PERMANOVA tests showed that 13.9% of the variation in plant composition between dung samples could be explained by the season (early/mid/late, p < .001), 8.8% by the herbivore species (cattle/horse, p < .001), and 5.7% by an interaction between herbivore species and season (p < .001) (Figures 2 and S2). The proportion of trees was higher in the diet of horses compared to cattle in early and late season (Wilcoxon tests: p < .01), but not in mid-season (Wilcoxon test: p > .05, Table S6, Figure 3). For shrubs, higher proportions were found in the diet of cattle in early season, compared to horses (Wilcoxon test: p < .01), but not in mid and late season (Wilcoxon tests: p > .05, Table S6, Figure 3). Higher proportions of forbs were found in mid-season, compared to early and late season, for both species (Wilcoxon tests: p < .001, Table S6, Figure 3), and for graminoids, higher proportions were found in the diet of horses, compared to cattle in all seasons (Wilcoxon tests: p < .01, Table S6, Figure 3). p-values from all performed tests are shown in Table S6.  (Skidmore, 1991). Thirty-three of the detected dungassociated arthropod species were predominantly (defined as ≥90% of all reads from the species in question) found in dung from cattle, while 18 species were predominantly found in dung from horses. A total of 20 species were predominantly found in dung from open habitats, while 11 species were predominantly found in dung from forest habitats ( Table 2). Read occurrence of selected species were found to match known phenology, habitat, and host preferences ( Figures S7 and S8). Relative read proportions of each sample and monthly means of the 15 most abundant species are shown in Figure S3.

| Richness and community composition of dung-associated arthropods
The multivariate model showed higher richness in cattle sam- Entomobryomorpha, while no difference between habitat types was found for any other orders (Figure 4). Higher richness was found between all seasons for the orders Mesostigmata and Diptera, whereas for Coleoptera, richness was only higher in mid-season compared to early and late season. For Poduromorpha, richness was higher in late season, compared to early season ( Figure 4). Only for the springtails in the order Symphypleona, higher richness was found in early season, compared to mid-season ( Figure 4). For all dung-associated arthropods, significantly higher species richness was found in cattle samples compared to horse samples in mid-season (Wilcoxon test: p < .001), but no difference was found between the two herbivores in early and late season (Wilcoxon tests: p > .05, Table S6, Figure 5a).
For both cattle and horse samples, higher arthropod richness was found in mid-season, compared to early and late season (Wilcoxon tests: p < .05), and for cattle, there was significantly higher richness in late season, compared to early season as well (Wilcoxon test: p < .001, Table S6, Figure 5a). For Coleoptera, there was no effect of herbivore species and habitat on species richness (Kruskal-Wallis tests: p > .05), but richness was significantly higher in mid-season, compared to early and late season for both cattle and horse samples (Wilcoxon tests: p < .05, Table S6, Figure 5b). This was consistent with the results from the multivariate analysis ( Figure 4) Figure 5c).

| Richness and community composition of nematodes
The multivariate analysis showed slightly higher richness in cattle samples compared to horse samples for Rhabditida, and higher richness in horse samples for Strongylida (Figure 4). No difference in richness was found between habitat types for any nematode orders, and for both Rhabditida and Strongylida, higher richness was found in mid-and late season, compared to early season ( Figure 4)

F I G U R E 3
Mean relative proportions of reads from different functional groups of plants (a,b) and different plant genera (c,d) for each month in dung samples from cattle (a,c) and horse (b,d). The 30 most abundant plant genera were assigned to four functional groups, while less abundant genera, and taxa that could not be sufficiently identified to confidently assign a functional group were included in the category "Unknown". In (c,d), the 10 most abundant genera are shown, and the rest are aggregated as "Other". Read proportions in individual samples are shown in Figure S1.

| DISCUSS ION
The extinction and extirpation of wild megafauna since the Late Pleistocene has resulted in terrestrial ecosystems worldwide being deprived of functions such as natural herbivory and physical disturbances (Dirzo et al., 2014;Sandom et al., 2014;Smith et al., 2018).
Introduction of large herbivores is therefore increasingly used in conservation contexts with the aim of restoring trophic interactions and disturbances, and promoting self-regulating systems (Pedersen et al., 2020;Svenning et al., 2016). However, there is an urgent need for more data on effects of such trophic rewilding in different areas, implementations, and contexts to enhance our understanding of the mechanisms of herbivore impacts, especially the impacts on invertebrates.
Using eDNA metabarcoding, we demonstrate contrasting patterns between cattle and horses in both diet and dung-associated biodiversity, suggesting niche differentiation and complementary vegetation and biodiversity effects of these two herbivores.
Our results further indicate that eDNA metabarcoding could be an efficient means for community characterisations in long-term monitoring of rewilding areas, and thereby for addressing the questions and uncertainties associated with future rewilding initiatives.

| Diet composition and richness
Pronounced differences in diversity and composition of diet were detected between cattle and horses. For cattle, diet richness was significantly higher in mid-season, compared to early and late season, whereas for horses, there was no difference between seasons. This suggests that the two species have different strategies for maintaining energy input during the low-productive seasons.
Our data suggest that horses feed opportunistically, as inferred by the high diet richness throughout the year (Figure 2), whereas the cattle were more selective (Figures 2 and 3). As an example, cattle increased the proportion of high-nutritious Rubus species (Kamler & Homolka, 2011;Poli et al., 1996) in their diet during winter, which might only be possible due to their digestive system as ruminants (Parish et al., 2022). We found relatively low amounts of F I G U R E 5 Mean species richness of all dung-associated arthropods (Arthropoda) (a), beetles (Coleoptera) (b), flies and allies (Diptera) (c), and nematodes (Nematoda) (d) in each month, for cattle samples (green) and horse samples (blue). Error bars show standard error of the mean. For each season, it is indicated whether there was a statistically significant difference (p < .001, (***)) in richness between the two species or no significant difference (p > .05, ns.). Graphical insets from https://openc lipart.org/ and https://commo ns.wikim edia.org/ with minor edits.
graminoids in the diet of both herbivores (Figure 3), unlike several studies showing a dominance of grasses in cattle and horse diets (Cromsigt et al., 2018;Scasta et al., 2016). Comparisons of eDNA and morphology-based analyses of feral horse dung have revealed potential preference of forbs over graminoids in eDNA analyses (King & Schoenecker, 2019). These disparities might be a result of the primers used having a reduced affinity to graminoid sequences, with artificially increased forbs to graminoids ratios as the result (see Supporting Information S5 and Figure S10). However, several of our samples showed high frequencies of graminoid sequences (e.g., AUG20, SEP3, SEP17), suggesting that grasses were amplified when present. In our study area, graminoids have been greatly reduced in distribution and abundance since the introduction of horses (Naturhistorisk Museum, 2020, 2021a). We speculate that heavy grazing before the onset of vegetation growth in spring, has caused a die-off of graminoids, leading to a dominance of forbs in the vegetation at the study site. This fits well with the suggested dominance of forbs on the Pleistocene steppe, where natural populations of large herbivores were present (Bråthen et al., 2021). Similar to the findings of Cromsigt et al. (2018), we detected a larger proportion of graminoids in the diet of the horses compared to the cattle, suggesting that graminoids are a relatively more important food resource for horses. It has recently been shown that horses can strongly reduce the dominance of coarse grasses such as Calamagrostis epigejos (Chodkiewicz, 2020;Henning et al., 2017), consistent with our findings in this study.
The plant species broom (Cytisus scoparius) and beach rose (Rosa rugosa), which have been a major management issue in the study area (Naturhistorisk Museum, 2021a). However, after the initiation of the rewilding programme, no control of these species has been performed. We found Cytisus sp. sequences in the horse samples, especially in the early season, which suggests that they to some ex- Horses have been shown to affect tree composition and colonization through browsing in rewilding areas (Garrido et al., 2021). Our data showed that the horses had a higher proportion of trees in their diet in early and late season compared to cattle (Figure 3, Table S6), which supports these findings. The most abundant genus of tree species in the diet of the horses was oak (Quercus), which is highly abundant in the study area. Previous investigations at the study site indicated that Scots pine (Pinus sylvestris) was one of the three most abundant plants in the horses' diet during winter (Hagstrup et al., 2020), but in this study we found very low read proportions from gymnosperms in horse dung samples. However, in the period immediately after the introduction of horses in 2016, they were ob-  (Somodi et al., 2008).

| Dung-associated communities
We found differences in dung-associated arthropod communities between the two herbivores ( Figure 6). Additionally, several species were predominantly found in samples from either cattle or horses, suggesting that the presence of both herbivores promotes biodiversity of dung-associated arthropods in rewilding systems. The dung from cattle was preferred by many flies, as well as some of the scarab beetle species (Scarabaidae: Aphodius, Onthophagus) whereas horse dung was preferred by two species of featherwing beetles (Ptiliidae: Baeocrara variolosa, Ptenidium nitidum), several Aphodius species, and some of the rove beetles (Staphylinidae) ( Table 2). Dung-associated beetles are able to detect dung from certain herbivores by olfactory cues (Dormont et al., 2004), and many Aphodius species are known to prefer certain dung types (Skidmore, 1991). This suggests that our detections accurately could indicate host preferences. Furthermore, we found seasonal differences in the dung-associated arthropod communities for both species (Figure 6), but the differences were more pronounced for cattle dung. This corresponds to the patterns found in the diet of the two herbivores (Figure 2), and we speculate whether a link could exist between herbivore diet composition and dung-associated arthropod communities. The digestive physiology of cattle allows for selection of more nutritious plants when available (Parish et al., 2022), which could potentially result in higher seasonal differences of available nutrients in cattle dung compared to horse dung, and dung-associated communities might respond to this. Additionally, we found a host shift for predatory beetles, which became highly abundant in cattle dung in March-April, after which they started to appear in dung from horses ( Figure S5), again F I G U R E 6 Non-metric multidimensional scaling (nMDS) plots and Venn diagrams for all dung-associated arthropods (Arthropoda) (a), beetles (Coleoptera) (b), flies and allies (Diptera) (c), and nematodes (Nematoda) (d) detected in the dung samples. First column: All samples, grouped by herbivore species. Second column: All samples grouped by habitat. Third column: Samples from cattle, grouped by season. Fourth column: Samples from horses, grouped by season.
suggesting complementary effects of the two herbivores. However, for investigating the direct effect of rewilding itself, long-term studies including samples taken pre-and post-implementation of rewilding at several independent sites should be performed.
Nematode richness and read proportions were consistently higher for horses compared to cattle (Figures 5d and S11). This might suggest that the parasitic load is higher for horses than cattle, which makes sense as horses are known to be parasitized by a large variety of strongyles (Family: Strongylidae) (Corning, 2009;Titcomb et al., 2022). This is also supported by our findings of higher richness of the order Strongylida in horse samples compared to cattle samples ( Figure 4). However, another contributing factor could be a better database coverage for horse parasites compared to cattle parasites. In this study, >50% of the obtained COI reads only had database matches with <80% similarity and were thus not identified. It has been shown that using 18S rDNA as a barcode resulted in far more identifications of free-living nematodes compared to COI (Ahmed et al., 2019), and more identifications could probably have been made by using this barcode. However, as the main goal of this study was to amplify dung-associated arthropods, we selected a marker that was better suited for arthropods due to higher database coverage (see Thomsen & Sigsgaard, 2019). To increase the possibility of identifying nematodes when using COI markers, we urge researchers working with parasites to publish COI sequences to reduce database deficiencies.

| Methodological considerations for eDNA metabarcoding of dung samples
Several considerations should be made when interpreting data obtained from eDNA metabarcoding (Alberdi et al., 2019;Deagle et al., 2019). First, the potential origins of the DNA sequences should be considered. It has been shown that eDNA from animals and plant pollen is transported by air (Clare et al., 2022;Johnson et al., 2019;Klepke et al., 2022;Lynggaard et al., 2022), and eDNA could also be transported to dung piles by visiting species. To address some of these issues, we included a field control to detect DNA input from the air, which contained no or very little DNA. However, these controls might only catch a fraction of the total diversity of airborne DNA at a site, which is expected to vary over very short spatial and temporal scales, due to the influence of abiotic factors on the transport of airborne particles (Burrows et al., 2009). Regardless, we expect that the relative amount of such background DNA is very low compared to DNA from organisms that have been consumed or have been in physical contact with the sample substrate recently. This is supported in our study by the detection of 94% of the highly abundant plant genera and 48% of the rare plant genera. Also, all detected plant genera have been recently observed in the area, except Callianthemum, which may represent DNA input from gardens nearby.
Additionally, it should be considered that no information about which plant parts were consumed is obtained from eDNA sequence data, which is an important limitation of this approach. Because some plant cell types are more easily lysed compared to others, the relationship between obtained DNA reads and ingested material is not the same across plant tissues . For example, it has been shown that leaves yield higher amounts of DNA compared to bark (Novaes et al., 2009 Stableton et al., 2022). This could be especially important when ruminants and nonruminants are compared as in this study, and to quantify such potential bias, controlled feeding trials should be performed. These challenges highlight that data obtained by molecular methods such as eDNA metabarcoding should be interpreted with care, and preferentially combined with data from traditional morphology-and observation-based methods (King & Schoenecker, 2019;Yoccoz, 2012).
Also, worth mentioning is the importance of high-quality reference databases, which was evident in this study. We found several problematic sequences in the GenBank nt database, which could potentially have reduced the taxonomic resolution to levels where biologically informative conclusions about diet and vegetation disturbance could not have been made (see Supporting Information S6).
We recommend that thorough quality checks of automated taxonomic assignments should always be performed to assess these issues. However, as sample sizes and sequencing output of eDNA studies increase, the number of ASVs is expected to increase accordingly, resulting in manual quality control becoming restrictively timeconsuming. This underlines the need for well-curated, high-quality reference databases as emphasized by multiple researchers in the field (Grant et al., 2021;Thomsen & Willerslev, 2015;Yoccoz, 2012).
Some initiatives have already been put in place, such as the sequence flagging system in BOLD (Milton et al., 2013, p. 25), or the BAGS program to assess the quality of reference libraries before taxonomic assignment (Fontes et al., 2020). We advocate that more resources should be invested in assessing the quality of reference databases, and that the tools already available should be used to ensure reliable taxonomic identification by automated approaches. An alternative approach is to construct geographically restricted and highly validated reference databases. This might be possible for diet assessments of herbivores in constrained, well-studied areas in Northern Europe, but not in remote areas, or for species-rich groups such as insects or nematodes.
Finally, bias introduced by the primers used could seriously impact results from eDNA metabarcoding studies. As an example, primers targeting COI have been shown to miss a relatively high number of target species of fish due to primer bias (Collins et al., 2019). Despite less problems with primer bias for other barcodes (e.g., 16 S rRNA), the low database coverage of these regions for the target groups in this study would have made species identifications difficult for many terrestrial arthropods (see Thomsen & Sigsgaard, 2019). Hence, we chose to use highly degenerated COI primers, which have shown consistent amplification and high detection rates of invertebrates, although these were designed for macroinvertebrate bulk samples and not eDNA (Elbrecht & Leese, 2017;Krehenwinkel et al., 2017).

| Conclusions and conservation perspectives
Trophic rewilding as a way of reintroducing key ecosystem functions is receiving rapidly growing interest and implementation as a For dung-associated arthropods and nematodes, important differences between the cattle and horses were also identified, illustrating the potential for increased biodiversity effects by using these two herbivore species. We thus suggest that grazing regimes with functionally different large herbivores should be prioritized in rewilding projects, due to their potential to increase vegetation heterogeneity and, in consequence, biodiversity. The highly accurate phenological patterns found among dung-associated arthropods ( Figures S7 and   S8) validate the eDNA approach, and warrant long-term studies to monitor effects before, during, and after implementation of rewilding initiatives. Future eDNA studies could also be used to reveal whether the composition of dung-associated arthropods is determined mainly by the type of herbivore (e.g., ruminant, or nonruminant), or by the composition of plants in the herbivore's diet, or to investigate patterns of invertebrate succession in herbivore dung. In conclusion, eDNA metabarcoding shows tremendous potential for biomonitoring of rewilding areas, and we encourage a wider implementation over larger spatial and temporal scales in comparison with other methods. Emil Ellegaard Thomassen led the writing of the manuscript with assistance and inputs from all the authors. All authors revised and approved the final manuscript.

ACK N O WLE D G E M ENTS
The authors would like to thank Trine Bech Søgaard and Annie Brandstrup for help with laboratory work, and employees at the Mols Laboratory for practical assistance. We also thank GenomeDK

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Raw sequencing data, files for demultiplexing, and the manually corrected taxonomy file have been made available on Dryad and can be accessed via: https://doi.org/10.5061/dryad.9p8cz 8wmb .