Molecular dietary analysis of two sympatric felids in the Mountains of Southwest China biodiversity hotspot and conservation implications

Dietary information is lacking in most of small to mid-sized carnivores due to their elusive predatory behaviour and versatile feeding habits. The leopard cat (LPC; Prionailurus bengalensis) and the Asiatic golden cat (AGC; Catopuma temminckii) are two important yet increasingly endangered carnivore species in the temperate mountain forest ecosystem in Southwest China, a global biodiversity hotspot and a significant reservoir of China’s endemic species. We investigated the vertebrate prey of the two sympatric felids using faecal DNA and a next-generation sequencing (NGS)/metabarcoding approach. Forty vertebrate prey taxa were identified from 93 LPC and 10 AGC faecal samples; 37 taxa were found in the LPC diet, and 20 were detected in the AGC diet. Prey included 27 mammalian taxa, 11 birds, one lizard and one fish, with 73% (29/40) of the taxa assigned to the species level. Rodents and pikas were the most dominant LPC prey categories, whereas rodents, pheasant, fowl and ungulates were the main AGC prey. We also analysed the seasonal and altitudinal variations in the LPC diet. Our results provide the most comprehensive dietary data for these felids and valuable information for their conservation planning.

high population densities relative to other carnivores in this ecosystem based on camera-trapping surveys 9,12 and molecular species identification of faecal samples (see Materials and Methods). However, very little data are available regarding the ecological functions of these cats or their trophic interactions with other species.
LPC is a small felid with an adult body weight range of 1.7-7.1 kg 13 . It is one of the most common wild cats with a wide distribution across Asia and is found in a variety of habitats 14 . The population status of the LPC varies in different regions. The species is listed as least concern on the IUCN Red List 15 , although some local populations and subspecies are endangered 16 . In China, LPC was extensively harvested for fur trade prior to the early 1990 s, and its populations in the South and Southwest China may have suffered significant reductions 17 . Current threats to the species are mainly habitat loss and fragmentation due to anthropogenic land use 8 . No systematic population survey of the LPC has been conducted in China. Although LPC is detected at a relatively high frequency by camera-trapping surveys in the study area 18 , the population status of the species is unclear. Analyses of faeces and stomach contents show that the main LPC prey items are rodents (mostly murids), supplemented with other small mammals, such as squirrels and shrews, birds, reptiles, amphibians and fish, across its range [19][20][21][22][23][24][25][26] . Little data are available on their dietary habits in most regions of China.
AGC is a medium-sized felid (adult body weight: 9-16 kg) mostly distributed in Southeast Asia and South and Central China, primarily found in forest habitats 17,27 . Due to habitat loss and declining prey abundance, AGC populations are decreasing in most of its range 28 and this species is listed as near threatened on the IUCN Red List 15 . AGC population status in China is unknown. In the study area, AGC occurs at a much lower density than LPC, according to camera-trapping records 9 . Data on AGC dietary habits are scarce. Morphological analysis of faeces or stomach contents of AGC in Thailand and Malaysia have identified small mammals, such as rodents and ground squirrels, small ungulates, birds, reptiles and amphibians 21,29,30 . However, AGC diets in China remain uninvestigated.
Conventional dietary analysis methods are primarily based on microscopic examination of hard parts of food remains in faeces or stomach/gut contents. These methods rely highly on undigested portions of prey, expert knowledge of their morphology, are generally time-consuming and often unable to provide precise identification of species 31 . Molecular identification of species based on sequence analysis of polymerase chain reaction (PCR) amplified fragments of a standardised DNA region (i.e. DNA barcoding 32 ) from faecal or gut samples provides a reliable and effective alternative to dietary analysis in a wide range of species. This type of analysis requires no a priori information on prey composition, is powerful for resolving closely-related species and highly complex diets, is highly efficient for processing large-scale samples and the results can be rechecked and refined at a later time with updated sequence information 31,33 . Recent technological advances in next-generation sequencing (NGS) combined with metabarcoding have made DNA-based dietary analysis even more efficient and affable 34 .
We previously analysed the LPC diet in this area using faecal DNA and a traditional cloning-sequencing approach and identified 16 vertebrate taxa in 25 faecal samples 35 . However, due to the small sample size and low throughput sequencing technique, the results are unlikely to represent the complete dietary profile of the species, nor demonstrate the seasonal or regional variations. The goals of this study were (1) to investigate vertebrate prey compositions of LPC and AGC and analyse potential seasonal and altitudinal differences in the diet in the mountain forests of Southwest China using an NGS and metabarcoding-based approach and (2) to inform conservation planning and management actions based on the dietary data.
We collected faecal samples in the Reserve from March 2013 to November 2014. The majority of samples were collected in spring (March-May) and autumn (September-November). We were unable to obtain samples during the summer season (July-August) or winter months (December-February) due to weather conditions. Mammalian faeces were collected along fixed transects at altitudes of 1,250-3,200 m ( Fig. 1) during weekly patrols, submerged in 95% ethanol and stored at − 20 °C. DNA extraction. We used the 2CTAB/PCI method to extract DNA from faecal samples 37 . We included approximately 100 mg of the external surface in each DNA extraction for molecular identification of species. Faecal samples that were confirmed to be from LPC or AGC were homogenised individually, and 100 mg of the homogenised faeces was used in DNA extraction for the molecular dietary analysis.

