A morphometric modeling approach to distinguishing among bobcat, coyote and gray fox scats

Coyotes Canis latrans, bobcats Lynx rufus and gray foxes Urocyon cinereoargenteus are all common mammalian mesopredators in coastal California and are found sympatrically in much of North America. Scats produced by these three animals are quite similar, but have historically been differentiated largely by morphology. I tested the efficacy of morphological classification of scat to species by building predictive models for species identification with a set of well-described, DNA-verified scats. I compiled a database of morphological, biogeochemical and contextual traits for a set of 122 DNA-verified bobcat, coyote and gray fox scats. I then took two different approaches to predictive modeling, using both discriminant function analysis and random forests to predict scats to species. I found significant differences among species in only three (diameter, mass and C:N ratio) of the 12 variables I considered. Linear discriminant analysis was only 71% predictive with the inclusion of a non-morphological variable in addition to morphological traits. Random forests similarly had only a 62% correct classification rate. Although scat morphology is not generally diagnostic to species for this set of mammalian mesopredators, these predictive morphometric models may still prove to be useful first-pass identification tools. The linear discriminant model in particular is able to identify scats with certain traits to species with a high degree of confidence, lending credence to the idea of ‘end member morphologies’ for scats produced by these different animals. I suggest that researchers take similar measurements to either use in the morphometric models presented here, or build similar models for their target species. These results also suggest that some previous studies using morphology-based scat identifications may have misrepresented or misinterpreted diets and space use by these sympatric mammals.

