Samples

We collected 276 faecal samples from 11 roosts of the little red flying-fox (Pteropus scapulatus), three roosts of the black flying-fox (Pteropus alecto), two of which overlapped with P. scapulatus, and one roost of the spectacled flying-fox (Pteropus conspicillatus) across Queensland between July 2017 and November 2019 (Fig. 1; Supplementary Table S1). Faecal samples were taken directly from captured bats using a sterile swab, and new plastic sheets were placed on the ground below roosts to collect additional faecal samples without handling animals. Faecal samples were generally collected before 12 noon to ensure adequate quantity of fresh sample material. However, this was not always possible in remote roosts. Faecal pellets and swabs were placed in sterile 1.5 mL microcentrifuge tubes with 95% ethanol to preserve plant DNA and to neutralise any pathogens. Negative controls were also taken in the field, and comprised a sterile swab placed into a sterile 1.5 mL tube of ethanol at the sample collection site. Full details of samples can be found in Supplementary Table S1.

Fig. 1. 

Faecal sample collection locations.

Leaf samples were collected from any tree species in flower or fruit near sampling sites, to obtain reference sequences of likely food species. Leaves were dried and preserved in sealed bags of silica gel. Species were identified by CSIRO botanists. Full details of leaf samples can be found in Supplementary Table S2.

DNA extraction

Samples were preprocessed by vortexing to dislodge as much of the sample from the swab as possible and to homogenise the sample, then centrifuged at 15000g for 2 min, after which the ethanol was set aside in a separate tube. Up to 200 mg of the sample was then transferred to a new tube containing 200 µL of ultrapure water, scraping the swab against the side of the tube to transfer as much of the sample as possible. The remaining sample was reunited with the original ethanol and stored at –80°C, to retain any residual DNA for future work.

Following a preliminary method comparison between two different kits (see Supplementary Methods and Results, Table S3), DNA was extracted using the Nucleospin DNA Food Kit (Macherey-Nagel, Düren, Germany), following the manufacturer’s protocol with the following modifications: samples were homogenised using new, sterile 3 mm tungsten carbide beads for 1 min in a TissueLyser II (Qiagen, Hilden, Germany). Following homogenisation, we then followed the manufacturer’s protocol for extraction of DNA using the NucleoSpin DNA Food Kit except the volume of Elution Buffer CE was decreased to 50 µL to increase the concentration of the final extraction. A negative control (no faecal sample added) was included in each batch of 12–24 DNA extractions. Leaf DNA extraction was conducted using the DNeasy Plant Kit (Qiagen). Dry leaf was homogenised in a TissueLyser II for 1 min with two 3 mm tungsten carbide beads, cleaned with 0.4m HCl in between samples. This step was repeated 2–3 times until the leaf was powdered. If the leaf was still not powdered after three repetitions, it was ground using a sterilised mortar and pestle. We extracted DNA from the dry, powdered leaf following the manufacturer’s protocol. We used only the DNA from the first 100 µL elution to maximise the concentration of DNA going into PCR. A negative control was included in each batch of DNA extractions.

To ensure that our methods were able to accurately detect and identify plant species in mixtures of similar species-richness to faecal samples, we created mock communities. Each mock community contained five plant species, from leaf DNA extractions, in approximately equal concentrations. Mock Community 1 comprised Ficus racemosa (Moraceae), Bauhinia hookeri (Fabaceae), Eucalyptus platyphylla (Myrtaceae), Corymbia clarksoniana (Myrtaceae), and Melaleuca leucadendra (Myrtaceae). Mock Community 2 comprised Syncarpia glomulifera subsp. glomulifera, Eucalyptus camaldulensis, Melaleuca bracteata, Melaleuca viridiflora, Melaleuca nervosa (all Myrtaceae). Mock Community 3 comprised Eucalyptus platyphylla (Myrtaceae), Corymbia torelliana (Myrtaceae), Brachychiton bidwillii (Malvaceae), Terminalia playtphylla (Combretaceae), and Melaleuca viridiflora (Myrtaceae). Mock Community 4 comprised Eucalyptus camaldulensis (Myrtaceae), Eucalyptus tereticornis (Myrtaceae), Alyxia ruscifolia (Apocynaceae), Syzygium suborbiculare (Myrtaceae), and Eucalyptus megasepala (Myrtaceae). All of these species were present in the reference databases for both rbcL and ITS2.

