Isotope analysis combined with DNA barcoding provide new insights into the dietary niche of khulan in the Mongolian Gobi

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.


Introduction
Dietary niche partitioning is a key feature facilitating coexistence in multi species assemblages of wild ungulates. Competition can be reduced by partitioning plant resources according to different taxa (generalists versus specialists), parts (only leaves, fruits, buds versus the whole plant), development stages, or functional types (grass versus browse) in a shared ecosystem [1][2][3]. However, with ever-increasing domestic ungulate populations, competition and avoidance are increasingly shaping resource availability for wild ungulates. When resources become limiting, wild ungulates are likely to be outcompeted by domestic ungulates, as the latter often greatly outnumber the former and are favoured by human interventions [4,5]. As a result, wild ungulates have to make do with depleted and altered rangelands or move into marginal habitats [6]. Shifts in their dietary niche are likely and ultimately can be expected to negatively affect their fitness [7]. Especially when dealing with threatened ungulate species, characterizing and monitoring the dietary niche is an important tool to detect anthropogenic impacts and guide conservation planning [8].
The wild ungulate community in Mongolia has remained almost complete, but plains species have become largely confined to the least productive or most remote parts of the country, likely in response to hunting, disturbance, and competition with livestock [9][10][11]. Nowadays, these remote regions, particularly the Gobi-Steppe ecosystem, constitute global refuges for several threatened ungulates, but are no longer isolated from global economic forces [12]. Mongolia's rural economy has been based on semi-nomadic livestock husbandry for millennia and it still forms the backbone for the majority of rural communities. However, socio-political changes, with a transition from centrally planned to free-market economy, have resulted in dramatic livestock increases [13,14]. This has led to growing concerns over rangeland health and impacts on threatened ungulates, including the Asiatic wild ass or khulan (Equus hemionus; [9]).
With an estimated 64,000 individuals, the Mongolian Gobi houses >80% of the global khulan population, and the situation in the Gobi determines the species' global status and trend [15,16]. Khulan are believed to be primarily grazers, but seem to be able to include a wide spectrum of grasses, forbs, and shrubs into their diet depending on region and season [17,18]. However, systematic, multi-seasonal studies at the landscape scale have been scarce [19]. Sequential stable isotope analysis of wild equid tail hair from the Dzungarian Gobi (DG) in south-western Mongolia suggested that khulan primarily graze in summer, but switch to a suboptimal mixed grass / shrub diet in winter, whereas sympatric domestic horses (Equus ferus caballus) and reintroduced Przewalski's horses (Equus ferus prezewalskii) don't show this switch and graze year round [20].
Winter is the time of high energy expenditure for thermoregulation and in the DG coincides with the return of herders and their livestock from their summer ranges in the Altai mountains onto the wild equid range on the Gobi plains. The most likely interpretation for the switch in the khulan's winter diet is displacement from more productive grass-dominated areas by humans and their livestock due to disturbance and persecution by humans directed at khulan, but not at charismatic Przewalski's horses [9,20,21]. In line with this hypothesis, one would expect higher proportions of low-quality browse in the diet in areas where khulan share the range with livestock camps year-round, or where the density of livestock camps is higher. This is the case in the South Gobi Region (SGR) in south-eastern Mongolia, where the absence of large mountain ranges negates the benefit of seasonal altitudinal movements while proximity to urban markets makes livestock herding more profitable.
Stable isotope analysis allows differentiation between isotopically contrasted food groups, but does not provide information at the plant species level [22]. Correct interpretation of results hence requires a priori knowledge of the most likely food plants in a given study area [22,23]. The plant community in the DG allows for a separation of grasses versus shrubs because C4 grasses (Poaceae, which consist of C4 as well as C3 species) are virtually absent, as they reach their northern distribution limit due to the cold and dry climatic conditions [24]. At the same time, C4 shrubs from the family Amaranthaceae persist and dominate the drier regions due to their high water retention and storage potential, deep roots, and hardy shrub lifeform [20,24]. In the SGR, the climate is slightly warmer and moister and allows for the additional growth of C4 grasses, potentially compromising the suitability of using stable isotope analysis to assess the dietary niche of grazers along the grazing / browsing continuum.
Recent developments in DNA barcoding-taxon identification using a conserved and standardized DNA region-and massively parallel sequencing have allowed for rapid and cost-efficient analysis of large numbers of fecal samples for diet analysis at the species, genus, or at least family level [25][26][27]. The DNA barcoding method overcomes many of the restrictions of classical behavioural and micro-histological studies and can provide unprecedented taxonomic details of diet items in fecal samples [28,29].
In this study, we combined isotope analysis of khulan tail hair (allowing inferences about multi-year diet seasonality and isotopic dietary niches) and fecal samples (providing a snapshot of the diet over the last 1-2 days; [30]) with DNA barcoding of fecal samples (allowing accurate plant taxa identification). Our overall goal was to achieve a more comprehensive understanding of the diet switch previously observed during the winter energetic bottleneck in one of two study areas in the Mongolian Gobi. Our two primary aims were to: i) verify previous results and assess whether the same analytical framework can be used for ungulates throughout the Gobi region and ii) look for support that the diet shift in winter in the DG is linked to seasonal differences in degrees of human and livestock presence. Specifically, we expected to find: 1. Based on stable isotope analysis of tail-hair samples-a seasonal shift from a grass dominated C3 summer diet to a shrub dominated C4 winter diet in the SGR in line with previous results from the DG; 2. Based on the DNA barcoding results-confirmation of the switch from a grass dominated summer diet to a mixed grass/browse diet in winter via taxonomic identification of plant remains in winter feces from the DG, but a less clear separation in winter feces from the SGR due to a climate less exclusive to C4 grasses; 3. Based on a comparison of both methods for fecal samples-high consistency between isotope signatures and DNA barcoding in respect to the proportion of C3/C4 plants identified in the diet; 4. Based on combined results from both methods-higher levels of C4 shrubs in the diet of khulan in the SGR due to the higher density and year-round presence of herder camps and the resulting higher potential for displacement by and competition with humans and their livestock.
We believe this study provides a valuable framework for future studies of ungulates in drylands and provides important empirical insights into the potentials and limitations of dietary analysis methods which are relevant for future research in general, and ungulate conservation in the Gobi specifically.

