Partitioning of diet between species and life history stages of sympatric and cryptic snappers (Lutjanidae) based on DNA metabarcoding

Lutjanus erythropterus and L. malabaricus are sympatric, sister taxa that are important to fisheries throughout the Indo-Pacific. Their juveniles are morphologically indistinguishable (i.e. cryptic). A DNA metabarcoding dietary study was undertaken to assess the diet composition and partitioning between the juvenile and adult life history stages of these two lutjanids. Major prey taxa were comprised of teleosts and crustaceans for all groups except adult L. erythropterus, which instead consumed soft bodied invertebrates (e.g. tunicates, comb jellies and medusae) as well as teleosts, with crustaceans being notably absent. Diet composition was significantly different among life history stages and species, which may be associated with niche habitat partitioning or differences in mouth morphology within adult life stages. This study provides the first evidence of diet partitioning between cryptic juveniles of overlapping lutjanid species, thus providing new insights into the ecological interactions, habitat associations, and the specialised adaptations required for the coexistence of closely related species. This study has improved our understanding of the differential contributions of the juvenile and adult diets of these sympatric species within food webs. The diet partitioning reported in this study was only revealed by the taxonomic resolution provided by the DNA metabarcoding approach and highlights the potential utility of this method to refine the dietary components of reef fishes more generally.

This file contains the following supplementary materials.
"Blocking primer development" Table S1. Fish16S assay and the host-specific blocking primer sequences (LEBP, Lutjanus erythropterus blocking primer; LMBP, L. malabaricus blocking primer).       Table S2. Taxonomic assignment resolutions (A) and the results of PERMANOVA (B), pairwise-PERMANOVA (C) using different percent identify match thresholds (0, 1 and 2%). Table S3. The list of prey taxa as lowest common ancestors (LCA) and number of ASV, species, genus and family which were assigned to the LCA. Table S4. The results of cross validation of diet compositions observations in the canonical analyses (CAP) ordination, following to the leave-one-out approach.

