Introduced species can exert disproportionately negative effects on island ecosystems, but their potential role as food for native consumers is poorly studied. Telfair’s skinks are endemic omnivores living on Round Island, Mauritius, a globally significant site of biodiversity conservation. We aimed to determine the dietary diversity and key trophic interactions of Telfair’s skinks, whether introduced species are frequently consumed, and if diet composition changes seasonally between male and female skinks. We used DNA metabarcoding of skink faecal samples to identify animals (COI) and plants (ITS2) consumed by skinks. There were 389 dietary presence counts belonging to 77 dietary taxa found across the 73 Telfair’s skink faecal samples. Introduced taxa were cumulatively consumed more frequently than other categories, accounting for 49.4% of all detections, compared to cryptogenic (20.6%), native (20.6%), and endemic taxa (9.5%). The most frequently consumed introduced species was the ant, Pheidole megacephala, present in 40% of samples. Blue latan palm, Latania loddigesii, was the most frequently consumed endemic species, present in 33% of samples but was only detected in the dry season, when fruits are produced. We found a strong seasonal difference in diet composition explained by the presence of certain plant species solely or primarily in one season and a marked increase in the consumption of animal prey in the dry season. Male and female skinks consumed several taxa at different frequencies. These results present a valuable perspective on the role of introduced species in the trophic network of their invaded ecosystem. Both native and introduced species provide nutritional resources for skinks and this may have management implications in the context of species conservation and island restoration.

Skink sampling on Round Island

Faecal samples were collected in March, June, July and December 2015 (Figure 1). Skinks were caught opportunistically by noose or hand after which defecation was induced using a gentle abdominal massage. The faecal samples were placed in polythene bags and dried over silica gel. Skinks were released unharmed within ten minutes of capture at the locations where they were caught. Faecal samples were collected from 196 individual Telfair’s skinks (identified by their sex, size that were recorded and distinguishing markings and body deformations, which were photographed) on Round Island, and previously underwent DNA metabarcoding to identify the floral component of skink diet (Moorhouse-Gann 2018). Of the samples successfully sequenced, 82 were randomly selected for the current study from both dry (40) and wet (42) seasons.

Primer selection

Animal primers were tested in silico with a broad range of vertebrate and invertebrate taxa using PrimerMiner (Elbrecht and Leese 2017) and in vitro with DNA extracted from animals sampled on Round Island. BerenF-LuthienR (Cuff et al. 2020) provided the most comprehensive coverage, amplifying all Round Island invertebrate DNA extracts tested. UniPlant primer pair (Moorhouse-Gann et al. 2018), was used to amplify the ITS2 DNA barcode in plants and successfully amplify almost all plant species found on Round Island.

DNA extraction, PCR amplification, and sequencing

DNA extraction from Telfair’s skink faecal samples and preparation of plant DNA for 250 bp paired-end Illumina MiSeq high-throughput sequencing followed Moorhouse-Gann et al. (2018; Supplementary Information S2 and Table S1).

We used the following procedure to identify animal prey in the diet of Telfair’s skinks. Polymerase Chain Reactions (PCR) used 25 μL reaction volumes containing 5 μL DNA template, 12.5 μL of multiplex PCR mix (Qiagen, Manchester, UK), 2.5 μL of both forward and reverse primers (0.2 μM each), and 2.5 μL of nuclease-free water (Qiagen, Manchester, UK). PCR conditions went as follows: 95^oC for 15 minutes, 35 cycles of 95^oC for 30 s, 54^oC for 90 s, and 72^oC for 90 s, and 72^oC for 10 minutes, as instructed by the manufacturer (Qiagen, Manchester, UK). Each sample incorporated a unique combination of molecular identification (MID) tags (Binladen et al. 2007) that allowed for each skink to be identified after pooling and sequencing as per Brown et al. (2014). These 10-bp fragments were added to both the forward and reverse primers for each sample and thus dietary taxon sequences could be assigned to individuals. PCR products were then run through a 2% agarose gel stained with SYBR®Safe (ThermoFisher Scientific, Paisley, UK). Twelve negatives were included in each PCR run, 10 PCR negatives and two extraction negatives. Additionally, two positive controls consisting of a standardised DNA concentration (4 ng / μL) of known invertebrate species likely absent from the study site (Supplementary Information S1) were used to control for tag-jumping between samples in the filtering steps detailed below. PCR products were run in a Qiagen QIAxcel Advanced System (Qiagen, Manchester, UK) to measure relative DNA concentrations and later measured individually using a Qubit Fluorometer (ThermoFisher Scientific, Paisley, UK) for more accurate determination of DNA concentrations. Each sample was then pooled based on the relative DNA concentrations of the amplicon of interest as measured by the QIAxcel Advanced System. Negative controls were pooled based on the average volume pooled for the skink samples. The pooling process involved adding a volume from each sample as a proportion of the sample with the highest concentration of DNA, to ensure approximate equimolarity of DNA from each sample. Each pool was cleaned using SPRIselect beads (Beckman Coulter, Brea, USA), with a left-side size selection using a 1:1 ratio. After final elution, the pool was run on a Qubit Fluorometer, to measure DNA concentration (=49.6 ng / μL), as well as an Agilent 2200 TapeStation with D1000 ScreenTape (Agilent Technologies, Waldbronn) to check for significant levels of primer dimer, which were not found. This pool of MID-tagged samples was then used for library preparation using the NEXTflex™ Rapid DNA-Seq Kit following the manufacturer’s instructions (Bioo Scientific Corp, Austin, TX, United States), which is suitable for pools with DNA concentrations of 1 ng – 1 μg. A final DNA concentration was measured for the prepared library using a Qubit Fluorometer (=11.7 ng / μL) and was then sequenced on an Illumina MiSeq desktop sequencer (Illumina, San Diego, CA, United States) with a Nano cartridge using 2 x 250 bp paired reads (expected reads ≤ 1,000,000).