DNA sequencing

The concentration of each faecal or leaf DNA extraction was determined using the Qubit 1X dsDNA HS Assay Kit on a Qubit 4 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). Samples with a concentration greater than 20 ng/µL were diluted to 10–20 ng/µL before amplification and sequencing.

Two DNA barcode markers were used, ITS2 and rbcL. Following a preliminary trial of two different primer pairs for each marker (Supplementary Methods and Results, Table S4), we conducted PCR amplification using the primers in Table 1. The two DNA barcode markers were amplified in separate reactions for each sample. We used an indexing strategy that gave each sample a unique combination of two indices, that were used for both DNA barcodes. The PCR reactions contained primers at a final concentration of 200 nm, 12.5 μL of KAPA HiFi ReadyMix (KAPA Biosystems, Boston, MA, USA), and 8.5 μL of the template DNA (up to 20 ng/μL) in a 25 μL reaction. To increase the chance of detecting all species in the mixture, each DNA extraction was included in three PCR reactions that were amplified separately, i.e. the reaction contents described above were divided between three PCR tubes. The PCR conditions included an initial period of heat activation for 3 min at 95°C; followed by 35 cycles of 30 s at 95°C, 30 s at 55°C, and 1 min at 72°C; followed by a final extension of 10 min at 72°C and then held at 10°C. For the short fragments, following a gradient PCR to determine optimal annealing temperature, PCR was conducted with an initial denaturing step of 95°C for 3 min; followed by 35 cycles of: 95°C for 30 s, 55°C for 30 s, 72°C for 1 min; and a final extension at 72°C for 10 min. A negative control (with ultrapure water as DNA template) was included in each PCR plate. We used our constructed mock communities as positive controls, and these were amplified in duplicate with different indices.

Table 1. 

PCR and sequencing primers used in this study

PCR product lengths include the fusion primers. Illumina tag is shown in bold; index, pad and linker sequences in plain text; and marker-specific primer in italics

After amplification, PCR products were purified with Agencourt AMPure XP magnetic beads (Beckman Coulter, Danvers, MA, USA) to remove any remaining primers and PCR artefacts. The presence of a PCR product of the correct size was assessed by running 2 µL on a QIAxcel (Qiagen), using the 15bp–3kb alignment marker and 50bp–800bp size marker. The DNA concentration was quantified using the Qubit 1X dsDNA HS Assay Kit on a Qubit 4 Fluorometer. Within each sequencing run, all PCRs were pooled at equal concentrations and run on a single flow cell of the Illumina MiSeq instrument by AGRF (Melbourne), using a 2 × 250-cycle paired-end run for the first sequencing run, and a 2 × 300-cycle paired-end run for subsequent batches. Each sample had an index combination that was unique across the entire project to avoid cross-contamination between sequencing runs. A subset of eight index combinations (two for each batch) was not used for a sample, and was demultiplexed by AGRF as a sequencing negative control.

Taxonomic identification

For each sequenced sample, the ITS2 and rbcL sequencing reads were sorted into separate fastq files using scripts provided by Anya Cutler (Emory University, USA) for analysis in R (R Core Team 2016). This script uses the first four base-pairs of the expected sequence to identify reads as either ITS2 or rbcL and sort them into separate files. The dada2 package in R (Callahan et al. 2016) was used to remove likely sequencing errors and generate a set of unique DNA sequences. For the subset of samples in the preliminary run, analysed from long amplicons with 250 bp paired-end sequences, the rbcL amplicon was too long for forward and reverse sequences to merge. Therefore, these samples were analysed based on the single reads without merging, and in subsequent runs we used 300 bp paired-end sequences. Sequences that were likely to be the result of contamination were flagged using the decontam package (Davis et al. 2018) in R, using the prevalence method which compares presence/absence of sequences in samples and negative controls. These were removed from further analyses. To determine if adequate sequencing depth had been used, species accumulation curves were generated for a subset of 10 samples from one sampling event (P. scapulatus, Charters Towers, September 2017) using the vegan package in R (Oksanen et al. 2019), with Amplicon Sequence Variants (ASVs; unique DNA sequences) in place of species and number of sequence reads as units.

