DOI: 10.5061/dryad.5x69p8d2d
title: Isotope analysis combined with DNA barcoding provide new insights into the dietary niche of khulan in the Mongolian Gobi

description: With increasing livestock numbers, competition and avoidance are
 increasingly shaping resource availability for wild ungulates. Shifts in
 the dietary niche of wild ungulates are likely and can be expected to
 negatively affect their fitness. The Mongolian Gobi constitutes the
 largest remaining refuge for several threatened ungulates, but
 unprecedentedly high livestock numbers are sparking growing concerns over
 rangeland health and impacts on threatened ungulates like the Asiatic wild
 ass (khulan). Previous stable isotope analysis of khulan tail hair from
 the Dzungarian Gobi suggested that they graze in summer but switch to a
 poorer mixed C3 grass / C4 shrub diet in winter, most likely in reaction
 to local herders and their livestock. Here we attempt to validate these
 findings with a different methodology, DNA metabarcoding. Further, we
 extend the scope of the original study to the South Gobi Region, where we
 expect higher proportions of low-quality browse in the khulan winter diet
 due to a higher human and livestock presence. Barcoding confirmed the
 assumptions behind the seasonal diet change observed in the Dzungarian
 Gobi isotope data, and new isotope analysis revealed a strong seasonal
 pattern and higher C4 plant intake in the South Gobi Region, in line with
 our expectations. However, DNA barcoding revealed C4 domination of winter
 diet was due to C4 grasses (rather than shrubs) for the South Gobi Region.
 Slight climatic differences result in regional shifts in the occurrence of
 C3 and C4 grasses and shrubs, which do not allow for an isotopic
 separation along the grazer-browser continuum over the entire Gobi. Our
 findings do not allow us to confirm human impacts upon dietary preferences
 in khulan as we lack seasonal samples from the South Gobi Region. However,
 these data provide novel insight into khulan diet, raise new questions
 about plant availability versus preference, and provide a cautionary tale
 about indirect analysis methods if used in isolation or extrapolated to
 the landscape level. Good concordance between relative read abundance of
 C4 genera from barcoding and proportion of C4 plants from isotope analysis
 adds to a growing body of evidence that barcoding is a promising
 quantitative tool to understand resource partitioning in ungulates.