Molecular species identification.
We designed the 16S-F/16S-R primer pair to amplify a ~350 bp fragment of the mitochondrial 16 S rRNA gene to unambiguously assign all known predator species in the study area 35 . The LPC and AGC sequences of this fragment differed at 15 nucleotides (GenBank Accession numbers LPC: JN392459.1 and AGC: KP202267.1). The PCR conditions and sequencing of the products were described previously 35 . The species was assigned when faecal DNA amplicons matched the sequences in the GenBank database with 98-100% identity.
We selected LPC and AGC faecal samples to be included in subsequent dietary analyses based on DNA quality and sample location information. Samples that failed to yield sufficient DNA for PCR amplification of food items were excluded. We also eliminated faecal samples collected in close proximity (< 50 m), which may represent pellets of a single defecation. A total of 94 LPC and 11 AGC faecal samples met the selection criteria and were included in the NGS run.
PCR for diet analysis. To amplify DNA of vertebrate species in the faecal DNA extracts by PCR, we used the 12SV5-F/12SV5-R primers targeting a ~100 bp fragment of the mitochondrial 12 S rRNA gene, which has been demonstrated to have high resolution power for identifying the genus and species across most vertebrate taxa 38 . A blocking oligonucleotide (PrioB) 39 designed to specifically inhibit amplification of the LPC sequences was also included in the PCR. We tested the blocking efficiency of PrioB with AGC samples and found that it also limited amplification of the AGC sequences at high efficiency (Appendix S1). The PrioB sequence diverged from small mammals (pikas, rodents and shrews), ungulates and birds that commonly occur in the study area (Appendix S2) and should not prevent amplification of these vertebrate groups.
PCR amplifications were conducted in a total volume of 20 μ L, including 1 × EasyTaq PCR SuperMix (TransGen Biotech, Beijing, China), 0.2 μ M 12SV5-F/R primers, 4 μ M PrioB, 0.4 mg/mL bovine serum albumin and 20-40 ng DNA extract. The PCR program started with an initial denaturation step of 5 min at 95 °C, followed by 35 cycles of 30 s at 95 °C, 30 s at 60 °C and no elongation. The 12SV5-F and 12SV5-R primers were tagged at the 5′ end with nine nucleotides that started with CC followed by seven variable nucleotides 40 . These tags differed by at least three nucleotides among the tags; hence, allowing a unique tag for each PCR to sort sequences according to samples following NGS. Seven negative PCR controls were also included in the amplifications to check for contamination.
Next-generation sequencing. PCR products were purified using a PCR purification kit (EasyPure PCR Purification Kit, TransGen Biotech). DNA concentrations were determined by agarose gel electrophoresis and using a NanoDrop 2000 spectrophotometer (ThermoScientific, Waltham, MA, USA) and were mixed in equimolar concentrations. Sequencing was carried out on the Illumina HiSeq 2500 (Illumina Inc., San Diego, CA, USA) at the NGS facility of the Biodynamic Optical Imaging Centre, Peking University, using paired-end sequencing following the manufacturer's instructions. A total of 100 nucleotides were sequenced on each end of the DNA fragment.
NGS data analysis. The sequence reads were first checked using the FastQC program to evaluate data quality, and low quality reads were filtered with the NGSQCToolkit v2.33.
We used the OBITools program 41 to analyse the sequence reads. The direct and reverse sequences were aligned using the illuminapairedend program, and aligned sequences with a quality score < 40 were removed using the obigrep program. Sequences with perfectly matched tags and a maximum of two mismatches in primers were identified using the ngsfilter program and kept for further analysis. Identical sequences were clustered into unique sequences using the obiuniq program. Sequences < 80 bp or with a total count in the whole dataset < 1,000 were removed using the obigrep program. PCR and sequencing errors were detected using the obiclean program 42 . A reference database was prepared by extracting vertebrate 12 S gene sequences from GenBank with the 12SV5 primers using the ecoPCR program 43 . Taxonomic identification was conducted using the ecoTag program 44 based on sequence similarity in this reference database. A unique taxon was assigned to each sequence. Sequence alignment and automatically assigned taxa were checked manually. We excluded sequences with < 94% identity to the query sequence to increase accuracy of the taxonomic assignments. We removed obvious human contaminants, sequences of the two felids and low-frequency sequences (< 0.1% or 50 reads in the sample or less than the count of the sequence in the negative control PCRs) that were the likely result of cross-contamination and/ or tag jumps 45 . The unique sequences were further collapsed into discrete taxa using a 2% sequence divergence threshold and relative abundance of the sequences. We refined the automatic taxonomic assignments with species distribution data from the study area 14,36,46 .
We used the following criteria for taxonomic assignments: (1) if the query sequence matched a single locally occurring species in the databases with ≥ 98% identity, a species was assigned; (2) if the query matched more than one species with ≥ 98% identity, a species was assigned based on knowledge of the distribution or the lowest taxonomic level that included all of the species was assigned if more than one species occurred locally; and (3) if the query matched the reference sequence with a maximum identity of < 98% but ≥ 94%, the query was assigned to the lowest taxonomic level that included all locally occurring species with the highest identity scores. If a single species showed the highest identity but was not known to inhabit the area, the most closely related local species within the same genus was assigned. Locally occurring species were identified using guides to the mammals 14 and birds 46 of China, and a survey of species in the Laohegou Nature Reserve 36 that included (but was not limited to) mid-to large-sized mammals, birds, amphibians, and fish. To avoid misidentification of taxa due to insufficient local species records, we kept to conservative taxonomic assignments and only excluded species that showed no occurrence in all of Southwest China.
Diet analysis. The LPC and AGC diets were quantified in two ways: (1) As a percent frequency of occurrence in faecal samples (%FC): where N i is the number of faecal samples containing the ith food resource, and N is the total number of faecal samples containing food; (2) As a proportion of occurrence of an individual food resource (%TX): The %TX of all food resources summed to one, but the sum of %FC of the food resources could be larger than one due to the presence of multiple food items in the faecal sample.
Calculation of the dietary parameters including the standardised Levins' measure of niche breadth (B A ), Shannon's diversity index (H), Peilou's J and Pianka's measure of niche overlap (O jk ) followed the equations shown in Appendix S3.
We analyzed the altitudinal variation in the LPC diet by dividing the entire altitudinal range (1250-2300) into four classes: low (< 1,500 m), lower-middle (1,500-2,000 m), higher-middle (2,000-2,500 m), and high (> 2,500 m). The definition of altitude classes was mainly based on vegetation type, as the distribution of forest types shows a strong altitudinal pattern. Mainly evergreen forests are observed at < 1,500 m, evergreen and deciduous mixed forests occur at 1,500-2,500 m, and deciduous and coniferous mixed forests and coniferous forests occur at 2,500-3,200 m. Since vegetation composition and structure also exhibits gradual changes along altitude gradients within each vegetation zone, we further divided the middle zone (1,500-2,500 m) into lower-(1,500-2,000 m) and higher-(2,000-2,500 m) classes to achieve relatively even ranges of the altitude classes.
We used the Wilcoxon signed-rank test to evaluate the seasonal difference of the main food categories in the LPC diet. Statistical significances of altitudinal variation in the LPC diet were assessed for each main prey taxon across the altitude classes with Fisher's exact tests. We estimated the 95% Wilson score intervals of the frequency of occurrence for the prey orders with the "Hmisc" package 47 in R software 48 .
Statistical analyses were conducted using R software 48 . All comparisons were two-tailed, with the significance level set to 0.05.