Sequences from leaf samples were separated from the dada2 output. Unique ITS2 sequences that were represented by at least 10% of the total sequence reads in the sample were used as reference sequences for taxonomic identification, as these may represent multiple paralogues, or alternate diploid alleles. Only one rbcL sequence from each leaf sample was retained as a reference sequence, as this locus is haploid and not expected to have paralogous copies. Remaining sequences were compared with the NCBI database ( www.ncbi.nlm.gov) using a BLAST search to confirm that they were the designated species and not contaminants. For rbcL, we appended our new sequences to the database of Bell et al. (2017b), and formatted the database for the dada2 functions assignSpecies and assignTaxonomy. The revised rbcL database is available from  https://doi.org/10.6084/m9.figshare.c.5504193.v1. A new reference database was generated containing all publicly available ITS2 sequences. Flowering plant ITS2 sequences were obtained from the NCBI database using the following search terms: ‘Magnoliopsida’[Organism] AND ITS2[All Fields]. The program ITSx was used to extract 5.8S+ITS2 sequences from longer sequences (Bengtsson-Palme et al. 2013). Publicly available reference databases, such as the NCBI database, rely on researchers to perform their own quality control before submitting sequences, and contaminating sequences are known to occur in the database (Breitwieser et al. 2019). To identify any erroneous sequences that could cause problems for our analysis, we conducted BLAST searches of the 50 most abundant ITS2 ASVs and 20 most abundant rbcL ASVs from our eDNA sequencing. Sequences were removed from our reference database if they were in the top 200 matches for one of these ASVs, but with taxonomy fields on NCBI that were outliers relative to all other sequences in the top 200 matches (ITS2: FM887018 Terminalia bellerica; rbcL: JX856776 Strychnos potatorium, KU564752 Actephilia sessifolia, KU564758 Antirhea putaminosa, MH598849 Datura stramonium, JX856673 Callistemon polandii, KU564753 Alectryon pubescens, KU564797 Eupotamia bennettii, KU564800 Freycinetia excelsa, KU564801 Freycinetia scandens, and KU564874 Sophora fraseri). There are likely other erroneous sequences remaining in our database, but these are not a close match with our high-frequency ASVs, and so are unlikely to influence our results. The database was then formatted for the dada2 functions assignSpecies and assignTaxonomy, and our reference leaf sequences were appended to this database. The ITS2 database is available from  https://doi.org/10.6084/m9.figshare.c.5203703.v1.

Matrices of DNA sequence by sample generated by dada2 analysis of each sequencing run were merged into a single data matrix for each marker (with the unmerged rbcL long reads from the preliminary sequencing run treated as a different marker, due to the inability to merge reads) following the dada2 instruction manual. Samples containing fewer than 1000 total sequencing reads were removed from further analysis. Identification to higher taxonomic level (up to genus) was conducted using the assignTaxonomy function in dada2. Species-level identification was determined using the assignSpecies function, which requires exact matches and can output multiple exact matches if present. Because ITS2 sequences from the same species typically vary by up to 1% (Chen et al. 2010), we also considered any ITS2 sequence on the NCBI database or among our reference sequences with at least 99% identity to be a possible species-level match. Only exact matches were considered for the less variable rbcL barcode, and where multiple exact matches occurred, all were considered plausible. Finally, for sequences with multiple plausible species-level identifications (multiple matches with at least 99% identity with ITS2 or multiple exact matches with rbcL), based on sequence similarity, species distributions and, where appropriate, flowering and fruiting records and observations were used to shorten the list to the most likely species. Specific observation of flowering and fruiting allowed us to confidently identify the plant part in most collections, particularly in woodland settings where species dominance is high. In rainforest settings with low species dominance, this is more difficult.

DNA sequencing results