Feces are rich sources of information: they are ubiquitous, easily located in the field and their collection generally causes little disturbance. Scat is therefore used frequently in wildlife ecology to evaluate the presence/absence of elusive animals (Palomares et al. 2002), estimate animal abundances (Hurlbert 1971, Kohn et al. 1999, Prugh et al. 2008, characterize diets (Reynolds and Aebischer 1991, Symondson 2002, Deagle et al. 2005, and investigate animal health or disease ecology (Gompper et al. 2003, Liccioli et al. 2012. Historically, scats have been identified to species by morphology, but the morphological distinctions among scats produced by mammalian mesopredators can be difficult to discern and improperly identified specimens can confound and even invalidate research (Harrington et al. 2009). Despite some successes with field based scat identifications (Zuercher et al. 2003, Prugh andRitland 2005), which are more reliable when made in conjunction with additional natural sign (e.g. tracks), a number of studies have now documented the pitfalls inherent in relying on morphology alone (Bulinski and McArthur 2000, Davison et al. 2002, Reed et al. 2004, Harrison 2006. Even in the wake of these findings, researchers continue to rely on morphology (generally diameter) to some degree for species identification, which often results in the conservative exclusion of many collected scat samples (Neale and Sacks 2001a, b, Grigione et al. 2011, Lukasik and Alexander 2012, Dowd and Gese 2012 and the potential for a relatively high proportion of scats to be misassigned to species. Molecular scatology can be employed to more reliably identify scat to species (Foran et al. 1997, Kohn et al. 1999, Bidlack et al. 2007), but it can be a prohibitively expensive tool when large numbers of scats need to be identified and processed. If scat morphology is distinctive, then careful documentation of the morphology of mitochondrial (mt) DNA-verified scats will make it possible for researchers to take scat measurements and confidently assign species identifications.
The goal of this paper is to evaluate whether a morphometric approach to scat identification is sufficient for distinguishing among scats produced by three of the most common mammalian mesopredators in coastal California: coyotes Canis latrans, bobcats Lynx rufus and gray foxes Urocyon cineroargenteus. To that end, I compiled a database of morphological, biogeochemical and contextual traits for a set of scats DNA-verified to species and built predictive models for species identification using two different methods: discriminant function analysis and random forests, a composite tree-based modeling approach. Reed et al. (2004) also used discriminant analysis to try and differentiate between scats produced by coyotes and Mexican gray wolves Canis lupus baileyi based on scat diameter, mass and length measurements, and they found that these three variables were not particularly predictive for these two species. This approach has not yet been taken for the suite of species and scat variables considered here. If morphology truly differs among coyote, bobcat and gray fox scats, then scats produced by these different carnivores should be statistically separable into groups according to differences in their morphologies (and biogeochemistry) and the predictive models should have low misclassification rates.  Fig. A1). Año Nuevo State Park is located about 30 km north of Santa Cruz, California. A number of different habitats occur in the park, including native dunes, coastal terrace prairie and mixed evergreen forest. Younger Lagoon Natural Reserve is part of the Univ. of California (UC) Natural Reserve System. It protects a remnant y-shaped lagoon on the north side of the town of Santa Cruz, California. Dense coastal shrub, willow thickets and coastal prairie occur within the reserve along with salt and freshwater marsh. Moore Creek Preserve is situated directly inland from Younger Lagoon on the opposite side of Highway One. A city park since 1998, the Preserve primarily protects managed coastal prairie habitat, which is periodically grazed. Although dogs are not allowed in any of these parks and preserves, when encountered, dog scat was easily avoidable because of its uniform texture and color.

Data collection
I collected mesopredator scats quarterly from 2011-2013 to capture changes in scat morphology due to seasonal dietary and weather differences. After first clearing transects of all scats, I returned one week later to collect the scats deposited during the intervening week. At the time of collection, I recorded scat locations with a GPS and placed scats in individually marked ziploc-bags with a desiccant to reduce moisture and enhance DNA preservation. I measured and recorded the following in the field: scat diameter, scat length, number of scat pieces, length of taper and degree of taper (Table 1). I also classified scats as segmented, ropey, both, or neither. In the lab, I recorded scat dry weight following both freeze-drying and oven-drying. Though odor is also often cited as a possible indicator of species (Stokes and Stokes 1986, Bang 2001, Halfpenny 2008), I did not include scat odor as a variable in this study because it is not easily quantifiable. I sent a subset of these well-characterized scat samples (n  122) to Wildlife Genetics International for mtDNAbased species identification. I collected scats under California Fish and Game permit SC-11995 and with the approval of the UC Santa Cruz IACUC (protocol Kochp1105).
In addition to recording the morphological traits of the scats, I also considered three non-morphological variables: the location of the scats on the trail or road, presence of other sign (e.g. tracks or scrape -a scratch mark left on the ground), and the carbon to nitrogen (C:N) ratio of the scat (Table 1). There are some known differences among the sign characteristics of these three species. It is reported that both bobcats and coyotes scratch the ground and leave scent marks to help maintain their territories (Bekoff 1977, Larivière andWalton 1997). Coyotes have scent glands between their toes and may scratch with one hind foot adjacent to feces or urine Ruff 1997, Elbroch et al. 2012). Bobcats will sometimes create more controlled repeated scrapes with both hind feet and deposit urine or scat at one end (Elbroch et al. 2012). Scrapes are not commonly associated with gray foxes, though they, like coyotes and bobcats, will often deposit scats in groups (Fritzell and Haroldson 1982). C:N ratios vary across ecosystems and through food webs, reflecting underlying organismal allocations to major molecules and chemical structures; terrestrial vascular plants tend to have high C:N ratios (Meyers 1994, Prahl et al. 1994 while animals tend to be much more nutrient rich and therefore have much lower C:N ratios (Sterner and Elser 2002). The C:N ratio of scat, then, should serve as a proxy for an animal's degree of carnivory; animals consuming a largely plant-based diet will produce scats with high C:N ratios and those consuming other animals will produce scats with low C:N ratios. Because they are obligate carnivores, bobcats should produce scats with relatively low C:N ratios. Coyotes and gray foxes both incorporate some amount of plant material, mainly fruit, into their diets, therefore their scats should generally have higher C:N ratios. I obtained C:N ratios of the samples as a by-product of stable carbon and nitrogen isotope analyses conducted as part of a different yet unpublished project. To prepare the samples for analysis, I extracted matrix samples from scats by gently breaking apart dried scats over a fine mesh sieve. Matrix material passed through the sieve while other scat components (e.g. fur, feathers or bones) were captured above. I cleaned the powdery matrix material for stable isotope analysis by placing it into filter paper cones and rinsing it first with MilliQ water and then with 0.1N HCl to remove possible CaCO 3 contaminants. After the sample was fully dry and homogenized, I weighed approximately 5 mg of scat matrix material into Sn boats. The samples were then combusted via Dumas combustion using an elemental analyser and analysed on a continuous flow isotope ratio mass spectrometer at the UCSC Stable Isotope Laboratory. I calculated the average analytical precision for the scat data as the SD of the C:N ratio of 41 replicates of an internationally calibrated in-house standard (PUGel); precision was 0.2.

Fecal genotyping
I conducted fecal genotyping in collaboration with Wildlife Genetics International (WGI). Before freeze-and oven-drying, I swabbed the scats with Q-tips, which were NA is there a scrape mark near the scat? 1  yes, 0  no C:N Ratio unitless ratio of carbon to nitrogen atoms in the scat, which is a proxy for the degree of carnivory of the animal then stored dried in unwaxed coin envelopes. DNA was extracted by clipping a small (∼3  3 mm) piece of each swab and processing the clippings as tissue samples using QIAGEN DNeasy Blood and Tissue Kits. The species test is a sequence-based analysis of the mitochondrial 16S rRNA gene (Johnson and O'Brien 1997). Two variants of this analysis were employed using either primers that amplify across all mammals or primers designed to preferentially amplify Carnivora sequences. WGI then compared the results to reference data from over 125 species of mammals.

Statistical analyses and predictive model construction
I performed one-way analysis of variance (ANOVA) to identify possible significant differences in the means of the traits of scats produced by the three different species. I first tested the assumption of normality with the Shapiro-Wilk test and that of homogeneity of variance with the Bartlett test. Two of the variables I considered (mass and C:N ratio) did not meet the basic assumptions of ANOVA, which I addressed with a log transformation. When I observed statistically significant differences among groups with ANOVA, I used the post hoc Tukey-test to determine which of the three species contributed to the differences. Means are reported  one standard deviation (SD) and I tested significance at the p  0.05 level. Some scats (a total of eight) lacked diameter, length, taper length, and/or taper index measurements not because the measurements were neglected, but because the scats had irregular morphologies, rendering those measurements irrelevant. These samples were coded as 'flat' and I excluded them from the calculation of trait summary statistics, but not from the morphology models (Supplementary material Appendix 2 Fig. A2). Irregularity was more common to gray foxes (of the eight irregular scats, one failed identification, one had mixed DNA, one was from a coyote, and the remaining five were from gray foxes), therefore I thought it important to include these scats in the RF models. All statistics were performed in R ( www.r-project.org ).
In an effort to build the most predictive model for identifying scats to species, I compared the results from two different approaches, multiple discriminant function analysis and random forests, each of which has different strengths and limitations. Discriminant function analysis (DFA) uses linear combinations of numerical predictor variables to optimally separate a categorical response variable -in this case 'species' -and subsequently predict group membership for samples of unknown origin. To achieve this, it identifies gradients of variation among groups such that differences between groups are maximized, while within group variation is minimized (McGarigal et al. 2000). The assumptions for DFA include multivariate normality, equality of variancecovariance matrices across groups, and independent observations (McGarigal et al. 2000). Generally, DFA is robust to violations of these assumptions, at least when sample sizes are large (Williams andTitus 1988, McGarigal et al. 2000). When the covariance matrices among classes are different, then quadratic discriminant function analysis (QDA) may be a more appropriate alternative; the quadratic discriminant function is very similar to linear discriminant function analysis (LDA), but accommodates a class-specific covariance structure (Kuhn 2013). Williams and Titus (1988) proposed a general rule of thumb for sample sizes, specifying that each group should have N  3P, where N is the number of occurrences in a group and P is the number of discriminating variables. My sample sizes are relatively small, so I limited the number of predictor variables in the final models to five or fewer. I standardized the variables and examined the covariance matrices for each group to check for their equality and inspected correlation matrices to check for possible collinearity among variables. The coyote, bobcat and gray fox covariance matrices are not equivalent, so I performed QDA in addition to LDA. I separated my data into distinct training and test sets (70% and 30%, respectively) and evaluated classification error using the test set. I also assessed classification error on the full data set using a jackknife classification method, which is a leave-one-out cross-validation procedure. I took a stepwise approach to DFA, beginning with an overfitted model that included all possible variables, and removed the least predictive variables one at a time. I identified the least predictive variables by examining the coefficients of the linear discriminant function, which communicate information about the relative importance of the predictor variables when the data have been centered and scaled (Kuhn 2013). I performed LDA and QDA first with the morphological variables alone and again with the inclusion of C:N ratios, most predictive morphology only LDA model included four variables: number of scat pieces, diameter, taper length and log mass. Diameter and log mass contributed most strongly to both the first and second linear discriminant (Table 2). There was a statistically significant difference among the centroids on the first linear discriminant (F 2,57  47.99, p  0.001); the post hoc Tukey-test indicated that all three species differ significantly from one another. The statistically significant difference among centroids on the second discriminant (F 2,57  16.79, p  0.001) was, however, driven solely by the coyote scat centroid. Only 60% of the scats in the training set were correctly identified to species and performance in the test set increased to 67% correctly classified. For the entire dataset, 62% of the jackknifed predictions were correct. There also were species-specific differences in model performance; bobcat scats were consistently classified correctly more often than scats from the other two species.
The inclusion of scat C:N ratio improved predictions slightly. The best performance was achieved by a model that included number of scat pieces, diameter, taper length, log mass and log C:N ratio. Log C:N ratio and log mass contribute most strongly to the first linear discriminant while diameter and log C:N ratio contribute most strongly to the second (Table 2). Statistically significant differences among centroids are driven by both coyote and gray fox scats along the first linear discriminant (F 2,57  47.99, p  0.001) and by coyote scats alone along the second (F 2,57  16.79, p  0.001). Seventy four percent of the scats in the training set were correctly classified, followed by 75% in the test set and 71% in the jackknifed full dataset, all of which are comparable to the morphology alone model, though slightly improved. Coyote and gray fox scats were identified correctly more often than bobcat scats in both the test set and jackknifed full dataset. I held the composition of the training and test sets constant across models to facilitate comparison. Because the data set is small, the exact composition of the training set does impact model performance to a certain degree, but overall predictivity in the full jackknifed dataset never exceeded 75%.
Predictions did not improve significantly with the use of the QDA models. The morphology only model (length, diameter, taper length and log mass) correctly identified 68% of the scats in the training set to species. Performance in the test set increased to 75% correctly classified. With the addition of log C:N ratio as a predictor variable, 75% of the scats in the training set were correctly identified to species and just 71% in the test set. Notably, both QDA models correctly identified coyote scats with greater frequency than the LDA models.

Random forests
Results from the RF models were similar to the DFA model results. The out-of-bag estimate of correct classification was 57% for the morphology only model, in which mass and diameter were by far the most important variables (Fig. 3). When all of the variables were included, the out-of-bag estimate of correct classification rate improved only slightly to 62%. Mass remained the most important variable followed by C:N ratio and diameter (Fig. 3). Other non-morphological variables, including location and scrape, gained greater the only non-categorical non-morphological variable. I performed DFA in R using the 'MASS' package (Venables and Ripley 2002).
The random forest (RF) approach to classification and regression is a non-parametric method developed in machine learning (Breiman 2001). It fits an ensemble of hundreds or thousands of classification trees to a dataset using recursive partitioning, a method that is particularly well suited to small datasets with many explanatory variables (Strobl et al. 2009b). There are no distributional assumptions for either the predictor or response variables in RF and it is capable of handling missing and categorical data, making it potentially better suited for ecological datasets (Cutler et al. 2007, Strobl et al. 2009b. Because I have different types of predictor variables, I used cforest in the 'party' package in R (Hothorn et al. 2006a, Strobl et al. 2007 with the default option controls  cforest_unbiased() and I used varimp() to evaluate variable importance (Strobl et al. 2007(Strobl et al. , 2009a. I evaluated model fit with the built in out-of-bag classification estimation. To aid in visualization, I also created a single classification tree using the ctree() command in the 'party' package (Hothorn et al. 2006b) using a stop criterion based on the univariate p-values.

Scat characteristics
Of the 122 scats I submitted to Wildlife Genetics for identification to species, nine failed identification and two contained mixed DNA. The remaining scats include 57 bobcat, 28 coyote, 25 gray fox and one spotted skunk. I excluded the failed, mixed, and non-target scats from all subsequent analyses, leaving a dataset of 110 positively identified bobcat, coyote and gray fox scats. Of these, 19 are missing values for one or more variable and were therefore excluded from the models and calculations of summary statistics for the missing variables, leaving a total of 91 scats for the LDA and QDA models, but still 110 for the RF models.
Of the morphologic traits I considered, only scat diameter and mass were significantly different among groups (F 2,101  14.59, p  0.001; F 2,106  16.98, p  0.001; Fig. 1). Results from the post hoc Tukey-tests suggested that coyote scats (16.45  6.4 mm) have different diameters from both gray fox scats (11.5  4.1 mm) and bobcat scats (13.3  5.6 mm), but gray fox and bobcat scats did not have significantly different diameters from one another. Gray fox scats also had significantly different masses (6.6  2.8 g) from both coyote scats (17.9  13.9 g) and bobcat scats (12.4  6.6 g), though bobcat and coyote scats had indistinguishable masses. Of the three non-morphologic variables I considered, scat C:N ratio was also significantly different (F 2,105  29.2, p  0.001) and the post hoc Tukey-test revealed that all three species have statistically significant scat C:N ratios (coyote  8.6  2.3, bobcat  6.7  1.1, gray fox  12.1  5.6).

Discriminant function analyses
Morphological traits could not reliably distinguish scats produced by coyotes, bobcats, and gray foxes (Fig. 2). The  importance than many of the variables considered in the morphology only model. Of the misclassified scats, just one was bobcat, 23 were coyote and 11 were gray fox, suggesting that the RF model is better able to identify bobcat scats than either coyote or gray fox. The single classification tree identified 5 significant splits in the data (Fig. 4): one in the C:N ratio, two in mass, and one each in diameter and taper length.

Discussion
Both the RF and DFA model results suggest that scats produced by coyotes, bobcats, and gray foxes can neither be reliably distinguished by morphology alone nor with the inclusion of some additional non-morphological variables. Neither modeling method achieved greater than a 75% accuracy, regardless of the variables considered. This is not surprising, given that I observed very few significant differences in scat traits among species. Only diameter, mass, and C:N ratio were statistically separable by species, and of these, only scat C:N ratios were distinct for all three species. In general, there is substantial morphological overlap among these species, an observation that has been made previously for scat diameters (Weaver and Fritts 1979, Danner and Dodd 1982, Farrell et al. 2000, Reed et al. 2004). Despite the substantial morphological overlap, there are still some morphological differences among mammalian mesopredator scats. The second LDA model identified two gray fox scats correctly with  90% confidence, and these scats shared similarly high C:N ratios ( 19.1), were not segmented, and weighed very little ( 3.4 g). The second QDA model performed even better, identifying four gray fox scats correctly with  90% confidence and the same general traits hold. The average gray fox is smaller and more omnivorous than either the average coyote or bobcat, which explains why these traits (low mass and high C:N ratio) would together be highly predictive of gray fox. Misclassified gray fox scats (LDA: n  8, QDA: n  7), on the other hand, all had lower C:N ratios (LDA: 8.59  4.0; QDA: 7.0  0.8) and tended to be segmented or to have larger diameters (LDA: 13.5  4.4 mm; QDA: 14.2  4.8 mm) than gray fox scats on average. Similarly, the second LDA model correctly identified three coyote scats with  87% confidence, all of which had large diameters (25.0  0.91 mm), long tapers (56.8  30 mm), and intermediate C:N ratios (8.93  2.2).
These results suggest that when certain scat characteristics are observed together, identification to species can be made confidently, and there is credence to the idea that there are end member morphologies -i.e. a suite of morphological traits that are characteristic of scats produced by one species -for scats produced by these mesopredators. Field guides typically describe coyote scats as twisted and tapered ropes, while bobcat scats are often described as tubular, smooth, and segmented with blunt ends (Halfpenny 2008, Elbroch et al. 2012. I found that 84% of the bobcat scats in my full dataset are segmented. Because bobcats are obligate carnivores, bobcat scats also tend to have lower C:N ratios, thus, a segmented scat with a low C:N ratio may therefore be attributable to a bobcat with a higher degree of confidence. I also only observed scrapes next to bobcat scats, suggesting that the presence of a scrape is more likely an indicator for bobcats. Both coyotes and bobcats are known to leave scrapes next to urine and/or feces (Larivière andWalton 1997, Bekoff 1977), though coyotes tend to do so only at the edges of their territories (Gese and Ruff 1997). Out of the 110 positively identified scats in this dataset, five were found adjacent to a scrape and all five of these scats were identified as bobcat. Scats exhibiting mixtures of typically 'bobcat' or 'coyote' traits are difficult for the DFA models to assign to species with a high degree of confidence (Fig. 5). For both the second LDA and second QDA models, the majority of species misidentifications occurred when the posterior probabilities rested between 50-60%, whereas only a handful of scats were assigned to the incorrect species with a high degree of confidence.
Because strong false positives are rare ( 5%), either the RF or the DFA models may still be useful tools for scat identification when used in conjunction with alternate identification methods. I recommend that researchers take scat measurements (all of those considered here or just those that proved to be important in the models) and enter them into a spreadsheet for use in R. Both a detailed protocol and the R code necessary to run the models can be found in the Supplementary Material (Supplementary material Appendix 3). Given that the inclusion of C:N ratios marginally improved model performance, researchers may choose to include this variable as well. These analyses typically run ∼$6-$10 per sample, which is still considerably less expensive than the ∼$25 per sample lab and labor costs required for mtDNA analyses for this project. Once a researcher has run their own scat data through the model, they can export the modelpredicted scat classifications as well as the posterior probabilities for each species assignment. Depending on the research question and target species, the researcher can pick a posterior probability cut-off of a certain value (e.g. 85 or 90%)  and scats with predictions to species that fall below that value would then require verification by another method. Scats assigned to species with posterior probabilities above that threshold can simply be reported as that species with the assigned probability. Researchers can be fairly confident in the designations above their cut-off value because strong false positives are so rare. Choice of model will depend on what the species of interest is; the RF model presented here was better able to accurately identify bobcat scats, while the QDA models performed better for coyote scats. If attempting to distinguish among scats produced by a different set of mammalian mesopredator species, researchers could take a similar approach by taking measurements, having a subset of samples DNA-verified to species, and building predictive morphometric models for the species in question. For example, Reed et al. (2004) used linear discriminant analysis to distinguish between coyote and Mexican gray wolf scats with results comparable to those reported here; they found that the model including diameter and mass provided the most accurate identification for both species (Mexican gray wolf, 65%; coyote, 86%). And these results seem to align well with conventional wisdom regarding coyote and wolf scat -that is that coyote and wolf scat are morphologically quite similar but different in size (Weaver and Fritts 1979, Halfpenny 2008, Elbroch et al. 2012. Differences among scats produced by other species may, however, be more pronounced, so it may be that this morphometric modeling approach would prove more valuable for a different species assemblage. The real strength of this tool is that it allows researchers to attach a quantitatively derived probability value to their species designations for scats.

Conclusions
These results demonstrate that bobcat, coyote, and gray fox scats are morphologically very similar and that contextual information such as location on trail or the presence or absence of a scrape only minimally increases the probability of correct identification. Additional information, such as proximity of fresh tracks, could improve identifications, but this information is often lacking, particularly when many scats are collected at once. The LDA and QDA models and classification tree presented here provide researchers with a first-pass tool to predict unknown scats to species with a pre-determined level of confidence and to then seek verification by other methods for scats predicted to species below the chosen confidence threshold. Either mtDNA analyses or scat-detecting dogs (Smith et al. 2003) could be used to provide this check. Alternatively, some researchers have also previously used the presence of a small number of guard hairs in scats to aid in successful identification to species (Miotto et al. 2007). Using the model to select scats needing further verification allows for at least some reduction in cost. Given the significant overlap in their scat morphology, it becomes critically important for researchers and wildlife managers investigating space use or possible dietary partitioning by these animals to assess morphology-based identifications with alternate methods. Some previous research on these animals that relied on scat morphology alone for species identification may benefit from a second glance. Even with the conservative exclusion of morphologically ambiguous scats, if morphological cut-offs were used to differentiate among species, it is likely that a greater proportion of scats were misidentified than previously thought.