description: <strong>Stable
 isotope analysis</strong> For multi-year diet analysis, we
 plucked several long tail hairs from adult khulan during capture and
 anaesthesia when GPS collars where attached to khulan as part of ongoing
 monitoring programs. For analysis, we selected individuals with the
 longest tail hairs for maximal temporal coverage, namely 3 females and 3
 males captured in July 2009 in the DG (Fig 1, S2 Table; (20)) and 4
 females and 4 males captured in August 2013 in the SGR. Hair samples were
 stored in labelled paper envelopes until arrival in the laboratory where
 they were cleaned and cut into 1 cm increments for sequential CN or H
 analysis, as reported in (39). For winter diet analysis, we
 collected fresh fecal samples from khulan in February 2016: 42 samples at
 6 different locations in the DG and 25 samples at 8 locations in the SGR
 (Fig 1, S3 Table). Samples were collected in plastic bags and kept cool
 before and during transport to the laboratory where outer layers were
 carefully removed to avoid contamination with other DNA caused by contact
 with the environment and to reduce endogenous DNA proportions. The
 remaining material was dried to constant weight at 50° C and ground using
 a mortar and pestle for isotope analysis and barcoding. The stable isotope analysis was
 conducted at the Leibniz Institute for Zoo and Wildlife Research, Berlin,
 Germany following methods described in ((20, 39) and S4 File). The
 precision of measurements was always better than 0.1 ‰ for
 <i>d</i><sup>13</sup>C and
 <i>d</i><sup>15</sup>N, and 1.0 ‰ for
 non-exchangeable <i>d</i><sup>2</sup>H values,
 based on repeated analysis of laboratory standards, calibrated with
 international standards. In the context of this study,
 <i>d</i><sup>2</sup>H values were primarily used
 to assign the <i>d</i><sup>13</sup>C profiles to
 the correct seasons (39). To convert raw hair isotope values
 (<i>d</i><sup>13</sup>C<sub>hair</sub>
 and
 <i>d</i><sup>15</sup>N<sub>hair</sub>)
 into diet values
 (<i>d</i><sup>13</sup>C<sub>diet</sub>
 and
 <i>d</i><sup>15</sup>N<sub>diet</sub>), we used the dual mixing model of (40). To calculate the % of C4 diet in the feces, we used a dual end-member mixing model from (41), where % C4 = (<i>d</i><sup>13</sup>C<sub>C3 plants </sub>+ ∆<i>d</i><sup>13</sup>C - <i>d</i><sup>13</sup>C<sub>feces</sub>) / (<i>d</i><sup>13</sup>C<sub>C3 plants</sub> -∆<i>d</i><sup>13</sup>C<sub>C4 plants</sub>). We used a Bayesian approach to estimate species specific core isotopic dietary niches based on bivariate, ellipse–based metrics (for <i>d</i><sup> 13</sup>C and <i>d</i><sup> 15</sup>N values) using SIBER implemented in the R package SIAR (Version 4.2); see S4 File in the original manuscript for more details on plant sampling, stable isotope analysis, and visualization. DNA analysis of fecal samples, species identification and plant metabarcoding DNA extraction and amplification Total genomic DNA was extracted from 100 mg of powdered fecal material using the FastStool kit (Qiagen, Germany) according to the manufacturer’s instructions but with some slight modifications. Specifically, the incubation with Buffer AL was performed at 95° C (rather than 70° C) for 10 minute and the final column elution step was carried out in 2 x 20 µl volumes for 10 minutes each. Extracts were stored at -20<sup>o</sup>C until required for downstream analyses. Dietary NGS library preparation followed the protocol of (25) using three primer sets to broadly cover all plant orders (trnL(UAA)g-h) while providing enough taxonomic resolution to differentiate genera from Asteraceae (ITS1-Ast) and from Poaceae<i> </i>(ITS1-Poa; S5 Table). For details on PCR and library purification and pooling (S6 File). Sequencing and data analysis Sequencing was performed as 125 basepair paired end (PE) on an Illumina HiSeq 2000 at the Vienna Biocenter Core Facility. Data were delivered as BAM files and converted to fastq format using the bam2fq command in SAMtools (42). Demultiplexing was performed with the command barcode_splitter (FastX toolkit, Hannon Lab). Paired reads were merged using the tool FLASH with the “innie” and “outie” options implemented to ensure that short reads that include sequencing adaptors are also incorporated. Sequencing adaptors and primers were then removed using the tool Flexbar before quality filtering to remove merged reads with greater than 20% of bases with quality scores less than Q30 (FastX toolkit, Hannon Lab). We also included two small experiments, sequencing 50:50 mixes of the five most important plants from the DG to test for amplification bias (S7 Fig.) and running a 2<sup>nd</sup> subsample for six feces from the SGR to get a feeling for within-fecal sample variation (S8 Fig.). Reads were processed using the UParse pipeline (43) to create clusters of operational taxonomic units (OTUs) showing &gt;97% similarity while simultaneously removing chimeric sequencing artefacts (all relevant scripts can be found in S13 File). Taxa was assigned to OTUs by blasting sequences against the NCBI Viridiplantae sequence database within the program CLC-Main Workbench (Qiagen, Germany). Settings were an Expect value of 10.0, Word size of 11, Match value 2, Mismatch value -3, gap costs of Existence 5 and Extension 2. Maximum number of hit sequences was set to 2. Taxonomy was accepted only if identity was &gt;95% and abundance was greater than 0.01% overall. All species sequence reads were collapsed down to genus level. We were not able to collect and sequence all potential plant species or genera over the huge expanse of the two study areas and did not expect all plant genera to be available in the NCBI Viridiplantae sequence database. Therefore, rather than <i>a priori</i> removing genera not native to Mongolia or the Gobi and potentially introducing a bias towards better known plant genera, we used a sequence read cut-off value of 0.5% over all genera reads by region (for read abundance RA) and individual scat (for frequency of occurrence FOO). The remaining sequence reads were checked for their presence in Mongolia and the Gobi based on (S10 Table). Our cut-off value removed all but one genus not found in Mongolia/Asia and two without documented records from the Gobi; all three were only of minor importance. We assume that there must be an unidentified closely related species in our study area or a misclassification sequence in the database (44). Whether or not we included these three genera made little difference to the overall results; the chose to include them which gave a slightly higher importance and diversity to the Asteraceae family.
description: Overview available files: 1) EXCEL file
 "2020_05_09_Barcoding_Khulan_read_counts_taxa_assigned_dryard.xlsx": Folder "DG_winter_scats":  Contains the number of barcoding reads for different plan taxa at genus or species level for 42 dung samples of khulan (<em>Equus hemionus</em>) from the Dzungarian Gobi in SW Mongolia collected in winter 2016. Folder "SGR_winter_scats": Contains the number of barcoding reads for different plan taxa at genus or species level for 24 dung samples of khulan (<em>Equus hemionus</em>) from the South Gobi Region in SE Mongolia collected in winter 2016. 2) EXCEL file "Dryard_isotope_data_BurnikSturm_Smith_etal_2021_dryard.xlsl": Folder "Raw_isotope_data_TAIL_HAIR": Contains the raw d13C<sub>hair</sub> (‰), d15N<sub>hair</sub> (‰), and d2H<sub>hair</sub> (‰) values of khulan tail hair aligned by date. Six of the khulan are from Dzungarian Gobi in SW Mongolia and were collected in 2009, and the remaining 8 khulan from the South Gobi Region in SE Mongolia and were collected in 2013. Folder "Raw_isotope_data_FECES": Contains raw d13C<sub>feces</sub> (‰) and d15N<sub>feces</sub> (‰) values and sample locations of 42 fecal samples of khulan in the Dzungarian Gobi in SW Mongolia and 24 dung samples of khulan from the South Gobi Region in SE Mongolia; all samples were collected in winter 2016. Folder "Raw_isotope_data_Plants": Contains raw d13C (‰), %C, d15N (‰), %N values and sample locations of 8 <em>Stipa</em> grass (a dominant C3 grass) samples and 4 <em>Haloxylon</em> (an important C4 shrub/tree) samples from the South Gobi Region in SE Mongolia , for comparison with values previously obtained for the two genera from the Dzungarian Gobi in SW Mongolia.