The preliminary run of Illumina MiSeq generated 2 304 061 paired-end 250 bp sequencing reads. Of these, 1 854 641 were from the 25 samples sequenced with long PCR products. After separating sequences into ITS2 and rbcL reads, and analysis via the dada2 bioinformatics pipeline, there remained 456 195 long ITS2 reads, 139 295 long rbcL forward reads, and 139 152 rbcL reverse reads. Two subsequent Illumina MiSeq runs of long ITS2 and long rbcL PCR products generated 11 608 382 and 11 490 757 paired-end 300 bp sequencing reads. After separation of ITS2 and rbcL sequences, filtering and cleaning with dada2, and contaminant removal with decontam, 3 129 419 and 3 261 146 ITS2 sequences and 1 973 615 and 3 579 131 rbcL sequences remained, for each sequencing run respectively (Supplementary Table S5). Of these, 99.3% of ITS2 sequences were able to be identified to family and 98.9% to genus with the assignTaxonomy function, while 37.9% were identified to exactly one species with the assignSpecies function. The rbcL sequences tended to have lower taxonomic resolution, with 99.4% able to be identified to family, 84.7% to genus, and 5.9% to exactly one species (Fig. 2).

Fig. 2. 

Taxonomic level of identification with ITS2 and rbcL. Percentage of Illumina MiSeq reads across all samples identified to each taxonomic level using long DNA sequences (350–450 bp of ITS2; ∼550 bp of rbcL), with taxonomic matching conducted using the assignTaxonomy and assignSpecies functions in dada2.

The decontam function identified 18 and 13 ITS2 sequences as contaminants in sequencing runs 2 and 3, respectively, representing 64 728 and 78 371 sequencing reads (2.1 and 2.4% of input reads). For rbcL, 4 and 6 sequences were identified as contaminants, representing 60 802 and 50 793 sequencing reads (3.1 and 1.4% of input reads), respectively. These sequences were removed from subsequent analyses, but, due to their low proportions in mixtures, they would have been unlikely to affect our results if they were retained.

Species accumulation curves of a subset of 10 samples, considering ASVs as species and sequencing reads as sampling units, showed most ASVs within a sample were detected with fewer than 5000 sequencing reads (Fig. 3). For many samples, this could be achieved with 1000 sequencing reads or fewer.

Fig. 3. 

Species accumulation curves of amplicon sequence variants (ASVs) from each sample, after quality control and filtering with dada2, as a function of the number of sequencing reads. (a) ITS2; (b) rbcL.

Taxonomic identification

Species known to be present in mock communities were almost always able to be detected with both ITS2 and rbcL, but were not always identified to the species level, despite all being present in reference databases (Fig. 4). Using ITS2, all species were detected except for Ficus racemosa (Moraceae) in Mock Community 1 and Melaleuca nervosa in Mock Community 2. We failed to amplify F. racemosa ITS2 on its own during the reference database supplementation step, which suggests there could be a primer mismatch. Sequence of ITS2 is available for this species on NCBI, so we were still able to include it in the reference database. Melaleuca nervosa was likely indistinguishable from M. viridiflora, which was also present, rather than being undetectable. Eucalyptus camaldulensis and E. tereticornis were indistinguishable but detected in all samples where present. Using rbcL, all species were detected except for Corymbia clarksoniana in Mock Community 1. Brachychiton bidwillii (Malvaceae) in Mock Community 3 was detected but misidentified as Ochroma sp. (Malvaceae). Eucalyptus camaldulensis, E. tereticornis and E. platyphylla were indistinguishable but detected where present. Likewise, M. leucadendra, M. nervosa and M. viridiflora were indistinguishable but detected where present.

Fig. 4. 

Taxonomic identifications of mock communities. Four mock communities of known species composition, as described in the text, were sequenced using long amplicons of ITS2 and rbcL and the sequences were identified using the assignTaxonomy and assignSpecies functions in dada2.

Species and genera from the family Myrtaceae were the predominant components identified in faecal samples. The most abundant ITS2 ASV in faecal samples was identified as Eucalyptus camaldulensis (Myrtaceae). Other abundant ASVs were identified as Eucalyptus tereticornis, Melaleuca viridiflora (Myrtaceae), Eucalyptus sp., and Melaleuca leucadendra (Table 2). The most abundant rbcL ASV was identified as Melaleuca sp. Other abundant ASVs were identified as Eucalyptus sp., Myrtaceae sp., Bauhinia sp. (Fabaceae), and Corymbia sp. (Myrtaceae). Some ASVs matched plant species that were wind-pollinated, mostly in the Poaceae and Casuarinaceae. These were likely contaminants from the field or laboratory and were usually present only in a single sample. The most frequently detected wind-pollinated species was Axonopus (Poaceae), which was identified from ITS2 in 11 samples and from rbcL in 14 samples. A full list of ITS2 and rbcL taxonomic identifications from faecal samples, can be found in the Supplementary Material (Tables S5 and S6, respectively). These data, along with data from additional faecal samples, is being used to understand flying-fox foraging behaviour and the results of this work will be published elsewhere.

