DOI: 10.5061/dryad.prr4xgxth
title: Data from: Density-dependent resource partitioning of temperate large herbivore populations under rewilding

description: 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. 
description: 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. <em>DNA
 extraction</em> 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.
 <em>DNA metabarcoding</em> 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).
 <em>Bioinformatic data analyses</em>
 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
 <em>in silico</em> 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.
description: # Density-dependent resource partitioning of temperate large herbivore
 populations under rewilding This README file was generated on 2024-02-07
 by Eduard Mas-Carrió GENERAL INFORMATION 1\. Title of Dataset: Data from:
 Density-dependent resource partitioning of temperate large herbivore
 populations under trophic rewilding 2\. Author Information Eduard
 Mas-Carrió1*, Perry Cornelissen2,5, Han Olff3,† and Luca Fumagalli1,4,†
 1\)Laboratory for Conservation Biology, Department of Ecology and
 Evolution, Biophore, University of Lausanne, 1015 Lausanne, Switzerland.
 2\)State Forestry Service, Amersfoort, The Netherlands. In 3\)Conservation
 Ecology Group, Groningen Institute for Evolutionary Life Sciences,
 University of Groningen, Groningen, Netherlands. 4\)Swiss Human Institute
 of Forensic Taphonomy, University Centre of Legal Medicine
 Lausanne-Geneva, Lausanne University Hospital and University of Lausanne,
 Ch. de la Vulliette 4, 1000 Lausanne 25, Switzerland. 5\)Institute for
 Biodiversity and Ecosystem Dynamics, University of Amsterdam, Netherlands.
   3\. Date of data collection: Nov 2018 – Nov 2019 – Nov 2020 – Nov 2021 
 4\. Geographic location of data collection: The Oostvaardersplassen, The
 Netherlands 5\. Funding sources that supported the collection of the data:
 FBM funding Vaud, Switzerland and SNF funding, Switzerland. 6\.
 Recommended citation for this dataset: Mas-Carrio et al. (2024), Data
 from: Density-dependent resource partitioning of temperate large herbivore
 populations under trophic rewilding, Dryad, Dataset DATA & FILE
 OVERVIEW 1\. Description of dataset 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, to put in perspective
 their dietary changes to the herbivore population in the
 Oostvaardersplassen. Sample collection 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. 2\.
 File List: Year 1 File name: Sper01_S1_L001_R1_001.ovp1.fastq.gz
 Description: Forward reads for 2018 sequences File name:
 Sper01_S1_L001_R2_001.ovp1.fastq.gz Description: Reverse reads for 2018
 sequences Year 2 File name: Sper_S1_L001_R1_001.ovp2.fastq.gz Description:
 Forward reads for 2019 sequences File name:
 Sper_S1_L001_R2_001.ovp2.fastq.gz Description: Reverse reads for 2019
 sequences Year 3 File name: Sper_S1_L001_R1_001.fastq.gz Description:
 Forward reads for 2020 sequences File name: Sper_S1_L001_R2_001.fastq.gz
 Description: Reverse reads for 2020 sequences Year 4          File name:
 ovp4-sper_S1_L001_R1_001.fastq.gz Description: Forward reads for 2021
 sequences File name: ovp4-sper_S1_L001_R2_001.fastq.gz Description:
 Reverse reads for 2021 sequences File name:
 ovp_plant.abundance_order_revised Description: metadata file with
 information associating plant species to their relative abundance in the
 Oostvaardersplassen. Includes plant species sequenced but not found in the
 grassland. File name: OVP_LH_trends Description: Metadata file with the
 herbivore countings for each species from 1983 until 2021 (correspond to
 counting on first of May every year, data used only for Figure 1). File
 name: herbivore_counts Description: Metadata file with the herbivore
 countings for each species from 2016 until 2021 (correspond to counting on
 first of October every year, data used for the modelling part). File name:
 OVP_env.full Description: Metadata file related to the scat samples
 collected. Informs on the exact GPS location, the species sampled, the
 quality of the scats, and further metadata. Comprises all 4 years
 together. Variables and units: \- Year: Year of sampling (1 =
 2018,2=2019,3=2020,4=2021) \- lat = latitude (GPS coordinates) \- lon =
 longitude (GPS coordinates) \- ele = elevation (in meters) \- time = time
 of sampling \- Sample = sample code \- species = herbivore species \-
 block = 1/2/3, describes the sub-regions within the OVP \- freshness =
 visual attribution of the freshness of the scats collected. Not used in
 the manuscript. \- day = day of sampling since the start of the sampling
 period for each year. \- extraction = extraction round that the sample was
 extracted within each year. \- sample weight = weight of the sample (in
 grams) File name: ngs_sper01_ovp1.txt Description: file to associate
 herbivore sequences to each individual sample using the forward and
 reverse unique tag combinations for year 2018 (OVP1) File name:
 ngs_sper01_ovp2.txt Description: file to associate herbivore sequences to
 each individual sample using the forward and reverse unique tag
 combinations for year 2019 (OVP2) File name: ngs_sper01_ovp3.txt
 Description: file to associate herbivore sequences to each individual
 sample using the forward and reverse unique tag combinations for year 2020
 (OVP3) File name: ngs_sper01_ovp4.txt Description: file to associate
 herbivore sequences to each individual sample using the forward and
 reverse unique tag combinations for year 2021 (OVP4) METHODOLOGICAL
 INFORMATION *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)](https://www.zotero.org/google-docs/?apsrQ4)), 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)](https://www.zotero.org/google-docs/?4p3jRS). 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)](https://www.zotero.org/google-docs/?HeUNw2). 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)](https://www.zotero.org/google-docs/?Hms03j). 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)](https://www.zotero.org/google-docs/?XyXhxy) 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)](https://www.zotero.org/google-docs/?KjGAdO) 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)](https://www.zotero.org/google-docs/?VW3axR). 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. DATA-SPECIFIC
 INFORMATION FOR: ovp_plant.abundance_order_revised 1\. Number of
 variables: 5 2\. Number of cases/rows: 164 3\. Variable List: Species:
 genus or species name of the plant taxon. Abundance: Relative abundance
 (in percentage) of that plant taxon in the OVP Genus: genus name of the
 plant taxon Order: order name of the plant taxon Ovp_presence: states if
 the plant taxon is found or not in the Oostvaardersplassen (OVP). 4\.
 Missing data codes:           None 5\. Abbreviations used:           NA;
 not applicable DATA-SPECIFIC INFORMATION FOR: OVP_LH_trends 1\. Number of
 variables: 6 2\. Number of cases/rows: 40 3\. Variable List: Year: Year
 Cattle: number of cattle individuals in the reserve Horse: number of horse
 individuals in the reserve Deer: number of red deer individuals in the
 reserve Barnacle geese: number of barnacle geese individuals in the
 reserve Greylag geese: number of greylag geese individuals in the reserve
 4\. Missing data codes:           None DATA-SPECIFIC INFORMATION FOR:
 herbivore_counts 1\. Number of variables: 8 2\. Number of cases/rows: 6
 3\. Variable List: Date: Month and year for the counting. Deer: number of
 red deer individuals in the reserve Cattle: number of cattle individuals
 in the reserve Horse: number of horse individuals in the reserve
 Barn.geese: number of barnacle geese individuals in the reserve
 Grey.geese: number of greylag geese individuals in the reserve Comments:
 Extra comments on the dynamics and context for the herbivore countings.
 4\. Missing data codes:           None 5\. Abbreviations used:          
 NA; not applicable DATA-SPECIFIC INFORMATION FOR: OVP_env.full 1\. Number
 of variables: 15 2\. Number of cases/rows: 284 3\. Variable List: Year:
 Year of sampling. Written as 1,2,3,4 (1 = 2018, 2 = 2019, 3= 2020, 4=2021)
 Lat=Latitude of the location of the scat sample Lon=Longitude of the
 location of the scat sample Ele= Elevation of the location of the scat
 smaple Time= Time of sampling of the scat sampling Sample= Sample code
 Species= herbivore species the scat belongs to Block= 1/2/3. Indicates
 which of the three subregions the sample was collected in, within the OVP.
 Freshness = visual identification of the freshness of the sample. +/- are
 used to subdivide the yes/no categorization. Day = indicates which day
 (starting in 1) of sampling the sample was collected at. Numbers for each
 year separatedly. Extraction = indicates which extraction run each sample
 was included in. Numbers for each year separatedly. sample_weight = weight
 of the sample before starting DNA extraction. If cell is empty, sample was
 not weighted. Data not used in the analysis. 4\. Missing data codes:
           None 5\. Abbreviations used:           NA; not applicable
 DATA-SPECIFIC INFORMATION FOR: ngs_sper01_ovp1.txt + ngs_sper01_ovp2.txt +
 ngs_sper01_ovp3.txt + ngs_sper01_ovp4.txt 1\. Number of variables: 7 2\.
 Number of rows: 384 3\. Variable list: Column1: Experiment name Column2:
 Replicate name Column3: Forward tag : Reverse tag Column4: Forward primer
 Column5: Reverse primer Column6: "F" - Needed to run Obitools
 Column7: Position in the plate \##################### INSTRUCTIONS TO
 PROCESS RAW SEQUENCING OUTPUT 4 sequencing files, two files for each
 sampling year (Each year was extracted amplified and sequenced
 independently). 1\. Unzip raw files 2\. Join paired ends to have aligned
 dataset 3\. Use ngs_files for each year respectively to assign sequences
 to each individual sample. 4\. Use the metadata files to assign samples to
 each individual, year and and other variables of interest.