Study area
The two study areas are located at the eastern and western end of the khulan distribution in the Mongolian Gobi (Fig 1). Towards the south and west, the khulan distribution range in the Gobi is limited by the fenced international border towards China, towards the east by the fenced Trans-Mongolian railway, and towards the north by increasing human and livestock density [31]. The western study area in the DG covers roughly 9,000 km 2 and falls entirely within Great Gobi B Strictly Protected Area which was extended from 9,000 km 2 to 18,000 km 2 in 2019. The eastern study area covers roughly 56,000 km 2 in the SGR. Both study areas fall primarily in the semi-desert zone, but vary somewhat in respect to topography, climate, vegetation, khulan population size, and lands use ( Table 1, S1 File).

Sample collection and preparation
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  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. Hair

Stable isotope analysis
The stable isotope analysis was conducted at the Leibniz Institute for Zoo and Wildlife Research, Berlin, Germany following methods described in ( [20,39] and S1 File). The precision of measurements was always better than 0.1 ‰ for δ 13 C and δ 15 N, and 1.0 ‰ for nonexchangeable δ 2 H values, based on repeated analysis of laboratory standards, calibrated with international standards. In the context of this study, δ 2 H values were primarily used to assign the δ 13 C profiles to the correct seasons [39].
We used a Bayesian approach to estimate species specific core isotopic dietary niches based on bivariate, ellipse-based metrics (for δ 13 C and δ 15 N values) using SIBER implemented in the R package SIAR (Version 4.2); see S1 File 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˚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 (ITS1-Poa; S4 Table). For details on PCR and library purification and pooling (S2 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 (S1 Fig) and running a 2 nd subsample for six feces from the SGR to get a feeling for within-fecal sample variation (S2 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 S3 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 (S5 Table). Our cutoff 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.
The method of calculating diet proportions influences the results of diet characterization. Initial studies suggest that using RA introduces biases particularly for common diet taxa [45] while using absolute or relative FOO of taxonomic units is sensitive to small samples sizes and likely overestimates the importance of rare items or abundant species accidentally ingested together with the target species [46]. To allow comparison with previous studies we chose to present diet for each region at the genus level as: 1) RA, representing the proportion of sequence reads belonging to a unique genus in a fecal sample, divided by the total number of sequence reads in that fecal sample and 2) FOO, representing the number or percentage of feces in which a unique sequence read for a genus was found with a frequency >0.005. We additionally derived relative frequency of occurrence (rFOO) by dividing the sum of the FOO of each genus by the sum of the FOO of all genera (and families for comparison with [47]). We used the Bray-Curtis index, scaled from 0 (complete diet overlap) to 1 (no diet overlap), calculated with vegdist in the R package vegan [48] to compare diet dissimilarities among samples and sampling sites.