References
Biosystems). Different amounts of LMBP were added to each of the PCR mixtures to test their efficacy; 1 µL of 100 µM LMBP (10 times higher than Fish16S assay), 1 µL of 200 µM LMBP (20 times higher than Fish16S assay), and no blocking primer. Ultrapure water was added to the PCR mixture to bring the reaction volume up to 23 µL per sample. qPCRs were performed in 25 µL reaction volumes containing 2 µL of DNA extract and 23 µL of PCR mastermix using a StepOnePlus Real-Time PCR System (Applied Biosystems), using the following cycling program: (i) 95°C for 5 minutes, (ii) 45 amplification cycles of 95°C for 30 seconds, (iii) three different annealing temperatures (50, 54, and 58 °C) for 30 seconds, and (iv) 72°C for 45 seconds, and (v) a final extension step at 72°C for 10 minutes. Cycle threshold (Ct) values were recorded from the amplification curves and compared to assess the efficiency of the blocking. primer. Twoway ANOVA was used to test the significance of the two factors (the amounts of LMBP added and annealing temperature) on the average CT values using the software RStudio (v.1.0.143, https://rstudio.com/) 10 . The mean CT values of L. malabaricus DNA template qPCR were significantly higher when LMBP was applied in the PCR mixture (two-way ANOVA; F(2,17) = 173.55, p < 0.0001) whereas there was no difference between the LMBP to Fish16S ratio of 10 to 1 and 20 to 1 (Fig. S2). This indicates that LMBP effectively reduced the amplification of L. malabaricus DNA but adding extra LMBP did not improve the effectiveness. The annealing temperature had no significant effect on CT values (two-way ANOVA ; F(2,17) = 4.39, p = 0.05), and there was no significant interaction between the two factors (two-way ANOVA ; F(4,17) = 1.69, p = 0.24) (Fig. S2). Based on these results, downstream PCRs were carried out with the PCR reaction containing 10 times more blocking primer than Fish16S assay (1 µL of 100 µM blocking primer to 1 µL of 10 µM Fish16S assay in each PCR reaction), with an annealing temperature of 58°C.
The second in-vitro pilot study was carried out on the mock samples containing host-and prey-DNA mixed at different ratios. Firstly, qPCR was carried out with the DNA extracted from the fin clips of L. erythropterus, L. malabaricus, and S. sagax using Fish16S assay without blocking primers, and the relative concentration of their DNA templates were estimated using their CT values. DNA extracts were added to create the mock samples with different ratios of prey-to host-DNA (1:1, 1:100, 1:1000, and 10:1). The mock samples were amplified with and without blocking primers and sequenced each with unique six to eight base pair multiplex identifier (MID) tags in duplicate. Unique combinations of forward and reverse MID tags allowed us to assign sequences to a sample after metabarcoding the pooled samples. Singleend sequencing was performed for the pooled amplicons using an Illumina MiSeq platform in the Trace and Environmental DNA (TrEnD) Laboratory at Curtin University in Western Australia, following manufacture's protocols. Version 2 reagent kit and either Standard or Nano flow cell were used for 300 -500 cycles. Both LEBP and LMBP effectively inhibited the amplification of host-DNA and increased the chance of detecting prey DNA when the prey-host DNA ratios were 1 to 1, 1 to 100, and 10 to 1 (Fig. S3). When the L. erythropterus DNA was 1000 times higher than the S. sagax DNA, LEBP was no longer able to suppress the amplification of host-DNA in order to detect the prey DNA (Fig. S3).
Following these results, downstream PCRs were carried out with the PCR reaction containing 10 times more blocking primer than Fish16S assay (1 µL of 100 µM blocking primer and 1 µL of 10 µM Fish16S assay in each PCR reaction), with the annealing temperature of 58°C. Table S1. Fish16S assay and the host-specific blocking primer sequences (LEBP, Lutjanus erythropterus blocking primer; LMBP, L. malabaricus blocking primer). Fish16S target teleost sequences at 16S rRNA mitochondrial gene region. The first 10 base pairs of the blocking primers (bold and italic font) overlap with the 3' end of the forward Fish16S assay, followed by 15 base pairs of the host specific sequences. The C3 spacer at the 3' end is a modified DNA oligonucleotide, which inhibits annealing.

Assessment of sampling and sequencing depth
The relationship between number of reads and prey taxa was examined using Pearson correlation in order to inspect whether the number of prey taxa detected was affected by sequencing depth. Analyses of variance (ANOVA) was also carried out to test the differences in mean number of reads and prey taxa detected between each group. The assumption of the samples coming from a normal distribution was tested and the number of prey taxa were transformed using a square-root transformation to meet the assumption. Rarefaction curves were plotted describing the diet of each species and life history stage. Random subsampling of sequences was conducted 1000 times at every 1000 reads for each sample, following to the approach explained by Colwell 11 , and the total number of ASV and prey taxa detected by subsampling were averaged within the samples of the same group (species and life history stage). Species accumulation curves were plotted to assess whether the number of samples were sufficient to capture the majority of their potential prey taxa consumed by each species and life history stage. The software RStudio (v.1.0.143, https://rstudio.com/) was used to carry out statistical analyses and subsampling, and produce plots 10 .  Table S2. Taxonomic assignment (A) as well as the results of PERMANOVA (B) and pairwise-PERMANOVA (C) using different percent identify match thresholds (0, 1 and 2%). The threshold defines the maximum difference between the percent identity matches of primary and nonprimary reference sequences allowed in the lowest common ancestor (LCA) assignment algorithm. When the difference between the percent identity matches of primary and nonprimary reference sequences was more than the threshold, the non-primary reference sequences were omitted prior to LCA assignment. PERMANOVA (B) and pairwise-PERMANOVA (C) examined the differences in diet composition of juvenile and adult Lutjanus erythropterus (LE) and L. malabaricus (LM), using 9999 permutations. Fullness of stomach was incorporated into the analyses as a covariate in order to test its effect on diet composition. The tests were based on Jaccard coefficient matric for presence and absence (PA) datasets.
A. Taxonomic Table S3. The list of prey taxa as lowest common ancestors (LCA) and number of ASV, species, genus and family which were assigned to the LCA. % match indicate the range of % similarity between each ASV and the reference sequences of assigned taxa. Where LCA taxa level is genus or higher, the species list contains more than one species with a common ancestor.