Elucidating the diet of the island flying fox (Pteropus hypomelanus) in Peninsular Malaysia through Illumina Next-Generation Sequencing

Study site

We conducted this study on Tioman Island (2°48′38″N, 104°10′38″E; 136 km²; Fig. 1A), located 32 km off the east coast of Peninsular Malaysia in the State of Pahang. This research was approved by the Economic Planning Unit of Malaysia (Permit number: 3242). Much of the island inland is still covered by primary tropical rainforest, which has been designated as Pulau Tioman Wildlife Reserve (83 km²). It has a hilly topography, with flat areas only along the coast (Abdul, 1999). The area designated as a wildlife reserve is composed of lowland mixed dipterocarp forest and hill dipterocarp forest. Most forested areas are still inaccessible due to the rugged topography, with many steep slopes and rocky outcrops (Latiff et al., 1999). The climate is tropical, uniformly warm and humid throughout the year (Hasan Basyri et al., 2001), but the island experiences the northeast monsoon from November to March (Bullock & Medway, 1966).

There are currently seven villages on the island, situated along the coastline (Fig. 1A). The majority of the local people are Muslim, and therefore due to religious dietary restrictions do not hunt bats for food or medicine (Aziz, 2016). As the island’s marine area is also a designated Marine Park and a popular tourist destination, many of the local people are heavily involved in the tourism industry (Abdul, 1999).

Study species

The island flying fox (Pteropus hypomelanus), also known as the variable flying fox and the small flying fox, roosts gregariously, forming colonies of up to 5,000 individuals. It is a widespread insular species, considered to be abundant throughout a distribution range that extends from the Maldives and Indian islands in the west to Melanesia in the east. Because of this, it is considered to be Least Concern on a global scale by the IUCN Red List; however, populations are decreasing (Francis et al., 2008; Olival, 2008), and the species is now listed as endangered on the Malaysian Red list (DWNP , 2010).

On Tioman, the island flying fox can be found roosting permanently in two villages: Tekek, on the west coast, and Juara, on the east coast (Fig. 1A), and forages throughout the island (Medway, 1966; Ong, 2000). Monthly roost counts conducted during March–October 2015 using a thermalscope (Pulsar Quantum HD38S) yielded estimated ranges of 2,178–5,385 individuals for the entire island (see Fig. S1 for further details). Local people have reported that the flying foxes do forage in other villages on the island.

Flying fox dropping collection