Results
NGS data analysis. The sequencing run generated 4.05 Gb of sequence data (40.1 million raw sequences), which was assembled to produce 18.6 million properly assembled (paired) sequences. After removing low quality sequences, primer or tag sequences and unidentified reads, we collapsed the sequences into 398,588 unique sequences. A further step to eliminate short (< 80 bp) or low-occurrence (< 1,000) sequences resulted in 638 sequences, which were reduced to 205 sequences after removing PCR and sequencing errors. The high-quality sequences were further reduced by removing low (< 94%) identity, low per-sample frequency and non-prey sequences and were collapsed into 40 discrete taxa. Based on the species distribution data, we assigned 29 taxa Scientific RepoRts | 7:41909 | DOI: 10.1038/srep41909 (73%) to species, four to genera, five to family or subfamily and two to order (Table 1). Of the 94 LPC and 11 AGC PCR products, 93 and 10 produced sequence reads, respectively and were included in the following analyses.
Dietary composition and diversity. Among the total 40 prey taxa, 37 were consumed by LPC and 20 by AGC, with 17 taxa shared by both species. The number of prey taxa per faecal sample varied from one to seven with a mean of 4.1 (SD 1.7) and 3.6 (SD 1.8) for LPC and AGC, respectively. The LPC diet contained 24 mammalian taxa (six orders), 11 bird taxa (three orders), one lizard and one fish, whereas the AGC diet was comprised   of 15 mammalian taxa (six orders) and five bird taxa (two orders) ( Table 1). The dietary niche analysis at the taxon level showed that LPC had a narrower diet breadth than that of the AGC, as indicated by Levins' standard- The most frequently occurring taxon in the LPC diet was pikas (Ochotona spp.), which was present in 76% of the faecal samples, followed by several rodents, including the Chinese white-bellied rat (Niviventer confucianus; %FC = 62%), South China field mouse (Apodemus draco; 44%) and Sichuan Chinese vole (Eothenomys chinensis; 40%). Another four taxa were detected in > 10% of the LPC samples: undetermined carnivorous species, 29%; large white-bellied rat, Niviventer excelsior, 27%; Père David's vole, Eothenomys melanogaster, 19%; and Temminck's tragopan, Tragopan temminckii, 16% (Table 1). Species in the order Rodentia were the most common prey of the LPC when the samples were analysed by the proportion of occurrence of individual prey taxa (%TX = 51%), followed by Lagomorpha (18%) (Fig. 2a).

