Otterly delicious: Spatiotemporal variation in the diet of a recovering population of Eurasian otters (Lutra lutra) revealed through DNA metabarcoding and morphological analysis of prey remains

Eurasian otters are apex predators of freshwater ecosystems and a recovering species across much of their European range; investigating the dietary variation of this predator over time and space therefore provides opportunities to identify changes in freshwater trophic interactions and factors influencing the conservation of otter populations. Here we sampled faeces from 300 dead otters across England and Wales between 2007 and 2016, conducting both morphological analysis of prey remains and dietary DNA metabarcoding. Comparison of these methods showed that greater taxonomic resolution and breadth could be achieved using DNA metabarcoding but combining data from both methodologies gave the most comprehensive dietary description. All otter demographics exploited a broad range of taxa and variation likely reflected changes in prey distributions and availability across the landscape. This study provides novel insights into the trophic generalism and adaptability of otters across Britain, which is likely to have aided their recent population recovery, and may increase their resilience to future environmental changes.Funding provided by: Knowledge Economy Skills Scholarships*Crossref Funder Registry ID: Award Number:Sample and data collection Samples and associated metadata were acquired from 300 otters collected between 2007 and 2016, obtained from the Cardiff University Otter Project, a national monitoring programme for dead otters sampled from across Great Britain (https://www.cardiff.ac.uk/otter-project). Most otters collected were killed by road traffic accidents, with a minority dying through drowning, being shot, starvation, or disease. Information on date (year and month) and location (as grid reference) of carcass collection were recorded at the site of collection. Grid references were used to plot data for spatial analysis. Detailed post-mortems were performed for each carcass during which biotic data were obtained (e.g., sex and size of individual). Faecal samples were collected from the rectum during post-mortem examination, wrapped in foil and stored at -20 °C. Following post-mortems, scaled mass index (SMI) was calculated for each individual otter using the following equation (Peig and Green 2009; Peig and Green 2010): SMI = Mi [ L0 / Li ] bSMA Mi is the body mass and Li is the length measurement of individual i, L0 is the mean length measurement for the entire study population and bSMA is the scaling exponent. Length was measured from nose to tail-tip to the nearest 5 mm. Mean length and the scaling exponent were both calculated from all otter data available as of January 2017 (n = 2477). The scaling exponent is the slope from the standard major axis regression of log-transformed values of mass against length.Otters were also classified into size categories based on their total length (nose to tail tip) using the 'bins' function in R (OneR v2.2 package; von Jouanne-Diedrich 2017), which applies a clustering method using Jenks natural breaks optimisation. Male and female otters were clustered separately into small (males 1131 mm, females > 1031 mm). Spatial data Spatial data describing proximity to the coast, urban land use, altitude, slope, and primary water habitat were collated using QGIS version 3.4.4 (QGIS Development Team 2019). Distance from the coast was calculated as the shortest distance (km) along a river from the location at which the otter was found to the low tide point of the mouth of the river (hereafter referred to as 'river distance'), using the package RivEX (Hornby 2020), because otters tend to travel along water courses rather than across land. As most otters were found as roadkill, and not all were adjacent to rivers, each otter was first assigned to the nearest river. Locations more than 1000 m from a river were checked, and if there was more than one river along which the otter might have travelled, then river distance was calculated for all rivers and a mean distance used. All otters within 1000 m of the coast were given a distance of zero if they were closer to the coastline than a river.Otter locations were mapped as points, and circular areas of 10 km diameter (hereafter referred to as 'buffers') mapped around each. Faecal samples typically reflect diet from the preceding 24-72 hours (in mammals; Deagle et al. 2005; Casper et al. 2007; Thalinger et al. 2016), during which time otters can travel up to 10 km (Chanin 2003), it was therefore deemed appropriate to use this distance to reflect the land used by otters within the sample timeframe. Buffers were used to calculate proportions of urban land-use (i.e., urban and suburban land use extracted from the 25 m resolution UK land cover map from 2007; Morton et al. 2011), mean altitude and mean slope (extracted from European Digital Elevation Model (EU-DEM) map; European Environment Agency 2011). We chose to focus on urban land-use as urbanisation may affect otter diet either through changes to prey assemblages or disturbance affecting an otter's ability to forage. Longitude, altitude, and slope were highly correlated, therefore longitude was used in further analyses as a representative for the three variables.Otters in England and Wales typically feed in freshwater river systems but will opportunistically feed in lakes or at the coast if these habitats are within range (Jędrzejewska et al. 2001; Clavero et al. 2004; Parry et al. 2011). Available prey differ between lakes, coasts, and river systems, as well as between different parts of the river network (e.g., tributary, main river channel). To assess whether water habitat type influenced dietary variation, we designated each otter to one of the following: transitional water (coastal and estuarine), lake, main river channel, or tributary (based on Water Framework Directive 2000/60/EC designations mapped using GIS shapefiles provided by Natural Resources Wales and Environment Agency). Otters within 2.5 km (half of a buffer's radius) from a lake or transitional water were assigned to that habitat, whilst those further away were assumed to be feeding primarily in the river network. The RivEX network map (Hornby 2020) was used to map all rivers, and individuals were further categorised according to whether their assumed habitat was primarily a main river or tributary. To do this, the total length of main river channels and tributaries was calculated within each 10 km buffer. The length of main channels was weighted 10 times greater to account for the greater cross-section of a main channel compared to tributaries (Benda et al. 2004) since waterways with greater areas are assumed to support more prey (Samarasin et al. 2014). The sum of weighted main river lengths and tributary lengths was calculated, and if more than 50 percent of each buffer was attributed to main river channel, the otter was assigned to the main river channel, otherwise it was assigned to tributary. Morphological analysis Each faecal sample was first thawed, homogenised by hand in a sterile container, and divided into subsamples; three samples weighing 200 mg each were collected for DNA analysis (one sample used for DNA extraction and the other two frozen as back-ups) and the remaining material was used for morphological analysis. Subsamples undergoing morphological analysis were then soaked in a solution of water and liquid biological detergent (water:detergent, 10:1) for 24 hours. Samples were passed through sieves with a 0.5 mm mesh and washed with water to ensure only hard parts remained which were air-dried for 24 hours. A record was made of any samples that did not contain any hard parts. Recognisable remains (bones, fish scales, feathers, fur) underwent microscopic identification using a range of keys (Libois and Hallet-Libois 1987; Coburn and Gaglione 1992; Prenda and Granado-Lorencio 1992; Prenda et al. 1997; Watt et al. 1997; Miranda and Escala 2002; Conroy et al. 2005; Tercerie et al. 2019; University of Nottingham 2020). Prey remains were identified to the finest possible taxonomic resolution and recorded as present within or absent from a sample. DNA metabarcoding analysis Faecal samples were processed for HTS, and subsequent bioinformatic analysis was conducted, as described in Drake et al. (2022). In summary, DNA was extracted from a subsample of faecal material and amplified using two metabarcoding primer pairs, designed to amplify regions of the 16S rRNA and cytochrome c oxidase subunit I (COI) genes, each primer having 10 base pair molecular identifier tags (MID tags) to facilitate post-bioinformatic sample identification.Two primer sets from different barcoding regions were selected to overcome biases associated with each region and broaden the range of taxa amplified: the 16S barcoding region selected for vertebrate DNA, and cytochrome c oxidase subunit I (COI) for invertebrate DNA. For 16S, the novel primer pair FN2199 (5'- yayaagacgagaagaccct -3') and R8B7 (5'- ttatccctrgggtarcthgg -3'; modified for this study from Deagle et al. 2009) were used, which targeted a 225-267 bp amplicon (including primers). For COI, the primer pair Mod_mCOIintF (5'- ggwacwggwtgaacwgtwtaycc -3'; modified for this study from Leray et al. 2013) and HCO-2198 (5'- taaacttcagggtgaccaaaaaatca -3'; Folmer et al. 1994) were used, which targeted a 365 bp amplicon (including primers). Primers underwent in silico testing using ecoPCR (Boyer et al. 2016) and were further tested in vitro. In silico and vitro tests showed primer sets amplified desired taxa. COI primers were also found to amplify a range of vertebrate taxa but did not cover the same range as the 16S primers, therefore justifying the use of both primer sets beyond the benefit of different primer pairs exhibiting different biases.Faecal samples were processed alongside extraction and PCR negative controls, repeat samples, and mock communities, which comprised standardised mixtures of DNA of marine species not previously detected in the diet of Eurasian otters, and some tag combinations were left unused to identify tag jumping. The resultant DNA libraries for each marker were sequenced on separate MiSeq V2 chips with 2 × 250 bp paired-end reads. Bioinformatic analyses were carried out using a custom pipeline, following which sequencing data underwent filtering steps to remove any remaining artefacts or contaminants in the data (Drake et al. 2022). Filtering involved removing taxa from each sample that contributed less than 0.5 % of a sample's total reads for 16S and 0.3 % for COI. Reads equal to or less than the maximum read count identified in unused MID-tag combinations or negative controls per taxon were also removed. This method was conservative but was selected to remove false positives which would otherwise overrepresent some prey groups present in some samples as contaminants (Drake et al. 2022).Reads were assigned to the finest possible taxonomic resolution and recorded as present within or absent from a sample, separately for 16S and COI. Reads assigned to non-food items remaining in the analysis were removed, these included taxa not assigned to the animal kingdom (e.g., fungi and bacteria, which were not considered pertinent to this study), those with poor taxonomic resolution (e.g., Eutheria, which includes all extant British mammals and thus was not useful for further analyses), reads from otters themselves (e.g. those assigned to Lutra lutra; Cuff, Kitson et al. 2022) and taxa with a maximum size < 3 mm (e.g. diatoms, assumed to be due to secondary or accidental predation). Following removal of non-food items, data from the two datasets were combined to give a more complete representation of the diet of otters through the complementarity of the separate taxonomic biases of the two primer pairs. If a taxon was present in either of the metabarcoding datasets, then that taxon was recorded as present in that sample. If a prey item was detected in a sample in both metabarcoding datasets, but at different levels of taxonomic resolution, only the presence with the finer taxonomic resolution was retained. Comparison of methods The frequency of occurrence for each prey item detected across the 300 otters screened was calculated for both morphological and metabarcoding datasets, allowing the two methods to be directly compared. Presences assigned to 'insect', 'beetle', 'mollusc' and 'snail' in morphological analysis were removed before comparing datasets; many identifications from these particular taxonomic groups were identified to a finer resolution through metabarcoding but removed as secondary predation or accidental consumption (Tercel et al. 2021), therefore these presences in the morphological analysis were also deemed likely to have occurred through secondary predation. Presences assigned to 'mammal' (identified from fur) in the morphological analysis were also removed before comparing datasets due to the uncertainty of fur coming from the otter grooming itself and metabarcoding identifying otter as the only mammal in these samples. Presence-absence matrices produced from each methodology were also combined in order to assess the overlap in data. Where both methods identified the same taxonomic group, the presence was assigned to the finest taxonomic resolution, whereas where there was ambiguity about whether the methods were detecting the same taxon or not, the presence was assigned to a coarser taxon. Combining datasets from each method revealed which data points were only detected by one method and which were detected by both (either at the same taxonomic level or at different resolutions). Binary matrices for prey detections were combined for the two data types, but each sample was represented separately for each method (i.e., not aggregated by sample). Statistical analysis The association between otter diet composition and biotic and abiotic drivers was explored using the combined data from morphological analysis and metabarcoding. Each taxon was assigned to a 'prey group' based on similarities in taxonomy, morphology, and ecological niche. A small number of prey identified to coarse taxonomic levels could not be assigned to a group, and were removed from subsequent analyses (prey presences removed: 'Salmo sp.', n = 5; 'Cyprinid', n = 2; 'Bird', n = 2). A prey group was recorded as present in an individual faecal sample if any one (or more) of the taxa assigned to that group were present. If a prey group occurred in less than three samples then the prey group was designated as rare and removed from subsequent analyses. Dietary composition was compared against biotic and abiotic drivers via multivariate generalized linear models in R [version 3.6.0] and R Studio [version 1.2.1335] (R Core Team 2019) using 'mvabund', and visualised using 'bipartitie' and 'boral' packages.The 'mvabund' package allows model-based analysis of multivariate data to test hypotheses regarding the effects of environmental variables on the composition of dietary data (Wang et al. 2012). Multivariate generalized linear models (MGLMs) are a robust method for detecting differences in communities with less abundant taxa and are less prone to misinterpretations due to mean-variance effects, compared to distance-based methods (Warton et al. 2012). The 'many.glm' function was used to create a MGLM using a binomial family and a 'cloglog' link function. The global models included the following fixed variables: sex, size of otter, SMI, year, season, distance from the coast (km), primary water habitat, percentage of urban land-use, latitude, and longitude. Interactions between sex and size of otter, distance from the coast and sex, distance from the coast and size, primary water habitat and sex, primary water habitat, and size, and between latitude and longitude were also included in the global model. Model assumptions were checked on the global model before conducting model selection via Akaike's Information Criterion (AIC) using the stepwise algorithm in the step function (Hastie and Pregibon 1992; Venables and Ripley 2002). The final model included the fixed variables longitude and distance from the coast. The significance of fixed variables on the overall diet and for specific prey groups was determined via likelihood ratio test using the 'anova.manyglm' function with Monte Carlo resampling and corrected univariate p-values for multiple testing.To complement the mvabund analysis, the boral package was used to plot significant variables. The boral package conducts Bayesian ordination and regression analysis on multivariate data (Hui 2016). Binomial models for boral analysis included the same fixed and response variables as in the final mvabund model. The number of latent variables was set as two. Model assumptions were checked and latent variable values extracted. Latent variables were plotted against significant fixed variables to visualise the individual samples and the indicator species that best described their position in a low-dimension ordination plot. Bipartite network plots were also created to visualise the structure and identity of otter trophic interactions against significant fixed variables using the plotweb function in the bipartite package (Dormann et al. 2008). Data generated by both metabarcoding and hard-parts analysis were visually compared using non-metric multidimensional scaling via the 'metaMDS' command with a Jaccard distance matrix and 999 tries in the 'vegan' package (Oksanen et al., 2016). Only samples for which both metabarcoding and hard-parts data were available were included in these plots. References Benda, L., Andras, K., Miller, D. and Bigelow, P. (2004). Confluence effects in rivers: Interactions of basin scale, network geometry, and disturbance regimes. Water Resources Research, 40(5), pp. 1–15. doi: 10.1029/2003WR002583. Boyer, F., Mercier, C., Bonin, A., Le Bras, Y., Taberlet, P. and Coissac, E., 2016. obitools: A unix‐inspired software package for DNA metabarcoding. Molecular ecology resources, 16(1), pp.176-182. Casper, R.M., Jarman, S.N., Deagle, B.E., Gales, N.J. and Hindell, M.A. (2007). Detecting prey from DNA in predator scats: A comparison with morphological analysis, using Arctocephalus seals fed a known diet. Journal of Experimental Marine Biology and Ecology, 347(1–2), pp. 144–154. doi: 10.1016/j.jembe.2007.04.002. Chanin, P. (2003). Ecology of the European Otter. In: Conserving Natura 2000 Rivers Ecology Series No. 10. English Nature, Peterborough, p. 64. Clavero, M., Prenda, J. and Delibes, M. (2004). Influence of spatial heterogeneity on coastal otter (Lutra lutra) prey consumption. (August), pp. 551–561. Coburn, M.M. and Gaglione, J.I., 1992. A comparative study of percid scales (Teleostei: Perciformes). Copeia, pp.986-1001. Conroy, J.W.H., Watt, J., Webb, J.B. and Jones, A. (2005). A Guide to the Identification of Prey Remains in Otter Spraint. 3rd edi. London, United Kingdom: The Mammal Society. Cuff J.P., Kitson, J.J.N., Hemprich-Bennett, D., Tercel, M.P.T.G., Browett, S.S., Evans, D.M. (2022). The predator problem and PCR primers in molecular dietary analysis: swamped or silenced; depth or breadth? Molecular Ecology Resources, Early View doi: 10.1111/1755-0998.13705 Deagle, B.E., Tollit, D.J., Jarman, S.N., Hindell, M.A., Trites, A.W. and Gales, N.J. (2005). Molecular scatology as a tool to study diet: Analysis of prey DNA in scats from captive Steller sea lions. Molecular Ecology, 14(6), pp. 1831–1842. doi: 10.1111/j.1365-294X.2005.02531.x. Deagle, B.E., Kirkwood, R. and Jarman, S.N. (2009). Analysis of Australian fur seal diet by pyrosequencing prey DNA in faeces. Molecular Ecology, 18(9), pp. 2022–2038. doi: 10.1111/j.1365-294X.2009.04158.x. Dormann, C., Gruber, B. and Fründ, J. (2008). Introducing the bipartite package: analysing ecological networks. Interaction, 1(0.2413793). Drake, L.E., Cuff, J.P., Young, R.E., Marchbank, A., Chadwick, E.A. and Symondson, W.O.C, 2022. An assessment of minimum sequence copy thresholds for identifying and reducing the prevalence of artefacts in dietary metabarcoding data. Methods in Ecology and Evolution, 13(3), pp.694-710. European Environment Agency (2011). European Digital Elevation Model (EU-DEM), version 1.1. Copernicus Programme, (Available at: http://land.copernicus.eu/pan-european/satellite-derived-products/eu-dem/eu-dem-v1.1/view) Folmer, O., Black, M., Hoeh, W., Lutz, R. and Vrijenhoek, R. (1994). DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Molecular Marine Biology and Biotechnology, 3(5), pp. 294–299. doi: 10.1371/journal.pone.0013102. Hastie, T.J. and Pregibon, D. (1992). Chapter 6: Generalized linear models. Statistical models in S: Wadsworth & Brooks/Cole Hornby, D. D. (2020). RivEX [Software]. Available from http://www.rivex.co.uk Hui, F.K.C. (2016). boral – Bayesian Ordination and Regression Analysis of Multivariate Abundance Data in r. Methods in Ecology and Evolution, 7(6), pp. 744–750. doi: 10.1111/2041-210X.12514. Jędrzejewska, B., Sidorovich, V.E., Pikulik, M.M. and Jędrzejewski, W. (2001). Feeding habits of the otter and the American mink in Białowieża Primeval Forest (Poland) compared to other Eurasian populations. Ecography, 24(2), pp. 165–180. Leray, M., Yang, J.Y., Meyer, C.P., Mills, S.C., Agudelo, N., Ranwez, V., Boehm, J.T. and Machida, R.J. (2013). A new versatile primer set targeting a short fragment of the mitochondrial COI region for metabarcoding metazoan diversity: Application for characterizing coral reef fish gut contents. Frontiers in Zoology, 10(1), pp. 1–14. doi: 10.1186/1742-9994-10-34. Libois, R. M., and Hallet-Libois, C. (1987). Eléments pour l'identificat