Collection of flying fox droppings took place once a month during March-October 2016 (i.e., eight months). Samples of droppings consisting of faeces and ejecta were collected for three mornings in the last week of each month from three separate day roosts in Juara (east coast) and two separate day roosts in Tekek (west coast). The number of roosts and sampling days were determined based on species accumulation curves of pollen morphospecies that were detected through preliminary microhistological analysis in June 2014. Program EstimateS (version 9.1.0; http://viceroy.eeb.uconn.edu/estimates/) indicated that sampling completeness (i.e., observed/estimated number of species; Soberón, Llorente & Onate, 2000) was around 97% using this sampling regime.

In Juara, three suitable roost trees (Fig. 1B) for sampling were selected based on accessibility and also on the highest/largest amount of faecal/ejecta splatter produced under the roost, in order to maximise sample yield. As flying foxes often shifted roosts or temporarily abandoned degraded roosts, this meant that sometimes different roosts were sampled in each location every month or even every morning, although most roosts were consistently sampled each month due to their constant high occupancy and best accessibility.

Plastic sheets measuring 0.8 × 1.0 m were placed under each roost after dark, once the bats had exited the roost to forage. The roosts were then visited the next morning for collection starting at 0700 h and ending at 1200 h (bats typically returned to the roosts around 0500–0600h); the plastic sheets were carefully moved away from the roost to a clear area for processing (Fig. 1C). As it was often difficult to differentiate faeces from ejecta, both were collected and analysed equally as ‘droppings’ (Fig. 1D). Droppings collected for processing were selected based on unique colour and texture, as this was assumed to be representative of plant diversity in the bats’ diet. Following the approach used by Stier & Mildenstein (2005) based on short gut-passage time for flying foxes (12–34 min; Tedman & Hall, 1985), we assumed that each bat voided its last meal once, and therefore each dropping represented a different individual’s food choice. We devised our own novel collection protocol where droppings were collected by swabbing them with a cotton bud, then placing each individual dropping into a 5 ml Eppendorf tube containing ∼1,000 μl of 95% ethanol. These tubes were then kept cool in the field, either by storing in a conventional freezer or by using a portable cooler box with ice packs, for 1–3 days before being transported off the island and then stored in a −80 °C freezer.

In order to simultaneously test the utility of NGS and compare it with conventional approaches, we collected two duplicate sets of 10 individual droppings from one single roost in Juara village during a single morning on 6 May 2015. One sample set was then kept in a conventional fridge for microscope analysis, whilst the other set was stored in the −80 °C freezer for molecular analysis.

Site-specific plant reference collection

We first checked a published list of genera of known food plants for Pteropus across its range (Marshall, 1985), cross-checked this against a preliminary checklist of seed plants for Tioman (Latiff et al., 1999), and also obtained information on possible flying fox food plants through talking to local people in Juara. We then searched for genera of similar plants in and around the two villages with the aid of a local plant expert. The botanical identification of plants (at least to genus) were subsequently verified by a trained botanist familiar with plants from the region. When we found an individual plant from one of these genera, we recorded its GPS location and collected pollen, fruit, and/or seeds if it was flowering or fruiting. We also collected 3–5 mature leaves for DNA extraction. The leaves were stored in Ziploc bags with silica gel under cool conditions to retard decomposition rates. Leaf samples from 19 different plant species were obtained for this purpose, constituting a preliminary library (Table S1).

Genomic DNA was extracted from approximately 25 mg of one leaf from each plant species using DNAeasy Plant Mini Kit (Qiagen, Halden, Germany) according to the manufacturer’s protocols. Primers targeted the rbcL gene, a protein-coding gene associated with the chloroplast genome of all living plants. DNA amplifications were performed in a mastermix containing 1 µL of DNA, 25 µL of OneTaq Quick-Load 2X Master Mix with Standard Buffer, (New England Biolab, Ipswich, MA, USA), 1 µL of 10mM forward primer rbcLaf-M13, 1 µL of 10 mM reverse primer rbcLa-revM13 (Table S2), and 22 µL of nuclease-free water. The PCR protocol was started with an initial denaturation step for 30 s at 94 °C, followed by 30 cycles of 30 s at 94 °C, 30 s at 48 °C, 40 s at 68 °C, and final elongation for 2 min at 68 °C. The PCR products were purified using 0.8× volume ratio of Agencourt Ampure XP beads (Beckman Coulter, Inc). The purified samples were sent to 1st BASE laboratories (http://www.base-asia.com/) for Sanger sequencing. The sequencing results were quality trimmed using CodonCode TraceViewer (http://www.codoncode.com/TraceViewer/) and aligned using MAFFT version 7.0 (Katoh & Standley, 2013).

A phylogenetic tree (Fig. S4) was constructed to support morphology-based identification of site-specific reference plants. Additional rbcL gene sequences that displayed high similarity to the sequenced plant specimens were mined from the BOLD/NCBI databases, aligned with MAFFT version 7.0 and trimmed using TrimAl version 1.9 optimized for maximum likelihood tree construction (Capella-Gutiérrez, Silla-Martínez & Gabaldón, 2009). The trimmed alignment was subsequently used to construct a maximum likelihood tree with FastTree using the –nt (nucleotide) and –gtr (generalised time-reversible model) setting.

Laboratory procedure

Based on the alignment of rbcL gene fragments from the site-specific reference plant DNA (Fig. S2), we designed a primer pair targeting 198 bp of the rbcL gene (minus primer sequence) using Primer3 (http://bioinfo.ut.ee/primer3-0.4.0/) on default settings (Table S2). Prior to primer synthesis, partial Illumina adapter sequences were added to the 5′ end of the designed primers, rbcL-357F and rbcL-556R, to allow barcoding and sequencing on the Illumina platform.

Individual droppings were pooled according to roost (n = 5, 2 in Tekek and 3 in Juara) and month (n = 8), creating 40 separate mixtures for analysis. The tubes containing the daily samples were first vortexed for 2 min to homogenise the content and subsequently, 1,000 µL of the sample was pipetted into another tube to form the mixture. Next, 100 µL of the mixture was used for DNA extraction similarly using DNAeasy Plant Mini Kit (Qiagen, Halden, Germany) instead of a stool-specific DNA extraction kit to improve the recovery of plant-derived DNA from faecal samples.

PCR reaction was performed using IlluM_rbcLF and IlluM_rbcLR. The 20 µL PCR cocktail consists of 10 µL Q5 Hot Start High-Fidelity 2× Master Mix (New England Biolab, Ipswich, MA, USA), 1 µL each of 10 µM forward and reverse primer, 1 µL gDNA and 7 µL nuclease-free water. All reactions were performed in a Veriti^® 96-Well Fast Thermal Cycler with the following protocol: initial denaturation for 30 s at 98 °C, 25 cycles of 10 s at 98 °C, 30 s at 55 °C and 10 s at 65 °C, with a final 1 min extension at 65 °C. The PCR product was purified using 0.8× vol. ratio Agencourt Ampure XP beads (Beckman Coulter, Inc). Then, 1 uL of Index 1 and Index 2 primers from Nextera XT kit were added to 3 uL of purified PCR product and combined with 5 uL of Q5 Hot Start High-Fidelity 2X M aster Mix (New England Biolabs, Ipswich, MA). The PCR protocol was as followed: initial denaturation for 30 s at 98 °C, 8 cycles of 10 s at 98 °C and 1 min at 65 °C, with a final 1 min extension at 65 °C.

Sequence analysis and taxon assignation

The purified amplicons containing the full length Illumina adapter and appropriate unique barcode were then quantified using KAPA Library Quantification kit (Kapa Biosystems, CapeTown, South Africa) on the EcoRealTime PCR system (Illumina, San Diego, CA, USA). Based on the qPCR data, the amplicons were normalised, pooled and subsequently sequenced on the MiSeq (2 × 250 bp paired-end run) located at the Monash University Malaysia Genomics Facility.

Illumina Nextera adapters and primer sequences of the reads were trimmed off using Trimmomatic v 0.33 and FastX trimmer, respectively (Bolger, Lohse & Usadel, 2014; http://hannonlab.cshl.edu/fastx_toolkit/). The trimmed paired-end reads were then merged using PEAR (Zhang et al., 2014) using default settings. Dereplication, singleton removal and Operationally Taxonomic Unit (OTU) clustering (default setting of 97% identity clustering; -id 0.97) were performed using the pipeline implemented in UPARSE (Edgar, 2013). The filtered OTUs were translated into protein sequence and manually inspected for translated sequence containing stop codon(s), which were removed from subsequent analyses. Taxonomic assignment of OTUs was conducted by matching sequences against those from the BOLD/NCBI non-redundant nucleotide database as of December 2016. The BLASTn output for each OTU was downloaded in “XML” format and imported into MEGAN6 (Huson et al., 2016) to calculate their Lowest Common Ancestor (LCA). LCA predicts the taxonomic rank of OTUs based on matches to a set of reference lineages such that the predicted taxonomic rank is shared by all matches within user pre-defined parameters. The implemented parameters for MEGAN6 were: Minimum Blast bit score = 200, Max expected E-value = 0.01, Min Percent Identity = 90.0, Top percent (the percentage of the best score a match needs to lie within) = 10 and Weighted LCA% = 75. To estimate the OTU relative abundance in each sample, reads were mapped to the filtered OTUs via USEARCH (97% similarity cut-off) and normalised to 10,000 reads (Caporaso et al., 2010).

Microhistological analysis

For the 10 dropping samples collected in May, we sent one set for NGS analysis (following the protocol above) and used another set for microscope analysis. For the latter, we first manually broke up the dropping contents in the tube to produce a relatively more representative liquid sample. We then droppped 1–3 drops of this liquid onto a microscope slide using a pipette. Fuchsin jelly was added to this in order to stain pollen grains within the dropping, a slip cover was placed on top, and the jelly was then melted over an open flame, sealing the slip cover to the slide. The slide was then cooled down in a conventional fridge in order to allow the jelly to solidify again before examination.

Once cooled, we placed the slide under a conventional light microscope (Leica DM E) and first examined it using 10/0.25 magnification in order to detect pollen grains and other plant parts. We used a self-made reference collection as well as photos from Start (1974), S. Bumrungsri (http://www.seabcru.org/seabcru-resources) and Mohamed (2014) to identify pollen and other plant parts. When necessary, we used higher magnification (40/0.65) to view the pollen grains. We used just ‘presence/absence’ to assess pollen to avoid quantification biases towards species that naturally produce greater amounts of pollen. Following Thomas (2009), we considered a species present if we found three or more pollen grains on a single slide.

Site-specific plant reference database

For the 19 specimens in the site-specific plant reference database, the phylogenetic tree (Fig. S4) based on a 689 bp aligned rbcL gene region recovered the monophyly of several plant genera with high SH-like local support (>0.90). However, the sister grouping among genera in general was not well-supported, possibly due to the lack of comprehensive plant taxon sampling, and the use of a single genetic marker. The rbcL-based phylogeny generally supported the morphology-based identification of reference specimens of potential flying fox food plants, with 15 specimens correctly identified at least to genus level as shown by their monophyletic clustering with plant species of the same genus (Table 1). Sample PTMN02 (identified as Streblus asper based on morphology) was an interesting exception as it formed a sister group with the Ficus clade (SH-like local support = 0.938), but not with other congeners. rbcL sequences representing the plant species collected on Tioman were scarce in the BOLD/NCBI databases, suggesting these databases have insufficient sequence representation that can affect taxonomic assignment of plant DNA in the droppings.

Feasibility of using NGS to study flying fox diet

With our newly designed rbcL primer, we were able to successfully extract, amplify, and subsequently identify plant DNA from all of the collected flying fox droppings. Initially, a total of 160 Operationally Taxonomic Units (OTUs) were recovered from the sequencing reads, of which 29 OTUs (Table 1) were retained after filtering based on cumulative relative abundance (>0.5%) and presence of stop codon(s) in reading frame. Using a conservative LCA approach, we identified at least three different plant genera and at least 14 plant families from the droppings (Table 1). In addition, 8 OTUs matched with specimens from the site-specific plant reference collection.

Based on sampling completeness (calculated using EstimateS) for OTU relative abundance data from five roosts (data pooled over three days) per month using Chao 1 species richness estimator (good for datasets skewed towards low abundance classes; Chao, 1984), sampling completeness was relatively high for the months March (99%), April (100%), June (100%), August (100%), September (88%), and October (96%). However, sampling completeness could be improved for May (55%) and July (79%). The month of May, which had the lowest overall sampling completeness, also had the highest number of droppings collected (Table S3).

The results from our NGS analysis of island flying fox droppings over eight months suggest that the diet at both Juara and Tekek during this time was dominated by four different plant taxa (Table 1; Fig. 2): OTU 1 (Ficus sp.), OTU 3 (likely to be Mangifera indica based on sequence match with site-specific plant reference collection) and OTUs 4 & 5 (Rubiaceae). Spatio-temporal patterns in the relative abundance of these four taxa in the diet were observed during the sampling period (Fig. 2). For example, OTU 5 appeared to be consumed in similar proportions at both Juara and Tekek across all months whereas OTU 4 was consistently consumed in low proportions in Tekek yet consumed irregularly in Juara over the same period (Fig. 3). Even between different roosts in the same site, spatio-temporal differences were observed, such as for OTU 7 (Fig. S5).

Performance of NGS vs. microhistological analysis

Microscope analysis identified two plant taxa in flying fox droppings (Table 2). Out of 10 individual droppings, three contained durian (Durio sp.) pollen. All the other droppings contained fig parts exclusively; no other plant parts were detected. Durian pollen occurred at extremely low abundance; in all cases, only 3–4 grains were detected per slide. No other pollen or plant parts were detected. NGS identified the same two plant taxa detected by microhistological analysis, and further identified an additional six plant taxa. However, Durio was not detected in the same samples as those identified via microscope.

NGS is a reliable tool to study flying fox diet

We have demonstrated that identification of plant taxa to family level is generally possible based on the partial sequence of rbcL using the LCA approach. In addition, some OTUs in our study were successfully assigned to the genus level. In order to be conservative, however, we avoided assigning most OTUs to species level, unless there were matches with BOLD/NCBI database sequences and site-specific reference plant sequences. Asspecies-level plant identification based solely on NGS is not straightforward, the use of a site-specific food plant DNA reference database is vital to help identify plant species in flying fox diets.

Other genes have been successfully used to identify plant species in animal diets. For example, Valentini et al. (2009a) found the trnL intronic region to be effective for Asian mammals, birds, and invertebrates, identifying 50% of the plant taxa found in the diets of these animals to species level. The same approach has been used for European bison (Kowalczyk et al., 2011), alpine chamois (Raye’ et al., 2011), and red-headed wood pigeons (Ando et al., 2013). The trnL intron has been previously reported to evolve more than three times faster than the protein-coding rbcL gene, thus potentially harbouring more variation and phylogenetic signal per base pair (Gielly & Taberlet, 1994). However, we chose the rbcL gene instead of the P6 loop of the trnL intron to study flying fox diet, as the rbcL gene is currently one of the two genes (the other being matK) that is increasingly being used for plant species identification. This is useful for studying plants originating from a less-studied region of high biodiversity such as Southeast Asia. Unlike the trnL database, the rbcL database is consistently growing in the BOLD database, and as a result is also likely to be represented by more sequences from properly identified and vouchered specimens. In addition, given that rbcL is a protein-coding gene, it enables the screening of erroneous OTU resulting from sequencing and/or amplification errors based on the presence of stop codon(s) in the translated reading frame. Instead of using recently developed and robustly tested primers (e.g., Little, 2014), we designed a new set of rbcL minibarcode primers due to the current lack of rbcL sequence representation for plant species from Malaysia or specifically, Tioman. One caveat of this approach would be the possibility of preferential primers binding to known diet items instead of unexpected ones. However, our results seem to indicate that our newly designed primers were, to a certain extent, capable of recovering OTUs belonging to a wide variety of plant families, in part due to the diverse representation of plant reference specimens that contributed to the expanded taxon coverage of our primers. Future studies aiming to achieve greater power of identification of plants in flying fox diets could consider using more than one target region (e.g., Hibert et al., 2013; Clare, 2014) coupled with an additional in silico PCR optimisation step using ecoprimers to improve the reliability and universality of the newly designed primers (Riaz et al., 2011). Also, in order to completely eliminate bias associated with PCR, metagenomic shotgun sequencing could be performed, albeit at a relatively higher cost depending on the required sequencing coverage. (e.g., Srivathsan et al., 2015; Srivathsan et al., 2016).

NGS can complement microhistological analysis

Our results showed that NGS can provide greater insights into the diet of flying foxes than conventional microhistological approaches by detecting a wider range of plant taxa, thus highlighting the utility and discriminatory potential of the newly designed rbcL primers to study flying fox diets. More importantly, the use of NGS allowed us to identify plant species even when no physical plant parts were found in the flying fox droppings. The plant genera and families detected from NGS have also been recorded by botanists as being present on Tioman, including the top four genera/families detected most abundantly in the droppings (Latiff et al., 1999; Mohd. Norfaizal et al., 2014).

In our study, attempts to use microscope analysis to identify plant parts in droppings proved to be challenging, as no pre-existing reference collection was available. Building our own microhistological reference collection for Tioman was time-consuming and labour-intensive—and the resulting collection often did not match up with the plant parts found in the flying fox droppings. However, obtaining DNA from plant specimens is still necessary to narrow down the identity of OTUs to species level. Indeed, 8 out of 29 OTUs had 100% matches to the sequences of plant specimens collected from the study site, highlighting the importance of building a comprehensive local sequence library beforehand, preferably specific to one’s particular study site.

It is important to note that NGS did not detect Durio in the same individual droppings as those identified via microscope. This is likely due to the low abundance of this plant taxon in the droppings affecting detection probability, especially since the NGS analysis used a more general primer that was not specific to Durio. This pollen detection probability is another caveat to be aware of; Scanlon et al. (2014) have cautioned that faecal subsampling methods can potentially lead to inaccurate detection of pollen in dietary studies, regardless of which method is used.

Caveats

Our sample collection method in the field, selecting only for droppings with unique colour and texture, may have introduced a bias that could result in underestimating the relative abundance of OTUs in the droppings. In particular, sampling completeness for the months May, June and July were relatively low, showing that more roosts and/or days needed to be sampled in order to obtain a complete representation of diet for these months. However, roost count data (Fig. S1) show that this high diversity in diet was not influenced by population abundance. Instead, perhaps diet choice and/or food resource diversity were relatively higher during this period compared to the rest of the year. Future studies should aim to collect every single dropping found underneath a roost to improve representativeness.

Given the potentially short flying fox gut passage times (Tedman & Hall, 1985), droppings collected from day roosts in the morning may also bias the analysis results towards food items that were consumed only at the end of the foraging period (Schmelitschek, French & Parry-Jones, 2009). Food plants that were consumed during the start or middle of the evening may not have been detected by our methods. Although Banack & Grant (2002) observed flying foxes returning to food resources that were foraged upon earlier, before then returning to day roosts, this is still a potential caveat to bear in mind. For example, primates are known to exhibit temporal patterning in diet choice, structuring their diet throughout their foraging period with different food items; it is believed that this is due to how different foods are processed, and give energy, at different rates, and therefore helps to ensure that the animals maintain high energy levels (Robinson, 1984; Ganzhorn & Wright, 1994; Chapman & Chapman, 1991). Given the size of Tioman, and the logistical challenges of observing flying foxes foraging, the best way to overcome this possible information gap is to conduct GPS tracking studies.

We also acknowledge that NGS approaches to diet identification are semi-quantitative because chloroplast abundance is variable in different plant species and different parts of the leaf. Ultimately, the ability of NGS to accurately identify food plants will always depend on sequence specificity of the primers. While the NGS approach has proven to be useful in elucidating the island flying fox’s varied diet on Tioman, for animals with such a diverse phytophagous diet, primer specificity will always be a limiting factor and there is a chance that unknown plant species will not be detected due to primer mispriming. Also, identical chloroplast DNA sequences can be present in different but related species, making it impossible to distinguish closely related plant species from each other in the diet. This could be one possible factor why several OTUs could only be identified to order/family/subfamily levels, suggesting that they require further phylogenetic investigation and/or may benefit from identification based on more rapidly-evolving plastid-coding genes such as matK. It is worth noting that the sequenced rbcL gene of some plant specimens collected in this study did not exhibit 100% identity matches to species in the BOLD/NCBI databases, which may be attributed to genetic diversity at the intra-species level or gaps in the database i.e., certain plant species consumed by the flying foxes may not yet have their corresponding sequences deposited in the database. Nevertheless, there is an urgent need for the BOLD/NCBI databases to have more representation of plant sequences from Peninsular Malaysia and Southeast Asia in general.

Another limitation of the NGS approach for generalist diets is that it does not identify which part of the plant was consumed. For specialised frugivores, nectarivores, or herbivores (e.g., folivorous), this may not be an issue. Flying foxes, however, are generalists which consume fruits, flowers, nectar, and even leaves (Marshall, 1985). It is this dietary plasticity which allows them to perform more than one ecological role in tropical landscapes. Therefore, identifying which plant parts are actually consumed is a crucial step towards identifying the ecosystem services that these bats provide. NGS can provide a first step towards identifying flying fox diet but should not be viewed as a replacement for microhistological analysis. Nevertheless, this approach has shed new light on flying fox diet by discovering plant taxa that were entirely missed out by the conventional approach. Ideally, studies using NGS should be combined with micro-histological analysis in order to fill in the gaps and broaden our understanding of pteropodid diet and foraging ecology. NGS can also be used in combination with comprehensive and long-term data on plant phenology, to observe which food resources are available at which time. Following on from this preliminary study, the identification of specific food plants via NGS can now help guide more in-depth plant sample collection and phenological observations.