Seasonal and altitudinal variations in
The LPC samples were collected at altitudes of 1,250-3,200 m. Species from the seven most common prey orders (Rodentia, Lagomorpha, Carnivora, Galliformes, Passeriformes, Soricomorpha, and Artiodactyla) were detected across the entire altitudinal range (low: < 1,500 m; lower-middle: 1,500-2,000 m; higher-middle: 2,000-2,500 m; and high: > 2,500 m) (Appendix S4). Fisher's exact tests showed that the frequency of occurrence (%FC) for Rodentia, Lagomorpha, and Carnivora varied significantly in different altitudinal zones (all p < 0.001), whereas altitudinal variations in the occurrence of Galliformes, Passeriformes, Soricomorpha, and Artiodactyla were insignificant. Calculation of 95% confidence intervals for Rodentia, Lagomorpha, and Carnivora indicated distinct distribution patterns for different prey orders (Fig. 4). For example, the frequency of occurrence of rodents was high at low and middle altitudes and reduced by half at altitudes > 2,500 m, whereas the occurrence of pikas and carnivores increased with altitude (Fig. 4). At the taxon level, various prey taxa also exhibited different patterns of occurrence (Fig. 5).  49 , raising questions about the prey selectivity of LPC. Additionally, pikas have been reported as LPC prey in our study area, but not in other regions, even where they are common food of other similar-sized carnivores 50 . Prey composition and availability along altitudinal gradients have been demonstrated to shape the diets of several carnivorous species 51,52 . Studies in Tangjiahe Nature Reserve showed that Niviventer and Apodemus spp. were the most abundant groups below 2,500 m in various forest types, but their numbers diminished quickly above this altitude; pikas were not detected below 1,800 m, occurred at a low number at 1,800-2,500 m, yet prevailed above 2,500 m in Alpine meadows and scrub habitats 53 . The altitudinal gradient in the occurrence of rodents and pikas is roughly correlated with these distribution data, despite that pikas were also frequent in the LPC diet at mid-elevations (1,500-2,500 m) in our analysis, suggesting a versatile generalist feeding tactic that may be influenced more by prey availability than a specific food preference. The above interpretation is based on the assumption that prey remains in faeces reflect predation occurring in the same altitudinal zone. However, LPC can move across large distances in short periods of time 21 , potentially compromising the association between faecal sampling and predation location. The impacts of altitude on LPC prey composition bear further investigation. Additionally, we detected several ungulates in the LPC diet. Thus, this small cat may possibly catch young wild boars and Chinese gorals, which are smaller-bodied ungulates, whereas the occurrence of takins in the diet may most likely be explained as scavenging behaviour. Our analyses also show that LPC consume a wide variety of birds, from small passerines to large ground-dwelling pheasants, demonstrating the predator's versatile hunting skills. We identified 20 vertebrate taxa in the AGC diet, representing the most comprehensive dietary analysis reported for this species. Similar to the LPC, AGC also frequently consumed rodents and pikas, but the proportions of larger prey, such as pheasants and ungulates, were higher than those in the LPC diet, possibly reflecting their larger body size compared to that of LPC. Due to the small sample size of AGC, the quantitative dietary analysis results regarding prey diversity remain tentative. AGC appeared to occupy only a portion of the LPC's range in the study area, as AGC faeces were found only within a relatively narrow altitudinal zone (2,300-2,700 m), whereas LPC faeces were collected throughout the reserve across all altitudes from 1,250 to 3,200 m. Additionally, the AGC population density is likely to be only a fraction of that of LPC in the study area 18 , which is in accordance with the different numbers of faecal samples collected for the two species. Further investigation of the AGC diet is required to better understand this species' food web hierarchy and niche interactions with other sympatric carnivores.