Seasonal isotope pattern of khulan tail hair and feces
Isotope values of khulan tail hairs from the SGR did not show the expected diet switch from high δ 13 C diet values in winter (up to 67% C4 plant intake) and low values (<10%) in summer as previously shown for the DG in [20]. Rather, values followed a bimodal pattern with highest δ 13 C diet values (high C4 plant intake of up to 84%) in mid-summer and mid-winter and lower values in fall (<40%) and spring (<20%; Fig 2).
Isotope values of khulan feces from winter showed only a slightly higher mean intake of C4 plants in the SGR (40.9%) as compared to the DG (37.5%, Table 2

Winter diet based on barcoding of feces
Using FOO, we identified 39 genera from 11 families, of which 27 genera from 6 families were present in feces from the DG and 29 genera from 10 families in feces from the SGR (Fig 3, S5  Table). Using RA the number of genera was reduced to 10 genera from 4 families for the DG (representing 96% of all reads) and 18 genera from 7 families for the SGR (representing 98% of all reads; Fig 3).
Diet calculations based on FOO and RA resulted in different proportions given to diet items at the genera level but agreed on the bigger picture. Poaceae, Amaranthaceae, Tamaricaceae, and Asteraceae were identified as the most important plant families in the winter diet of khulan in the DG and the SGR. The importance of Tamaricaceae (represented in the diet by the single genus Reaumuria with one species in the Gobi) is undisputed but may be somewhat overestimated because of an amplification bias in favour of the Reaumuria (for details see S1 Fig). Six additional plant families were of minor importance, with four only found in the SGR and one only found in the DG (Figs 3, 4, S5 Table).
The diet in the SGR was more diverse and, contrary to our expectation, it was dominated by Poaceae, while the diet from the DG was a more even mix of Amaranthaceae, Poaceae, Tamaricaceae, and Asteraceae. Reaumuria, Stipa (Poaceae), and Haloxylon (Amaranthaceae) were the most important genera in the DG and Eragrostis (Poaceae) and Reaumuria the most important genera in the SGR (Fig 3).
Genera with a C4 photosynthetic pathway were more important in the SGR as compared to the DG (Table 3, Fig 3). However, contrary to our expectation, C4 genera in the SGR were primarily recruited from Poaceae rather than from Amaranthaceae like in the DG. In general, FOO suggested a higher degree of similarity in diet composition between the two study areas than RA (Table 3, Fig 3).
Despite a larger number of genera found in the feces, diet dissimilarities among feces based on RA were only slightly higher in the SGR as compared to the DG (mean Bray-Curtis dissimilarity = 0.35 and 0.31 respectively; low values = high overlap or low dissimilarity). Feces were more similar within sampling location than among all locations pooled, with the effect being stronger in the SGR as compared to the DG (weighted mean Bray-Curtis dissimilarity = 0.11 and 0.25 respectively ; Fig 4).
There was an overall good agreement in the proportions of C4 plants in the winter diet of khulan based on isotope analysis (tail hair and feces) and barcoding (feces) with values being consistently higher for the SGR as compared to the DG (Table 2).

Concordance of barcoding results with past isotope analysis
Our barcoding results revealed an unexpected diet composition in the SGR, but they also largely confirmed the interpretation of previous isotope analysis of tail hair from the DG [20]. Khulan in the DG feed on a mixed grass / shrub diet in winter, the consumed grasses are largely composed of C3 grasses from the genus Stipa, and a large proportion of consumed shrubs is composed of C4 shrubs, namely Haloxylon ammodendrum, Anabasis brevifolia, and Salsola from the family Amaranthaceae. New insights highlight the importance of genera in the Asteraceae family, providing further evidence that the importance of forbs may have been underestimated with conventional diet analysis, at least for equids living under semi-desert or desert conditions [49].
Genera identified with barcoding for the DG roughly matched results of micro-histological analysis by Xu et al. [9] in respect to the high importance of Stipa, Haloxylon, and Reaumuria, but differed greatly in respect to Anabasis (S6 Table). For rarer genera, differences vary considerably, and interestingly, [19] failed to detect any C4 grasses, which may point to the difficulty of finding rare diet items or distinguishing between anatomically similar species [29]. On the other hand, feces only provide a snapshot of a khulan's diet over the last 1-2 days [30] and fecal composition varies within and, even more so, between sample locations. Given the relatively small sample sizes of studies available so far and the method differences between barcoding and micro-histological studies, discrepancies for diet items of lower importance are hardly surprising and should not be overinterpreted [29,50].

Khulan diet seasonality in the South Gobi Region
At first glance, our isotope data seems to support the hypothesis that khulan in the SGR are forced to shift to a more C4 browse-dominated diet due to the higher density and year-round presence of herder camps [20]. The isotope analysis confirmed higher levels of C4 plants and identified a bi-modal seasonal diet shift in the SGR with C4 peaks in winter and summer. However, subsequent barcoding contested the validity of this conclusion as it identified C4 grasses as the main source of C4 plants in the winter diet of khulan in the SGR. As it turns out, genera with C4 grasses like Eragrostis, Cleistogenes, and Aristidia dominated the winter diet, rather than only complementing the consumption of the widespread and common C3 Stipa grass, which constitutes a staple for khulan in the DG. Furthermore, the shrubs consumed were composed of genera that included C4 and C3 species instead of being dominated by C4 species as in the DG. This mix of genera with C3 and C4 pathways in grasses and shrubs makes the unambiguous isotopic separation along the grazer-browser spectrum impossible for the SGR. Barcoding allowed us to unambiguously identify C4 grasses as the main source of C4 plants in the diet of khulan in the SGR in winter. We currently lack this information for the C4 summer peak, though C4 grasses seem the most likely candidate [47]. The C3 lows in spring are likely a consequence of the first emergence of C3 spring forbs like Allium and cold robust C3 grasses before increasing temperatures allow the re-emergence of C4 grasses [51]. The same may be happening with C3 and C4 grasses in fall, with the former being photosynthetically active longer than the latter and therefore preferred by grazers like khulan. However, additional seasonal studies are required to confirm an overall, year-round high use of C4 grasses. Moreover, if confirmed, further field investigations will be needed to test whether the high C4 grass use reflects their higher availability in areas where khulan roam or is the result of a preference over common C3 grasses like Stipa.

Unresolved question of diet shift in response to humans
Climatic and biogeographic changes from the western to the eastern Gobi result in such profound diet differences in khulan that predictions from one region cannot be tested in the other region within a simple comparative framework. Our results, highlighting a grass-dominated winter diet in the SGR, offers the tantalising assumption that this occurs despite disturbance by and competition with local herders and their livestock, and that one would expect that the summer diet in the SGR consists entirely of abundant grasses and forbs and excludes shrubs altogether. Unfortunately, we lack the seasonal barcoding data to confirm this hypothesis. Barcoding results from khulan summer feces, collected and analysed by [47] suggests that use of the shrub dominated plant families Amaranthaceae and Tamaricaceae is reduced, use of Poaceae remains high, and use of forbs increases and diversifies. The summer diet includes several new families with genera typical for spring and summer and absent in winter like Alium (Amaryllidaceae) and Iris (Iridaceae; S4 Fig). However, while there is some evidence of a diet shift in khulan from summer to winter, quantitative comparisons with our results are not possible because [47] used the trnL primer alone which results in a lower taxonomic resolution. Re-calculating diet composition for comparison as FOO at the family level, results in a great loss of taxonomic detail and hence is unsuitable for detecting more subtle differences or gradual shifts in diet. With 23-47% (depending on study area and method of calculation) the proportion of grasses in the winter diet of khulan in the Mongolian Gobi is generally low compared to other wild and free-ranging equids for which grasses are the primary diet items, typically amounting to 70-100% [18,25,49]. However, without a fully comparative framework which includes sympatric domestic species and nutritional analysis, it is difficult to draw inferences about whether these changes solely reflect seasonal availability of preferred plants or point towards competition with, or displacement by, livestock.

Methodological considerations
Our barcoding results showed good consistency with isotope analysis in respect to the proportion of C4 plants in the diet and also with the micro-histological study by [19]. Our study therefore adds to a growing body of evidence that barcoding and the use of RA is a promising quantitative tool, especially for comparative diet studies within the same bio-geographical region, where the main interest is diet dissimilarity or diet trends, rather than diet composition per se [25,46]. The latter is more difficult to obtain as results are influenced by primer selection, amplification bias (as shown also in our own small experiment), potential interactions with other food items, differences in the density of target DNA among plants and between different structural parts of the same plant, and different sensibility of plant matter to digestive breakdown [28,50]. However, alternative methods have their own challenges [29].
The comparative nature of our study and the use of an identical analysis framework accounts for most of the aforementioned constraints. We are therefore confident that our results allow for a valid comparison of the winter diets between the two regions of the Mongolian Gobi and provide sufficient taxonomic resolution for a more comprehensive understanding of the diet switch previously observed during the winter energetic bottleneck in the DZ. In addition, we identified a large spectrum of genera utilized by khulan, several of which had not been previously identified. The unexpected high use of C4 grasses in the SGR raises the question of whether khulan, which are generally considered bulk grazers, may actually be much more selective than previously assumed.

Conclusions
Our analysis approach of combining isotope analysis of tail hair and barcoding of feces 1) provided unprecedented taxonomic detail and documented previously unknown differences in the winter diet of khulan in two different parts of the Mongolian Gobi, 2) established that longitudinal isotope analysis of long tail hair in an ungulate is a suitable tool to detect and compare population level diet shifts, 3) confirmed that the specific climatic conditions in the DZ allow for an isotopic differentiation along the grazer-browser continuum, and 4) revealed that the same conditions are not met in the SGR.
Our results once again highlight the importance of parallel analysis approaches when using indirect analysis methods, especially when moving from intensively studied to less intensively studied or larger study areas. The fact that khulan feces were more similar within sampling locations than among sampling location, highlights: 1) the need for fecal sampling to be spread across a region to account for the influence of landscape heterogeneity, 2) a potential for using the habitat characteristics of the wider sampling location as covariates to explain diet differences.