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.

Stable isotope analysis

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 d¹³C and d¹⁵N, and 1.0 ‰ for non-exchangeable d²H values, based on repeated analysis of laboratory standards, calibrated with international standards. In the context of this study, d²H values were primarily used to assign the d¹³C profiles to the correct seasons (39).

To convert raw hair isotope values (d¹³C_hair and d¹⁵N_hair) into diet values (d¹³C_diet and d¹⁵N_diet), 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 = (d¹³C_{C3 plants} + ∆d¹³C - d¹³C_feces) / (d¹³C_{C3 plants} -∆d¹³C_{C4 plants}).

We used a Bayesian approach to estimate species specific core isotopic dietary niches based on bivariate, ellipse–based metrics (for d ¹³C and d ¹⁵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^oC 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 (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^nd 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 >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 >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 a priori 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.

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 (Equus hemionus) 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 (Equus hemionus) 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_hair (‰), d15N_hair (‰), and d2H_hair (‰) 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_feces (‰) and d15N_feces (‰) 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 Stipa grass (a dominant C3 grass) samples and 4 Haloxylon (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.