Bioinformatics

The Illumina Nano cartridge run generated 750,645 reads. High-throughput sequencing data for the animal component of Telfair’s skink diet followed the bioinformatic process of Drake et al. (2021): FastP (Chen et al. 2018) was used to check the quality of reads, discard poor quality reads (<Q30, <125bp long or too many unqualified bases, denoted by “N”), trim reads to a minimum length of 300 bp and merge read pairs from Miseq files (R1 and R2). Read pairs were assigned to samples and demultiplexed using Mothur v1.39.5 (Schloss et al. 2009), after which MID-tag and primer ends were removed. Unoise3 (Edgar 2010) was used to remove replicates, denoise the sequences, and group identical sequences into zero-radius operational taxonomic units (ZOTUs, which are clustered without % identity to avoid multiple species being nested within an OTU). Processed sequences were given taxonomic information from GenBank using BLASTn v2.7.1 (Camacho et al. 2009) with a 93% identity threshold. This threshold was chosen to capture the wide variety of invertebrates on Round Island to genus- or family-level, most of which have not been barcoded or formally described. When more than one taxon was assigned to a sequence, we manually checked the feasibility for the presence of each taxon on Round Island by searching published articles, unpublished reports, and personal observations of species accounts. If these manual checks were inconclusive, we assigned the sequence to a higher taxonomic level (genus, family, order, etc.). MEGAN Community Edition v6.18.9 (Huson et al. 2016) was used to analyse the BLAST output and assign taxonomic identities to each ZOTU. Using the lowest e-value (a value estimating the number of hits “expected” by chance when searching a database of a given size - in this instance anything less than 0.00001) the top hit was assigned to each sequence. Where top hits were taxonomic levels higher than species, these were manually checked and assigned to a feasible taxon or deleted from the analysis if erroneous. ZOTUs that were assigned to the same taxon were aggregated.

Data were cleaned for statistical analysis following the methods set out by Drake et al. (2021): the combined removal of the maximum read count in blanks and negative controls, and reads not meeting a pre-defined per sample threshold, removes both erroneous reads (laboratory contaminants and sequencing errors) that are likely to occur in low abundances mitigates tag-jumping and bleeding of over-represented taxa into other samples, whilst utilising a per sample threshold and those arising through tag-jumping and bleeding of over-represented taxa into other samples removes erroneous reads (laboratory contaminants and sequencing errors) that are likely to occur in low abundances. The maximum read count of known contaminants and other obviously erroneous ZOTUs across the dataset was calculated as a percentage of their respective total sample read count, and any read counts less than this were removed. For this, a threshold of 0.3% was applied, removing low-frequency laboratory contaminants and sequencing errors. Following this, the highest read count within a blank or negative per ZOTU was calculated and any ZOTU reads below this value were removed. In addition, we established an extra per-ZOTU filtering step, which removed remaining erroneous taxa. The per-ZOTU threshold was set to 0.74%. After these filters were applied, read counts were converted to presence-absence data for each sample. Nine samples were removed due to the absence of any dietary detections, leaving 73 samples to be taken forward for statistical analyses. Bioinformatic analysis for plant sequencing data followed Moorhouse-Gann et al. (2018) (Supplementary Information S2).

After animal ZOTUs were given taxonomic information, status of each taxon relative to Round Island was determined for each by manually searching for relevant data in published articles, unpublished reports, and personal species accounts, and then classified as “cryptogenic”, “endemic”, “introduced” or “native”. Cryptogenic species were defined as species that had no clear status, either because of poor taxonomic resolution, or because they may be known natives of the Indian Ocean islands but their history on Round Island is unknown. Plant status was taken from Moorhouse-Gann (2018).