Molecular dietary analysis.
Most previous studies of the LPC diet were based on microscopic examination of food remains in faeces where generally < 10 taxa were identified, often only to the family or higher taxonomic levels 20,21,23 , demonstrating the superior resolution power of the faecal DNA-based method. The number of vertebrate prey categories identified in the present study (37 taxa) considerably exceeds not only reports based on a traditional cloning-Sanger sequencing approach for the LPC diet in the same area (16 taxa) 35 and in South Korea (13 taxa) 25 , but also a study based on NGS analysis of LPC faeces collected in northern Pakistan (18 taxa) 50 . Technical advances represented by the NGS approach offer remarkably higher throughput and economy compared to traditional morphological analysis and sequencing methods, demonstrating the high detection sensitivity and assignment power of molecular-metabarcoding/NGS methods. The difference in the number of food taxa between the two NGS-based studies can be explained by the higher species diversity of potential prey in the However, species resolution and accuracy using molecular identification are critically dependent on the completeness and quality of the reference databases, and limited coverage of the potential prey can strongly impede species-level identification. We found that the pika sequences in LPC faeces matched the Ochotona spp. sequences with 94-99% identity, suggesting that there might be more than one pika species in the diet.  14 , only the 12 S rRNA sequence of O. cuzoniae is available in the GenBank database. Therefore, we only assigned this taxon to the genus. Pikas were the most frequently occurring food item at the individual taxon level; yet when analysed at the ordinal level, they showed a lower proportion of occurrence in the diet than rodents due to the lack of species-level resolution within the genus Ochotona. One promise of the molecular dietary analysis is that taxonomic assignments can be refined as reference databases improve. The species resolution of molecular identification approaches is also restricted by the extent of sequence divergence among closely-related species and can show varying degrees of discriminating power for different taxa. The DNA barcode used in our study demonstrated high species resolution in most mammalian groups, yet low power differentiating avian species (particularly babblers and ducks). Limited species resolution can be particularly challenging in fine-scale dietary analysis when the study area has high species diversity, such as ours, where nearly 200 bird species have been recorded 36 . Thus, a combination of DNA barcodes targeting specific food groups has been recommended to improve species-level identification 54 .