Table 2. 

Taxonomic identification of commonly occurring amplicon sequence variants (ASVs; unique DNA sequence variants) in flying-fox faecal samples

Taxonomic identifications were based on a combination of the following method(s): (1) the taxonomic classifications generated by dada2 using the assignTaxonomy and assignSpecies functions, (2) 99% sequence identity at the species level to sequences on NCBI and our reference sequences, and/or (3) knowledge of geographic range and flowering phenology of the species

The taxonomic composition of a subset of samples is shown in Fig. 5. Within the same samples, a similar taxonomic composition was determined with both DNA barcodes. In most cases, the same genera were detected with both ITS2 and rbcL, but with ITS2 taxonomic identification was more likely to be to species-level, or a small subset of the species within a genus. Occasionally a taxon was detected with only one or the other marker. For example, Musa sp. (Musaceae) in 1709005 and Vitex sp. (Lamiaceae) in 1709010 were identified only from rbcL sequences, and Shotia brachypetala (Fabaceae) in 1709010 was identified only with ITS2 sequences.

Fig. 5. 

Taxonomic composition of a subset of flying-fox faecal samples. Faecal samples were sequenced using long amplicons of ITS2 and rbcL and the sequences were identified using the assignTaxonomy and assignSpecies functions in dada2.

Performance of standard pollen DNA metabarcoding methods

We adapted our method from a protocol previously developed for the identification of plant species in bee-collected pollen (Bell et al. 2017a). Faecal samples can contain enzymes involved in digestion, which may inhibit PCR. For this reason, we trialled an alternative DNA extraction kit, specifically designed for faecal samples. We extracted higher quantities of DNA with the method previously optimised for extraction of DNA from pollen, with no PCR failure. Faecal samples may also be expected to contain degraded DNA due to the digestion process. At least one previous study using DNA metabarcoding for diet analysis of a flying-fox used primers specifically designed to sequence a shorter DNA fragment than typical DNA barcodes (Aziz et al. 2017). However, we found that full-length DNA barcodes could be successfully amplified and sequenced. This is consistent with the results of Lim et al. (Lim et al. 2018a, 2018b) and Chan et al. (2021), who were able to successfully amplify and sequence faecal samples from nectivorous and frugivorous bats with full-length DNA barcodes. The gut passage time of flying-foxes is usually fast, although some seeds can be retained for up to 24 h (Westcott et al. 2001; Aziz et al. 2021). Fast gut passage times may mean that the DNA is less degraded than faecal DNA of other species. Lack of degradation could also be because the majority of DNA in the faecal samples was present within pollen grains.

We found that, in combination, ITS and rbcL have been able to detect all species in a mixture. This is consistent with the results of Bell et al. (2019), where mixtures of known pollen species composition were all able to be identified using these two DNA barcodes. This combination of barcodes has not been previously tested in north-eastern Australian study systems. The standard chloroplast markers rbcL and matK, in addition to trnH-psbA as a supplementary marker (also from the chloroplast genome) have been successfully used for DNA barcoding in south-east Queensland (Shapcott et al. 2015). As a study system dominated by Myrtaceae species, particularly Eucalyptus and related genera, we may have expected some difficulties with DNA barcoding. Previous studies have shown that distantly related Eucalyptus species can share chloroplast genomes due to hybridisation and introgression (Schuster et al. 2018). For this reason, we supplemented the more traditional chloroplast barcode rbcL, with the nuclear ribosomal barcode ITS2. However, DNA barcoding via nuclear rDNA could also be complicated, due to paralogous gene copies. Combining the DNA barcoding data with field observations of flowering, and fruiting where appropriate, and a knowledge of plant species distributions, allowed us to overcome these issues. For example, Eucalyptus tereticornis and E. camaldulensis partially overlap in range and flowering season and share ITS2 sequences. For most samples, assignment to species was possible based on the distribution of each species, then our flowering records.

