Density-dependent resource partitioning of temperate large herbivore populations under rewilding

In tropical grazer assemblies with abundant large predators, smaller herbivores have been shown to be limited by predation and food quality, while the larger species are regulated by food abundance. Much less is known on herbivore resource partitioning in temperate grazing ecosystems, where humans are typically the regulators. In the Oostvaardersplassen ecosystem in The Netherlands, a unique multispecies assemblage of cattle, horses, red deer and geese developed after initial introduction of a few individuals in 1983. During the first 35 years, this herbivore assemblage without predation or human regulation gradually changed into increasing dominance of the smaller herbivore species. Carrying capacity was reached around 2008, after which numbers started fluctuating depending on winter conditions. A population crash, especially of red deer, in winter 2018 led to heavy societal debate around animal welfare, after which active population regulation was introduced. This suggests strong niche overlap and competition between these very different-sized herbivores, possibly due to their homogenising effect on vegetation composition and structure at high densities. We used eDNA metabarcoding of dung to quantify the diet composition of cattle, horse, red deer and geese, annually in early winter from 2018-2021 and calculated their niche overlap. Overall, we found strong diet overlap between species. The diet of horse and cattle remained mostly unaltered and it was the one of red deer that changed the most across the years. Niche overlap decreased with increasing red deer population size, the most abundant species. When calculated as total energy expenditure, we found niche overlap was more linked to the shifts in red deer than to the total herbivore energy fluctuation. We suggest red deer changed their diet mainly in response to their own population size, reducing their niche overlap with increasing red deer population. In this case, resource competition translated into shorter vegetation height, reducing resource availability and forcing herbivores to consume different plant taxa. We conclude that in this temperate ecosystem, inter- and intraspecific resource competition are key factors structuring community composition and dynamics from small to large herbivores, with a competitive advantage of the smaller species, but with also opportunities for resource partitioning. 

Methods
Dung samples were collected across the grassland part of the OVP, with collections divided in three sub-areas (see Suppl. Figure 1A) during November 2018, 2019, 2020 and 2021 for the four main herbivore species, i.e., cattle, horse, red deer and geese (Barnacle geese and Greylag geese combined as the species were not identified from the dung shape). Per species and year, 15 scat samples were collected (5 per sub-area, Suppl. Figure 1A), leading to a total of 60 scat per year. Samples were spaced by at least 10 m to reduce the chance of re-sampling the same individual, and GPS coordinates were taken for each sample. Only freshly deposited dung samples were collected, and samples were then stored in dried silica beads at room temperature, in order to dry and preserve them, without need for freezing, until DNA extraction could be done at the University of Lausanne, Switzerland.

DNA extraction

We used between 0.5 and 1 g of dry dung as the starting point for the extraction. Extractions were performed using the NucleoSpin Soil Kit (Macherey-Nagel, Düren, Germany) following the manufacturer protocol. A subset of the extractions was tested for inhibitors with quantitative real-time PCR (qPCR) applying different dilutions (2x, 10x and 50x) in triplicates. qPCR reagents and conditions were the same as in DNA metabarcoding PCR reactions (see below), with the addition of 10,000-fold diluted SybrGreen (Thermo Fisher Scientific, USA). Following these analyses, all samples were diluted 5-fold before PCR amplification. All extractions were performed in a laboratory restricted to low DNA-content analyses.

DNA metabarcoding

DNA extracts were amplified using a generalist plant primer pair (Sper01, (Taberlet et al., 2018)), targeting all vascular plant taxon (Spermatophyta). Sper01 targets a 10-220 bp gene fragment of the P6 loop of trnL intron, chloroplast DNA. To assign the DNA sequences to each sample, primers were tagged with eight variable nucleotides added to their 5’-end with at least five differences between tags. The PCR reactions were performed in a final volume of 20 µL. The mixture contained 1 U AmpliTaq® Gold 360 mix (Thermo Fisher Scientific, USA), 0.04 µg of bovine serum albumin (Roche Diagnostics, Basel, Switzerland), 0.2 µM of tagged forward and reverse primers and 2 µL of 5-fold diluted template DNA. PCR cycling conditions were denaturation for 10 minutes at 95 °C, followed by 40 cycles of 30 s at 95 °C, 30 s at 52 °C and 1 min at 72 °C, with a final elongation step of 7 min at 72 °C. Amplifications were performed separately for each species and in replicates (4 per sample) in PCR plates each containing 60 DNA extracts, 12 blanks as well as 8 extraction, 8 negative and 8 positive PCR controls (DNA assembly of 10 plant species with increasing relative concentrations). The use of blanks allows estimating the proportion of tag switches (i.e., false combination of tags, generating chimeric sequences) during library preparation (Schnell et al., 2015). Amplification success and fragment sizes were confirmed on agarose gel. PCR products were subsequently pooled per PCR plate. Amplicons were purified using a MinElute PCR Purification Kit (Qiagen, Hilden, Germany). Purified pools were quantified using a Qubit® 2.0 Fluorometer (Life Technology Corporation, USA). Library preparation was done following the TagSteady Protocol (Carøe & Bohmann, 2020). After adapter ligation, libraries were validated on a fragment analyzer (Advanced Analytical Technologies, USA). Final libraries were quantified, normalised and pooled before 150 paired-end sequencing on an Illumina MiniSeq sequencing system with a Mid Output Kit (Illumina, San Diego, CA, USA).

Bioinformatic data analyses

The bioinformatic processing of the raw sequence output was performed using the OBITools package (Boyer et al., 2016). Initially, forward and reverse reads were assembled with a minimum quality score of 40. The joined sequences were assigned to samples based on unique tags combinations. Assigned sequences were then de-replicated, retaining only unique sequences. All sequences with less than 100 reads per library were discarded as well as those not fitting the range of metabarcode lengths. This was followed by two different clustering methods. First, pairwise dissimilarities between reads were computed and lesser abundant sequences with single nucleotide dissimilarity were clustered into the most abundant ones. Second, we used the Sumaclust algorithm (Mercier C, 2013) to further refine the resulting clusters based on a sequence similarity of 97 %. It uses the same clustering algorithm as UCLUST (Prasad, D.V., 2015) and it is mainly used to identify erroneous sequences produced during amplification and sequencing, derived from its main (centroid) sequence. Remaining sequences were assigned to taxa using a reference database. We built a database for Sper01 by running an in silico PCR based on all the plant sequences available in the EMBL database (European Molecular Biology Laboratory). We kept a single sequence per taxonomic id that was annotated at least to genus level.

Further data cleaning and filtering was done in R (version 4.0.2) using the metabaR package (Zinger et al., 2021). Sequences that were more abundant in extraction and PCR controls than in samples were considered as contamination and removed. Operational taxonomic units (OTUs) with similarity to the reference sequence lower than 97 % were also eliminated from the dataset. Removal of tag-leaked sequences was done independently for each library. This approach allowed us to discard single OTUs instead of whole PCR replicates. However, PCR replicates with too small reads count were also discarded. 

Remaining PCR replicates were merged by individual, keeping the mean relative read abundance (RRA), frequency of occurrence (FOO) and presence-absence.