Ecosystem function of mesopredators and conservation implications. Removal of apex predators
from an ecosystem can lead to loss of top-down control and subsequent outbreaks of mesopredator populations, an effect known as "mesopredator release 55 ", which may result in declines in prey populations and eventual "ecological meltdown 56 ". Despite overall extirpation of former apex predators in many mountain forests in Southwest China, neither explosive growth of mesopredator populations, nor depletion of the main prey groups have been observed 18,36 . Habitat loss and fragmentation, hunting by humans, intra-guild competition and limited prey resources may account for the lack of apparent increase in mesopredators in this ecosystem. Ecological processes in such a complex ecosystem can be multivariate and more intricate, and the dynamics and consequences following natural and human disturbances may not fully agree with predications often deduced from analyses of comparatively simple ecosystems 56 . Therefore, a closer examination of the many environmental, biological (intraspecific and interspecific) and anthropogenic factors is required in a particular ecosystem to gain a better understanding of the population dynamics and to predict the patterns of ecological response.
Our results indicate that rodents and pikas constitute the majority of the felids' diets (particularly the LPC) in the forest ecosystem of Southwest China. These small mammals are generally considered pests, as they degrade forests and grasslands and have been targets of control by poisoning and execution in many regions of China 57,58 . However, small mammals assist with seed dispersal and improve plant diversity 59 , provide important prey for many carnivores and predatory birds 60,61 and may serve as keystone species to enhance ecosystem functioning 62,63 . Pest control by poisoning has been demonstrated to be ineffective and often incurs undesirable consequences on the environment and untargeted species 64 . In contrast, conserving predator communities is an ecologically safe, economic and sustainable strategy to control small mammalian populations and maintain ecosystem equilibrium, stability and resilience 65 .
In conclusion, our results reveal the highly diverse dietary habits and versatile foraging skills of two felids. The other mammalian carnivores (e.g., Asiatic black bears, hog badgers, Chinese ferret badgers, yellow-throated martens, and masked palm civets) in this ecosystem are largely omnivorous, but LPC and AGC are highly carnivorous and may play pivotal roles as predators 66 . Large proportions of the LPC and AGC diets were composed of pikas and rodents, suggesting that these predators play an important role controlling their populations in the mountain forests of Southwest China. Therefore, caution should be taken in pest control as it may pose serious threats to predators and ecosystem sustainability. The faecal DNA and NGS-based approach enabled a very large-scale dietary analysis and fine species resolution. The remarkable taxonomic discriminating power offered by the molecular approach was particularly powerful to disentangle the food web structure in a highly complex ecosystem such as our study area, where closely related and morphologically similar prey species coexist.