Relative performance of ITS2 and rbcL

The taxonomic resolution was found to be higher for ITS2 than rbcL, with many ITS2 sequences being identified to species level (37.9%), and almost all (98.9%) identified to genus. On the other hand, only 5.9% of rbcL sequences were identified to species, and only 84.7% identified to genus. A higher rate of species-level identification with ITS2 compared with rbcL has been observed elsewhere in DNA metabarcoding of environmental samples, such as bee pollen (Bell et al. 2017a), and soil (Fahner et al. 2016). For single-species DNA barcoding, ITS2 has also been shown to have higher taxonomic resolution than rbcL (Chen et al. 2010; Li et al. 2011), especially if multiple sequences per species are included (Kolter et al. 2021), and even with shorter amplicons species-level resolution has been shown to be 86.1% (Moorhouse-Gann et al. 2018). For this reason, we recommend that both markers be used, or if only one marker is used that ITS2 be preferred over rbcL if high-level taxonomic resolution is required.

Food plants of Queensland flying-foxes

Species and genera from the family Myrtaceae were the predominant components identified in faecal samples. Frequently occurring species included E. camaldulensis, E. tereticornis, M. viridiflora, M. leucadendra, other Eucalyptus and Melaleuca spp., and Corymbia spp. Other studies have also shown flying-fox species to frequently include flowers of Myrtaceae species in their diets. For example, Parry-Jones and Augee (2001) found Myrtaceae pollen in the faeces of Pteropus poliocephalus from Sydney in all months of the year. This is likely due to the prevalence of Myrtaceae species in Australian vegetation, as well as the nectar and pollen resources provided by these species. The food plants identified in this study are examined in more detail in Bradford et al. (2020).

Shortcomings and priorities for future research

We validated our method using mock communities of known species composition, made from DNA extracted from leaf material, of species that were present in reference databases. Using both DNA barcodes in combination, all species could be detected. However, occasionally species were missed with one barcode or the other. If analysing only ITS2 data, it is possible that some species such as Ficus sp. may be missed, but this can be resolved by combining the ITS2 with rbcL sequences, or redesigning primers for ITS2 to improve amplification. Our mock community analysis demonstrated that if DNA is present in a faecal sample, then it can be detected with either or both of ITS2 and rbcL. However, we have not confirmed that there will be adequate DNA in faeces for all species eaten to be detected. When gut passage time is short, this could mean that only species consumed within the last few hours are detected. Finally, eDNA analysis alone cannot determine which plant part is being consumed. Combining eDNA analysis with microscopic analysis of faecal samples could assist with determining which plant part is consumed, although if plant parts are eaten in liquid form, such as fruit juice or nectar, these may not be detectable in the faeces (Aziz et al. 2017). Combining eDNA analysis with knowledge of plant phenology should provide insight on whether flowers or fruits are being eaten, and behavioural observation would determine whether feeding could result in pollination or seed dispersal.

Both ITS2 and rbcL are unable to detect and identify all plant taxa to species level. However, they are currently our best options for taxonomic identification in species-rich ecosystems, based on their reasonable taxonomic resolution and reasonably complete reference databases. Several prospective methods for eDNA-based identifications of species mixtures, including faecal and pollen samples, may become more feasible in the future. Methods that sequence the whole genome or a reduced-representation of the whole genome, or that skim organellar (chloroplast and mitochondrion) sequences from whole genome data have demonstrated improvement in taxonomic resolution over single DNA barcodes as well as quantification of the proportions of species in mixtures (Gómez-Rodríguez et al. 2015; Bista et al. 2018; Lang et al. 2019; Ji et al. 2020). This improved resolution could even allow identification at intraspecific levels, allowing analyses of pollen dispersal between populations. However, for diet analysis in species-rich ecosystems, particularly where the target taxa are highly mobile, the amount of DNA sequence required for reference databases could make such studies prohibitively expensive. There are currently several programs underway aimed at sequencing near-complete genomes of a large representation of the planet’s biodiversity (Lewin et al. 2018). As more genomic data become publicly available, and DNA sequencing technologies continue to improve, the use of whole genome methods for diet analysis will become more feasible, providing a powerful method for fine-scale taxonomic identification and quantification.