Evaluation of genetic markers for the metabarcoding of Australian marsupials from predator scats

Selection of candidate metabarcoding primers

Candidate marsupial metabarcoding primers were identified from the literature on the basis of amplicon length (<300 bp) and organism specificity (Supplementary Table S1). In total, we evaluated five sets of primers that have previously been shown to amplify mammalian DNA (Table S1). These included two sets of primers for each of the 12S rRNA and 16S rRNA genes, and one set of primers for the COI gene. For the 12S rRNA gene, the 12SV5 primers were developed for use with degraded DNA and tested on carnivore scats (Riaz et al. 2011) and have been used on Tasmanian devil scats to identify diet items (McLennan et al. 2022), whereas the primers L1085 and H1259 were designed on the basis of human sequences and have been tested across multiple mammalian taxa with high success (Kitano et al. 2007). Both the 16SMam primers (Taylor 1996) and the L2513 and H2714 primers (Kitano et al. 2007) for the 16S rRNA gene were developed for use with ancient or degraded DNA. Most published sets of COI primers amplify longer fragments, but the uni-minbar primers (Meusnier et al. 2008) were specifically designed to amplify shorter fragments of degraded DNA.

Evaluation of primer utility for amplification of mammal DNA

We used two analyses to determine the likely effectiveness of each candidate primer pair for reliably amplifying mammalian DNA from unknown samples. First, we used the R package PrimerTree (https://github.com/MVesuviusC/primerTree, Cannon et al. 2016) to evaluate the potential utility of each primer set to amplify DNA from, and distinguish among, taxa within the class Mammalia. We queried each candidate primer set against GenBank, using the BLAST function within PrimerTree (with the default settings), to identify all vertebrate sequences that could be amplified. In this context, sequences were considered to be amplifiable if they included complementary binding sites for both the forward and reverse primers for a given marker. Results were visualised by generating phylogenetic trees showing all amplifiable sequences for each marker at the level of class (noting that these trees do not necessarily represent true species relationships owing to the very short nature of the amplified region). Branch lengths on these trees relate directly to the level of resolution provided by a given marker for particular taxa; short branch lengths indicate that the marker is likely to have very low resolution, whereas long branch lengths indicate a high capacity to distinguish among species (Cannon et al. 2016). We estimated the likelihood that each primer set would bind to a high proportion of sequences within the target group, and whether they could also be expected to bind to non-target groups.

In our second analysis, we evaluated each of the five primer sets at the infraclass level, to determine their suitability for application to different groups of mammals, especially marsupials. For each of the three candidate genes, all available marsupial DNA sequences were downloaded from GenBank (12S, 16S and COI) and BOLD v3.0 (COI only) (accessed February 2017) and sorted into operational taxonomic units (OTUs) by using the R package PrimerMiner-0.12 (https://github.com/VascoElbrecht/PrimerMiner, Elbrecht and Leese 2016). OTUs were aligned in Geneious v10.0.5 by using a MAFFT alignment (Katoh et al. 2002) with the default settings. Primers were evaluated against OTU alignments following the recommended PrimerMiner protocol. PrimerMiner evaluates single primers or primer pairs, considering the adjacency, position and type of each mismatch between primers and the alignment. Each primer pair was tested at the following different thresholds: 50, 100, 150 and 200. A higher threshold allows for more base-pair mismatches in the primer binding region. The higher the required threshold, the less suitable the primer is likely to be for a given OTU. Therefore, primer pairs with lower thresholds assigned for an OTU would be preferred for further analysis and use in a metabarcoding study.

Construction of locally appropriate reference-DNA databases for the Tasmanian Diprotodontia

On the basis of our initial evaluation (see Results), we selected two markers, 12SV5 and 16SMam, for further consideration, because these showed the greatest potential for broad amplification of mammal DNA, including marsupials. To allow evaluation of the diagnostic ability of these markers in the Tasmanian Diprotodontia, a group that was depauperate in genomic resources, we generated custom reference DNA-sequence databases for the 12S rRNA and 16S rRNA genes. We sequenced representatives of most extant species of the order Diprotodontia that currently occur in Tasmania. For each target species, samples were included from across the geographic range where possible, including outside Tasmania (Table S1). Tissue and skin samples were sourced from museums, researchers, and through the collection of road-killed animals. DNA was extracted from tissue samples following a salting-out method, as follows: for each sample, a piece of tissue approximately 2 mm² was incubated overnight at 55°C in 330 μL of tissue extraction buffer (40 mM Tris–HCl; 100 mM NaCl; 20 mM EDTA, pH 7.2; 0.6% SDS; 0.5 mg/mL Proteinase K). After incubation, 150 μL of 7.5 M ammonium acetate was added. Samples were chilled at −80°C for 25 min before centrifugation for 20 min at 19,000g and 4°C. Supernatant was transferred to a new tube and DNA was ethanol precipitated. DNA pellets were resuspended in 50 μL of 1 × TE (10 mM Tris-HCl; 1 mM EDTA) (MacDonald et al. 2011). In addition, following tissue lysis and incubation at 55°C, 3 μL of 4 ng/μL RNase A (Sigma-Aldrich) was added to each sample and incubated at 55°C for 30 min. DNA was extracted from skin samples using the Qiagen DNeasy® blood and tissue kit, following the manufacturer’s recommended protocol, with the addition of an extra 100 μL of ATL buffer (Qiagen) and 10 μL proteinase K during the lysis step, with incubation for 24 h.

An 800 bp region of the mitochondrial 12S rRNA gene was amplified using primers 12C (5′ AAAGCAAAKCACTGAAAATG 3′) and 12GG (5′ TRGGTGTARGCTRRRTGCTTT 3′) (Springer et al. 1995). PCRs of 25 μL contained 50 ng of DNA, 1 × MyTaq HS Red (Bioline), and 0.4 μM of each forward and reverse primer. PCR conditions consisted of denaturation at 95°C for 5 min, followed by 15 cycles of 95°C for 20 s, a touchdown of 65–50°C for 20 s and 72°C for 60 s, followed by 30 cycles of 95°C for 20 s, 50°C for 20 s, and 72°C for 60 s, followed by a final extension at 72°C for 4 min.

A 900 bp region of the mitochondrial 16S rRNA gene was amplified using primers Mammal_16SF (5′CGAGCCTGGTGATAGCTGGTT3′) (Modave 2017) and U16S-REV (5′ TCCGGTCTGAACTCAGATCAC 3′) (Yang et al. 2014). PCRs of 25 μL contained 50 ng of DNA, 1 × MyTaq HS Red (Bioline), and 0.4 μM of each forward and reverse primer. PCR conditions consisted of denaturation at 95°C for 5 min, followed by 15 cycles of 95°C for 20 s, a touchdown of 65–50°C for 20 s and 72°C for 60 s, followed by 30 cycles of 95°C for 20 s, 50°C for 20 s, and 72°C for 60 s, followed by a final extension at 72°C for 4 min.

PCR products were purified using Diffinity RapidTips (Sigma-Aldrich). Sequencing reactions of 20 μL included 1 μL of purified PCR product, 3.5 μL BigDye 3.1 (Applied Biosystems), 1 × BigDye sequencing buffer and 2 μM primer. The reaction conditions involved an initial denaturation at 94°C for 5 min, followed by 30 cycles of 96°C for 10 s, 50°C for 5 s, and 60°C for 4 min. Sequencing reactions were purified using a standard ethanol/EDTA precipitation method (Applied Biosystems 2009). Sequencing was performed on an AB 3730 × l DNA Analyser at the ACRF Biomolecular Resource Facility (Australian National University), following the manufacturer’s protocol. Forward and reverse raw sequences were aligned in Geneious v10.0.5 and subsequently checked and edited manually, before a consensus sequence was generated for each sample.

Our custom reference databases were supplemented with all publicly available 12S rRNA and 16S rRNA sequences from GenBank (accessed March 2017 for the target species; Table S1). The final DNA reference database for each gene was constructed in Geneious by aligning all sequences (consensus sequences generated in this study and sequences obtained from GenBank) using the MUSCLE option (Edgar 2004). The 12SV5 and 16SMam primers were then mapped to their respective custom reference database. For each marker in turn, the full gene alignment was trimmed to include only the sequence region between the primer binding sites, and these trimmed marker-specific alignments were used for subsequent analyses.

Genetic distance-based evaluation of the candidate metabarcoding markers

The R package SPIDER v1.1-2 (https://cran.r-project.org/web/packages/spider/index.html, Brown et al. 2012) was used to determine the diagnostic abilities of the 12SV5 and 16SMam markers, and the risks of false positive and false negative identifications associated with each. Singleton species (species for which only a single sequence was available) were removed from the two marker-specific reference databases, because intraspecific genetic differences could not be calculated for these species. Three analyses, the nearNeighbour, bestCloseMatch and threshID functions (included within the R package SPIDER), were then undertaken for each candidate marker (using the same code and parameters as described in (MacDonald and Sarre 2017)) to determine which sequences, if any, were likely to be ambiguously identified or misidentified. We evaluated genetic-distance thresholds of 1%, 3% and 5% for the bestCloseMatch and threshID analyses. Pairwise genetic distances were calculated using a raw or uncorrected model.

Phylogenetic evaluation of the candidate metabarcoding markers

We generated a maximum-likelihood phylogenetic tree for each marker, by using only the short gene regions amplified by each set of primers, to determine whether the 12SV5 and 16SMam markers have sufficient phylogenetic resolution to discriminate among the taxonomic groups of interest. The evolutionary history was inferred for the 92 bp 16SMam region by using the maximum-likelihood method based on the Tamura 3-parameter model (Tamura and Nei 1993). A discrete gamma distribution was used to model evolutionary rate differences among sites (+G, parameter = 0.7925) and was determined using the Akaike information criterion (AIC) calculations by using MEGA v7.0 (https://www.megasoftware.net/home, Kumar et al. 2016). The evolutionary history was inferred for the 100 bp 12SV5 fragment by using the maximum-likelihood method based on the Kimura 2-parameter model (Kimura 1980). A discrete gamma distribution was used to model evolutionary rate differences among sites (+G, parameter = 0.1723) and was determined using the AIC calculations by using MEGA v7.0 (https://www.megasoftware.net/home, Kumar et al. 2016). Phylogenetic trees were prepared using garli v2.01 (https://github.com/zwickl/garli, Zwickl 2006) with 1000 bootstraps. The trees were rooted using the genus Dasyurus as an outgroup, because the Dasyuromorphia represents a sister taxon to the Diprotodontia (Mitchell et al. 2014).

Field testing of candidate metabarcoding markers

We conducted a proof-of-concept analysis by using samples collected in one 3 km × 3 km grid square in north-eastern Tasmania as part of a broadscale predator-scat survey (Ramsey et al. 2018). Ten predator scats were collected on 1 day in the field in 2008, dried, and sent to the University of Canberra Trace DNA Laboratory for analysis. Laboratory and bioinformatic methods are outlined below.

Custom reference-DNA database for all known prey species in Tasmania

To enable interpretation of DNA metabarcoding data from predator scats, the two custom reference databases developed for Tasmanian Diprotodontia were augmented by the addition of sequences from other Tasmanian mammals. These included all extant species of Muridae, Dasyuridae, Felis and Leporidae that occur within Tasmania. We used Obitools v1.01 (https://pythonhosted.org/OBITools/introduction.html, Boyer et al. 2016) to generate metabarcoding reference databases including either 12SV5 or 16SMam sequences as appropriate, from all mammals represented in Genbank (Benson et al. 2013), to ensure that we represented as many potential predator and prey species likely to be encountered in Tasmania as possible. The final 12S rRNA custom DNA reference database constructed for prey analysis consisted of 321 sequences representing five orders, 10 families, 17 genera and 41 species. The 16S rRNA custom DNA reference database consisted of 204 sequences comprising five orders, 10 families, 17 genera and 35 species.

Scat survey

A series of predator-scat surveys were conducted across eastern Tasmania by the Department of Natural Resources and Environment, Tasmania, and the Invasive Animals Cooperative Research Centre (IACRC). Standardised methods were used, in which 3 km² survey units were searched for a specified duration of 10 person-h, following linear features (Ramsey et al. 2018). These surveys focused on three zones, namely, north-eastern Tasmania, north-central Tasmania and south-eastern Tasmania, in 2008, 2009 and 2010 respectively (Sarre et al. 2013). Subsequently, the north-eastern zone was re-surveyed, by using the same methods and survey unit sites, in autumn 2014. All mammal-predator scats encountered, including those from Tasmanian devils, spotted-tailed and eastern quolls, cats, dogs, and red foxes, were collected, stored in individual paper bags, and air dried in preparation for DNA analysis. Each scat was also swabbed in situ before collection during the 2014 survey and air dried. Here, we focus on the scats collected from one survey unit (unit O101) in the north-east during 2008.

Scat DNA extraction and evaluation for metabarcoding

To obtain DNA for metabarcoding, DNA was extracted from 10 scats (0.1–0.2 g/scat) plus one extraction negative by using the Invitrogen™ ChargeSwitch® forensic DNA purification kit following the manufacturer’s protocol. DNA extraction and PCR set-up were conducted in a dedicated trace DNA laboratory that was located in a PCR product-free building and subjected to strict quarantine protocols to minimise the risk of contamination.

We used quantitative PCR (qPCR) to evaluate amplification success from serial dilutions of DNA from each sample, to provide an indication of expected metabarcoding success and the likely level of PCR inhibition in each sample (Ramsey et al. 2015). We used the ViiA7 Real-Time PCR System (Thermo Fisher Scientific) to screen three replicate PCRs (by using undiluted DNA extract, 1/10 dilution, and 1/100 dilution as the templates) of each DNA sample. Each reaction had a total volume of 25 μL, including 2.5 μL of 10 × Gold buffer, 50 μM MgCl₂, 10 mg/μL bovine serum albumin (BSA), 10 μM each of the relevant forward and reverse (either 12SV5 or 16SMam) primers, 0.6 μL of 1/2000 dilution of SYBR Green, 6.3 μM dNTPs, 0.2 μL Amplitaq Gold, H₂O and 2 μL of the appropriate DNA dilution. PCRs were denatured at 95°C for 5 min and cycled through 40 cycles of 95°C for 30 s, 57°C for 30 s, 72°C for 30 s, followed by denaturing at 95°C for 15 s and a melt curve from 60°C to 95°C at a rate of +1°C/5 s. We did not use a species-specific blocking primer targeted to block the predator DNA because we did not know, a priori, for any given scat, which species of predator each scat came from. Three negative controls containing PCR water instead of DNA were included in each run. We conducted a comparative CT analysis using the ViiA7 Real-Time PCR System’s software v1.2.4, (Thermo Fisher Scientific, https://www.thermofisher.com/au/en/home/life-science/pcr/real-time-pcr/real-time-pcr-instruments/viia-7-real-time-pcr-system/viia-7-software.html), with a threshold of 40,000 ΔRn determined automatically. We calculated the CT value for each dilution of each DNA sample and used this to select the most appropriate DNA dilution to use for generation of DNA metabarcoding libraries. Samples were excluded from further analysis if all replicates failed to amplify or if all replicates had a CT value higher than 39.

Metabarcoding-library preparation

Two metabarcoding libraries, namely one for the 12SV5 marker and one for the 16SMam marker, were constructed for each DNA extract and for each DNA extraction negative control. Library construction involved a single qPCR step using fusion-tagged primers for the relevant markers. For each marker, forward primers included, from 5′ to 3′, the Illumina P5 adapter, a custom MiSeq sequencing primer, a multiplex identifier (MID) sequence tag, and the gene-specific forward primer. Reverse primers included, from 5′ to 3′, the Illumina P7 adapter, a multiplex identifier (MID) sequence tag, and the gene-specific reverse primer. Molecular identification tags were used to identify libraries from individual scats, so that these could later be pooled. Reactions of 25 μL total volume included 2.5 μL of 10 × Amplitaq Gold buffer (Applied Biosystems), 50 μM MgCl₂ (Applied Biosystems), 10 mg/μL BSA (Invitrogen), 20 μM each of the appropriate forward and reverse primers (Integrated DNA Technologies), 0.6 μL of 1/2000 dilution of SYBR Green I (Invitrogen), 6.5 μM dNTPs (Thermo Scientific), 0.2 μL Amplitaq Gold (Applied Biosystems), nuclease-free water (Integrated DNA Technologies) and 3 μL of DNA (using either undiluted DNA extract, 1/10 dilution, or 1/100 dilution as determined during the evaluation described above). qPCR conditions consisted of denaturation at 95°C for 5 min, 40 cycles of 95°C for 30 s, 57°C for 30 s and 72°C for 2 min, extension at 72°C for 10 min, denaturation at 95°C for 15 s, then a melt curve from 60°C to 95°C at a rate of +1°C/5 s, with a final hold at 4°C. We used the ViiA7 Real-Time PCR System (Thermo Fisher Scientific) and conducted comparative CT analyses by using the ViiA7 Real-Time PCR System’s software v1.2.4, with a threshold of 40,000 ΔRn determined automatically. We calculated the CT value for each library and submitted libraries for sequencing if they had a CT value lower than 39. Library construction was repeated for samples where amplification initially failed at only one gene. Samples, including DNA extraction negative controls, for which amplification failed at both genes were excluded from further analysis.

Amplicon libraries were pooled and cleaned in two stages. Initially, we pooled equal volumes of up to 10 libraries of similar concentration (as indicated by similar CT values); at this stage libraries from the two markers were not mixed. Each of these amplicon pools was cleaned using Ampure XP beads (Beckman Coulter), as described by (Murray et al. 2013), before being visualised on a 2% Tris–acetate–EDTA (TAE) agarose (Sigma Aldrich) gel stained with SYBR Safe (Invitrogen). A proportion of each library pool was then added to a final super pool, normalised according to its relative strength on the agarose gel. At this stage, libraries from the 12SV5 and 16SMam markers were mixed. The final super pool was cleaned using Ampure XP beads by using the same protocol, before equencing using the Illumina MiSeq platform, conducted at the Ramaciotti Centre for Genomics at UNSW.

Analysis of scat DNA metabarcoding data

We used Trimmomatic v0.36 (https://github.com/usadellab/Trimmomatic, Bolger et al. 2014) to remove adapter sequences, to trim poor-quality reads by using the sliding-window approach and default thresholds. We used Obitools v1.01 (https://pythonhosted.org/OBITools/introduction.html, Boyer et al. 2016) to cluster and identify reads. The ngsfilter command was used to demultiplex reads and assign each read to a gene and sample of origin according to the unique combination of MID tags and primers. We used the obistat command to remove possible chimeric sequences that are caused by PCR and sequencing errors (Schnell et al. 2015). We applied a threshold on sequence length taking only 12S rRNA sequences of greater than 80 bp and 16S rRNA sequences that were greater than 100 bp. Sequences were then compared with those in our Genbank reference-DNA database and our custom reference database of DNA sequences, by using the default settings in Obitools, and assigned to family, species and/or genus as appropriate. All reads that could be assigned to a family, genus or species of origin were compiled into a gene table for each marker, and then reviewed manually. Any sequence with only one read per scat was removed at this stage.

For each scat sample, the expected predator of origin was determined on the basis of results for each marker independently and then compared. Where results were ambiguous at the species level, the predator was assigned at the family or genus level as appropriate. Species inferences were based on the percentage of scats in which the prey was found by using presence or absence as the criterion because other methods based on the number of sequence reads do not necessarily reflect the intake of prey (Pompanon et al. 2012). Sequence reads from presumed prey items were identified to species level where possible. When these reads could be assigned to order or family only at one marker, results were cross-checked against reads observed at the other marker and against similar reads on GenBank, in case the identification could be resolved further. Results from the genetic distance-based evaluations described above suggest that a 1% genetic distance threshold would be sufficient to identify unknown sequences to species level for Diprotodontia at the 16S and 12S markers used here. Consequently, in basic local alignment search tool (BLAST) analyses, any unassigned scat DNA sequence that matched a GenBank sequence with an identity score of ≥99% and an E-value of ≤10–5, was assigned to that species, whereas the genus- or family-level identifications were retained for scat sequences that matched with identity scores of <99%. This reduced the probability that inaccurate species assignments remained in the dataset, because an inaccurate sequence is not expected to match a reference species at ≥99% by chance (Yamamoto et al. 2017). Finally, the tables for each gene were combined to produce one set of predator and prey sequences for each scat. Predators were assigned according to the highest number of reads in a sample and then removed from the prey column. Where reads from more than one predator species were identified from the same sample, the species with the highest number of reads was assumed to be the predator responsible for the scat, and species with lower read counts were included for analysis as potential prey.

Ethics statement

A license to collect road-kill in Tasmania was issued by Tasmania’s Department of Primary Industries, Parks, Water and Environment (Permit No. 13240).

Evaluation of primer utility for amplification of mammal DNA

The COI primers uni-minbar were developed as universal eukaryote primers and have previously been tested on eutherian mammals and three species of South American, but not Australian, marsupials (Meusnier et al. 2008). Here, we used PrimerTree to identify the presence of uni-minibar primer binding sites in different classes of animal, by comparison to all published sequences from the COI gene. We observed that these primers have the potential to amplify sequences from a broad range of taxa, including some mammals and many invertebrates, but were unsuitable for application to birds, with few primer binding sites identified from COI sequences from the class Aves (Fig. 1a).

Fig. 1.

Phylogenetic trees generated using the R package PrimerTree (Cannon et al. 2016), showing: (a) COI sequences with suitable binding sites for the uni-minibar primers; (b) 12S rRNA sequences with suitable binding sites for the 12SV5 primers; (c) 16S rRNA sequences with suitable binding sites for the 16SMam primers; and (d) 16S rRNA sequences with suitable binding sites for the L2513–H2714 primers. Sequences are resolved at the level of class. These trees include only those sequences that are likely to be amplified using the relevant primers.

The PrimerTree analysis provides an indication of the breadth of taxa likely to be amplified by a candidate primer set, but it provides no indication of the proportion of species that are likely to be amplified within each taxonomic group. We used PrimerMiner-0.12 to evaluate expected amplification success for each candidate marker within a desired group, in this case the infraclass Marsupialia, by evaluating the conservation of primer binding sites in all published marsupial sequences. This allowed us to identify any potential taxonomic bias in amplification for each marker. Conservation of primer sites in a high proportion of species by using a lower threshold was desirable because higher thresholds allow a greater number of base-pair mismatches within the primer binding region.

These results show that the COI uni-minbar marker is a poor marker for marsupial metabarcoding (Fig. 2a). We observed very low levels of primer binding-site conservation among published sequences from marsupial species, with only 6% of sequences from the order Diprotodontia, 10% from Didelphimorphia and 50% from Peramelemorphia considered to be likely to be amplified by this marker by using the least-stringent threshold of 200. Further, it appears that these primers would be unlikely to successfully amplify any sequences from the Dasyuromorphia, Microbiotheria, Paucituberculata or Notoryctemorphia. For these reasons, we did not consider the COI uni-minibar marker further in our analyses.

Fig. 2.

Bar plots showing the proportion of OTUs within each marsupial order for which primer binding sites are conserved for four candidate mitochondrial markers: (a) the COI uni-minibar primer sites showed extremely low conservation in any marsupial group; (b) primer sites were highly conserved across all marsupial orders for the 12SV5 marker; (c) primer sites were highly conserved across all marsupial orders for the 16SMam markers; (d) primer sites were conserved for some marsupial orders for the L2513–H2714 primers, but these would be unlikely to amplify DNA from Microbiotheria and Paucituberculata. Higher thresholds allow a greater number of nucleotide mismatches between primer and target sequences. Plot is not shown for the L1085–H1259 12S rRNA primers because we were unable to detect any published marsupial sequences that are likely to be amplified by these primers.

Only one of the 12S rRNA marker tested, 12SV5 (which was designed to be a universal marker for vertebrate species), was identified as broadly applicable by using PrimerTree (Fig. 1b). This analysis demonstrated that the 12SV5 primers are likely to amplify DNA from Mammalia, Amphibia, Aves and Actinopterygii, and therefore identified these primers as suitable for vertebrate metabarcoding. We were unable to generate a reliable tree for 12S rRNA primers L1085–H1259, even though initial testing in Geneious v10.0.5 suggested that the marker could be a suitable candidate for amplification of marsupial DNA.

We used PrimerMiner to evaluate both sets of 12S rRNA primers for all orders in the infraclass Marsupialia. Primer binding sites were highly conserved across all marsupial orders for the 12SV5 marker (Fig. 2b), with amplification expected from >98% of sequences from Dasyuromorphia and Diprotodontia, and 100% of sequences from all other marsupial orders by using the most stringent threshold of 50. This indicated that the 12SV5 marker is a suitable candidate for eDNA metabarcoding studies of Marsupialia. Despite initial indications that the L1085–H1259 marker would be a suitable candidate for amplification of marsupial DNA, we were unable to detect any published marsupial sequences that are likely to be amplified by these primers (Fig. 2c). For this reason, the L1085–H1259 marker was excluded from further consideration.

Our PrimerTree analysis demonstrated that both primer sets targeting the 16S rRNA gene are likely to reliably amplify mammal DNA. Binding sites for the L2513–H2714 primers were conserved in DNA sequences originating from Mammalia, Amphibia and Actinopterygii (Fig. 1d), whereas binding sites for the 16SMam primers were observed in DNA sequences from Amphibia and Mammalia (Fig. 1c).

PrimerMiner analysis of the 16SMam marker showed that primer binding sites were highly conserved across all orders in the infraclass Marsupialia (Fig. 2d), with amplification expected from >95% of sequences from Diprotodontia, >85% of sequences from Dasyuromorphia, and 100% of sequences from all other marsupial orders by using the most stringent threshold of 50. Similarly, the L2513–2714 primer binding sites were highly conserved within four of the seven marsupial orders (Fig. 2e), with potential for amplification of 100% of sequences at the least stringent threshold. This marker is also likely to be amplified from over 98% of sequences from the order Dasyuromorphia by using a threshold of 100, and from over 90% of sequences by using the most stringent threshold of 50. However, two of the three South American marsupial orders, Microbiotheria and Paucituberculata, exhibited poor conservation of primer binding sites and would be unlikely to amplify using the L2513–2714 marker. This result suggests that the marker L2513–2714 would be suitable for application to Australian but not South American marsupials. We elected to exclude this marker from further evaluation in this study, primarily because the longer amplicon length relative to the 16SMam marker makes this marker less suitable for analysis of degraded or eDNA. We selected the 16SMam marker for further evaluation because it has a shorter amplicon length (the primers were developed specifically for application to highly degraded DNA), it has broad applicability within Marsupialia, and it displays a higher utility for mammal DNA overall.

Genetic distance-based evaluation of the candidate metabarcoding markers

We constructed custom reference-DNA databases for all Tasmanian Diprotodontia and some mainland Australian Diprotodontia at the 12SV5 and 16SMam markers. We sequenced 58 tissue samples at the 16S rRNA gene and 164 tissue samples at the 12S rRNA gene. Gaps in these databases were filled with published sequences obtained from GenBank for both the 16S rRNA (n = 15) and 12S rRNA (n = 35) genes (Table S1). Singleton species were then removed, leaving a total of 194 sequences from 19 different species and 10 genera in the 12SV5 database, and 64 sequences from 11 different species and eight genera in the 16SMam database.

We used the R package SPIDER to conduct genetic distance-based analyses, which demonstrated that the 12SV5 and 16SMam markers can both identify marsupial sequences at a high level of resolution. All 64 sequences in the 16SMam database were correctly identified using a nearNeighbour analysis. In the bestCloseMatch analysis, 59 (using a 1% genetic distance threshold) or 63 (3% and 5% threshold) of the 64 sequences were correctly identified (Table S2). The remaining sequences could not be identified because no other sequences in the database were within the specified genetic-distance threshold. However, on closer inspection by using higher genetic-distance thresholds, the nearest neighbour to each of these sequences originated from a conspecific. This makes clear the importance of using a reference database that fully captures local intra-specific variation. In the threshID analysis, we again observed 59 correct identifications by using a 1% genetic-distance threshold. The same five sequences that could not be identified in the bestCloseMatch analysis remained unidentified in the threshold-identity analysis. Only 49 or 47 sequences were correctly identified using thresholds of 3% and 5% respectively (Table S3). No sequences were incorrectly identified and, as before, one sequence could not be identified because all other sequences in the database were outside the selected genetic-distance threshold.

We observed greater opportunity for ambiguity in identification when using the threshID analysis, which considers all sequences within the genetic distance threshold as potential matches, rather than just the closest match. This occurs when there is no clear barcoding gap between closely related species, typically because sequences from conspecifics and from at least one other species (usually a sister taxon) are within the specified-distance threshold. Unfortunately, location data were not available for several sequences obtained from GenBank, making it difficult to determine whether samples that could not be identified, or were ambiguously identified, may represent genetically distinct populations or subspecies, as might be expected for Tasmanian versus mainland Australian specimens. In these circumstances, and if data interpretation also considered gaps in coverage of intra-specific variation in the reference database, we considered it appropriate to recommend a bestCloseMatch analysis with a 3% threshold, for the 16SMam marker for Diprotodontia.

All 194 sequences in the 12SV5 database were correctly identified using a nearNeighbour analysis. In the bestCloseMatch analysis, 186 (using a 1% genetic distance threshold), 192 (3% threshold) or 193 (5% threshold) sequences were correctly identified (Table S4). Eight sequences (1% threshold) and two sequences (3% threshold) could not be identified because no other sequences in the database were within these specified-distance thresholds. Furthermore, a single sequence from Potorous tridactylus (sequence AY245617) was ambiguously identified (by using a 5% threshold) as being within a 5% genetic distance of sequences from other Potorous species as well as sequences from Bettongia penicillata and Thylogale thetis. No information is available regarding the origin of the sample from which this sequence was obtained (Westerman et al. 2004). This species is found in Tasmania and on mainland Australia, and some genetic differentiation has been observed among individuals from different geographic regions (Frankham et al. 2012). Consequently, this ambiguity could reflect regional genetic variation or mistaken identity. In the threshID analysis of the 12SV5 database, by using a 1% threshold, 182 sequences were correctly identified, and eight sequences could not be identified (these were the same eight sequences that could not be identified in the corresponding bestcloseMatch analysis). The remaining four sequences (two from Bettongia lesueur and two from Cercartetus nanus) were listed as incorrectly identified. However, we found that in each case no other sequence was within the 1% genetic distance threshold, indicating that ‘no identification’ would be a more appropriate result for these samples by using these parameters. Using a 3% threshold, each of these sequences clustered only with sequences from conspecifics, which suggests that incorrect identification of these samples would be unlikely. As observed for the 16SMam analysis, sampling location information was not available for some GenBank sequences, making it difficult to determine whether these sequences represent isolated populations or different subspecies. On the basis of our results, it again seems most appropriate to use the bestcloseMatch analysis, but this time with a 5% threshold, for the 12SV5 marker for Diprotodontia. However, care must be taken when considering identifications within the Macropodiformes, particularly the genus Potorous, which appears to display lower genetic differentiation at this marker. In this case, the combined use of two genetic markers for DNA metabarcoding would improve confidence in identification of species known to be at risk of misidentification or ambiguity.

Phylogenetic evaluation of the candidate metabarcoding markers

The phylogenetic tree for the 92 bp 16SMam amplicon demonstrated that all species included could be distinguished and grouped as expected according to their recognised taxonomic relationships (Fig. 3). All but two bootstrap values for differentiation within a genus were strongly supported (>90), with the exceptions being Notamacropus rufogriseus (50) and Macropus giganteus (49). All other bootstrap values between genera were well supported (>70), with the majority (seven) having very strong support (>90). Previous research identified two clades within Thylogale billardierii in Tasmania (Macqueen et al. 2010), and our results are consistent with this despite the much shorter sequence length used in this study. Published data at the 16S rRNA gene were limited for most extant species of Diprotodontia that have distributions covering both Tasmania and mainland Australia. This made it difficult to determine whether Tasmanian animals could be reliably distinguished from mainland Australian animals by using the 16SMam marker. In contrast to our results for the 16SMam amplicon, phylogenetic resolution was poor using the 12SV5 marker (Fig. 4). Specifically, the tree topology did not reflect recognised taxonomy and some bootstrap values were very low.

Fig. 3.

Maximum-likelihood tree for Tasmanian and some mainland Diprotodontia, with Dasyrus as an outgroup, at the mitochondrial 16S rRNA gene using the 16SMam marker. The evolutionary history was inferred by using the maximum- likelihood method based on the Tamura 3-parameter model. A discrete gamma distribution was used to model evolutionary rate differences among sites. The tree is drawn to scale, bootstraps shown at branch nodes.

Fig. 4.

Maximum-likelihood tree for Tasmanian and some mainland Diprotodontia, with Dasyurus as an outgroup, at the mitochondrial 12S rRNA gene using the 12SV5 marker. The evolutionary history was inferred by using the maximum-likelihood method based on the Kimura 2-parameter model. A discrete gamma distribution was used to model evolutionary rate differences among sites. The tree is drawn to scale, with bootstrap values at branch nodes.

Field testing of candidate metabarcoding markers

We applied our markers to 10 scats collected from a single survey unit (O101) from the north-east of Tasmania in 2008. Our metabarcoding analysis using the 12SV5 and 16SMam markers identified DNA from the following three predators among the scats (Fig. 5): Tasmanian devil (Sarcophilus harrisii, n = 4), feral cat (Felis catus, n = 3) and spotted-tail quoll (Dasyurus maculatus, n = 3). We also identified a total of 13 discrete prey taxa, including six species of marsupial, three species of bird, a lizard from the Scincidae family, and three common introduced animals of cats, dogs, and the black rat (Rattus rattus). However, we note that the methods used were unable to distinguish consumption of cats and dogs by devils from other potential interactions, for example, if a cat or dog were to overmark or otherwise come into direct contact with a devil scat. There were no obvious omissions of common taxa in our samples; however, further surveys would be required to determine the detectability of specific rare taxa. The resolution of prey species was high for Diprotodontia and Rodentia, but because our reference database was deficient in Tasmanian bird and reptile species, we were unable to assign all sequences to species level. We also note that spotted-tail quolls have been reported to eat large numbers of invertebrates at times (Fancourt et al. 2018). Our markers are unlikely to detect members of this class.

Fig. 5.

Prey species detected in scats collected from Grid square O101 in north-eastern Tasmania in 2008. Bar graph indicates the percentage of scats in which each putative prey species was detected for each predator.