We conducted the study at the Ramat Hanadiv Nature Park, located on the southern edge of the Mount Carmel ridge, Israel. The landscape is covered by a Mediterranean woodland (maquis), with steep, rocky slopes and patches of shallow soil. The climate is Mediterranean, characterized by hot, dry summers and mild, rainy winters with a mean annual rainfall of 600 mm (which falls between October and April). The growing season is closely associated with rainfall distribution and green annual herbaceous vegetation covers the soil patches from January to May. The vegetation in the park is dominated by low trees (mainly Phillyrea latifolia L.) and tall shrubs (Pistacia lentiscus L. and Calicotome villosa Poir. Link). Bushes of Ephedra foemina Forskk., Asparagus stipularis Forskk. and Sarcopoterium spinosum L. Spach grow between the trees. Around 500 herbaceous species can be found in the park, mostly Cyclamen persicum, Brachypodium distachyon, Urospermum picroides, Asphodelus ramosus, Anagalis arvensis, and Convolvulus pentapetaloides (Hadar et al., 1999).

Large areas of the park are grazed by a resident herd of goats containing 150 adult females (Alpine, Damascus and Mamber breeds). The goats forage the park daily, for approximately 4 h, and receive supplementary feed (hay and processed pellets including 16% crude protein, vitamins and minerals) upon their return to the pen.

To trace individual and seasonal variation in ingestion of PDA, we collected fecal samples from a group of 18 marked individual goats, once a month from February to June (end of winter to the beginning of summer) when arthropod abundance is at its highest. We collected the feces (~5 g) directly from the goats’ rectum using sanitary gloves (replaced for each goat to avoid contamination) upon their return from the pasture. The samples were placed in plastic zip-lock bags and kept in a cooler box with ice packs. Upon return to the lab (within 3 h), we pulverized each sample using a mortar and pestle, thoroughly mixed (homogenized) it and transferred it to-80° C until DNA extraction.

DNA extraction. We extracted DNA from the fecal samples using a QIAamp Fast DNA Stool Mini Kit (ID: 51604, Qiagen, Germany) according to the manual of “isolation of DNA from stool for pathogen detection,” with the following modifications: we conducted the lysis with InhibitEX buffer at 95° C and conducted the incubation with proteinase K and AL buffer at 70°C for 20 min. We eluted the DNA twice with 50 μl of AE buffer (total 100 μl) and stored it at-20°C.

PCR amplification. We PCR-amplified DNA using the Zeale et al. (2011) primer set (ZBJ-ArtF1c: 5’-AGATATTGGAACWTTATATTTTATTTTTGG-3′ and ZBJ-ArtR2c: 5’-WACTAATCAATTWCCAAATCCTCC-3′) which targets a 157 bp segment of the mitochondrial cytochrome c oxidase subunit I (COI) gene in arthropods. These COI markers were chosen, despite possible taxonomic biases (Clarke et al., 2014; Deagle et al., 2014; Elbrecht and Leese, 2017) due to their comprehensive reference database (Ratnasingham and Hebert, 2013) and ability to effectively discriminate at the species level (Deagle et al., 2014). Furthermore, they have been successfully applied in many dietary studies involving arthropods (Kaunisto et al., 2017; Alberdi et al., 2018; Vesterinen et al., 2018; Rytkönen et al., 2019).

For PCR amplification we used 12.5 μl EmeraldAmp MAX HS PCR Master Mix (Takara bio Inc., Otsu, Shiga, Japan), 1 μl of forward and reverse primers (0.5 ng/μl), 5 μl of DNA template (10–100 ng), 0.2 μl of BSA (0.1 μg/μl) and 5.3 μl of sterile water for a total reaction volume of 25 μl per sample. The PCR program included an initial denaturation step at 98°C for 2 min, followed by 35 cycles of 10 s at 90°C; 30 s at 52°C; and 30 s at 72°C, and a final elongation step of 1 min at 72° C. Negative controls were carried out for each PCR assay. We verified the amplification by agarose gel electrophoresis and stored the PCR products at-20° C.

Library construction. COI amplicon sequencing was performed on an Illumina MiniSeq instrument (Illumina, CA, USA) at the Sequencing Core (UICSQC), University of Illinois, Chicago, using a two-stage amplification protocol as described by Naqib et al. (2018). PhiX DNA was used as a spike-in control. A total of 79 fecal samples were sequenced based on gel verification (Supplementary Table S1): 16 samples from February, 15 from March, 17 from April, 18 from May and 13 from June. Sequencing results (pair-ends, 2 × 150 bp) were delivered as FASTQ files.

Sequence Analysis. We processed the Illumina FASTQ files in R using DADA2 (Callahan et al., 2016) to clean them of PCR and Illumina-sequencing errors. We trimmed and filtered the FASTQ reads for quality using the command “filterAndTrim” with the parameters maxN = 0, maxEE = c(5,5), trimLeft = c(30,15). Error rate estimation was carried out using the “learnerror” command with default parameters (randomize parameter was set to “TRUE” to randomly sample nucleotides and reads across all samples). We then implemented the DADA2 algorithm for error correction and merged forward and reverse reads using the “mergePairs” command (minimum overlap of 6 bp). We removed likely chimeras using the command “removeBimeraDenovo” with default parameters. This process generated an amplicon sequence variants (ASV) count table containing 1,313,939 high-quality reads. ASVs with sequences longer or shorter than the expected barcode length (>170 bp or < 150 bp) were eliminated, retaining 409,596 reads binned in 666 ASVs.

Taxonomic identification. We aligned the ASV sequences against the NCBI GenBank nt-database and imported the BLAST output files into “MEGAN community edition” [v.6.21.11 (Huson et al., 2016)] for taxonomic classification. This was done using the Lowest Common Ancestor (LCA) approach with parameters set to: bitscore >100, E value <10^7, min support =1, top percent =10, LCA algorithm = naïve and percent to cover = 100. Any unclassified sequences were examined against the BOLD database (Ratnasingham and Hebert, 2007). We retained only sequences belonging to insects and arachnids. Moreover, we removed any taxa found in a single sample alone (listed in Supplementary Table S2), to provide a more accurate estimate of the arthropod community. This resulted in 76 samples with 342,435 reads and 223 ASVs (Supplementary Table S3). The raw sequence data is available at the NCBI database under BioProject accession number PRJNA579572.

To avoid sequencing depth bias, we rarefied the data to 913 reads per sample using the VEGAN R package command “rarefy” (Oksanen et al., 2007), as the curves of most samples reached a plateau at this point (Supplementary Data Sheet 2 Supplementary Figure S1). This process resulted in 57 goat fecal samples (9 samples from February, 10 from March, 13 from April, 13 from May and 12 from June), binned in 215 ASVs (Supplementary Table S4).

Occurrence data (presence/absence) and relative read abundance (RRA) are two approaches currently used to analyze sequence data in DNA metabarcoding dietary studies. Occurrence data is considered to be more cautious, yet RRA may be more accurate at representing population-level diets (Deagle et al., 2018). Indeed, some recent DNA metabarcoding studies of insectivore diets focus on RRA analysis (Vesterinen et al., 2018; Pringle et al., 2019). For these reasons, we chose to use RRA as the main method for analyzing our data (calculated as the proportion of reads of each ASV in each fecal sample). To avoid potential biases of the RRA analysis (which can arise from differential digestibility, amplification biases, etc.), we also performed a supporting analysis of percent of occurrence (POO, calculated as the number of samples containing a given food item, rescaled to 100% across all food items; see Deagle et al., 2018; Supplementary Table S5). In both analysis we removed any sequences with an abundance of <1% (i.e., true presence threshold, see Deagle et al., 2018; Supplementary Table S6). No assumptions were made regarding the actual proportions or biomass of arthropods.

We calculated diversity (Shannon H′) and richness (Fisher’s alpha) indices of the ASVs using PAST (Hammer et al., 2001). Since not all individual goats had a valid sample for each month, we selected goats that had at least three consecutive samples, resulting in six goats with samples from April–May (hereafter “six-goat group”). This group was used to compare the diversity indices among sampling months using a Friedman’s test (post hoc Wilcoxon signed-rank test, Bonferroni corrected for multiple comparisons). The number of arthropod families ingested among sampling months and among individual goats was also compared within the six-goat group using a repeated measures ANOVA (with the individual goats and sampling months as the repeated-measures factors). Both tests were performed using IBM SPSS software v.25.

We examined the similarity of the overall arthropod community (ASVs) among February–June and among all individual goats using a non-metric multidimensional scaling analysis (NMDS) in R (RRA: Bray–Curtis similarity matrix; POO: Jaccard similarity matrix). To test whether the community significantly differed among sampling months and individual goats (and the six-goat group independently) we performed a permutational multivariate analysis of variance (PERMANOVA) with 999 permutations using the R-VEGAN function ADONIS (RRA: Bray-Curtis; POO: Jaccard) for each factor separately. When considering all individual goats, goat ID was set as strata. In the six-goat group, a nested ADONIS function was implemented (goat ID as nested factor).

Finally, we assigned all arthropods to functional feeding groups (Supplementary Table S7, see Supplementary Data Sheet 1 for full details): herbivores, predators, parasites, dung-associated arthropods, aquatic arthropods or unknown. We also determined herbivore feeding niche based on the literature: (1) exophages – feed on the plant surface; (2) endophages – feed within plant parts (such as gall formers and arthropods boring or mining leaves, stems, fruit and seeds); (3) unknown. We compared the mean RRA of the functional feeding groups among sampling months using a Friedman’s test (post hoc Wilcoxon signed-rank test, Bonferroni corrected for multiple comparisons). The number of genera per functional feeding group was compared among sampling months using the Pearson chi-square test. Both tests were performed using IBM SPSS software v.25.

While grazing, the goats ingested a wide variety of arthropods, mostly plant-dwelling taxa (88% of both RRA and POO, Figure 1A and Supplementary Data Sheet 3: Supplementary Figure S1A) and a few aquatic taxa ingested while drinking (2% of both RRA and POO, Figure 1A and Supplementary Data Sheet 3: Supplementary Figure S1A). A small proportion of the RRA belonged to dung-associated arthropods (2% RRA, Figure 1A; 3% POO, Supplementary Data Sheet 3: Supplementary Figure S1A), which might have contaminated the feces during collection, and to arthropods whos’ functional feeding group could not be determined (8% RRA, Figure 1A; 7% POO, Supplementary Data Sheet 3: Supplementary Figure S1A).

The plant-dwelling taxa ingested by goats were comprised of a complex food-chain of herbivores, predators and parasites (Figure 1B). In total, the goats ingested species from 48 families of herbivores (93% RRA of PDA, Figure 1B; 90% POO of PDA, Supplementary Data Sheet 3: Supplementary Figure S1B), mainly of the orders Hemiptera, Diptera, Lepidoptera and Coleoptera (Figure 1B). The majority of these herbivores were exophages (77% RRA, Figure 1B; 74% POO, Supplementary Data Sheet 3: Supplementary Figure S1B). The goats also ingested species from nine families of predators (5% RRA of PDA, Figure 1B; 6% POO of PDA, Supplementary Data Sheet 3: Supplementary Figure S1B), mostly spiders (Araneae), and four families of Hymenopteran parasitoids (2% RRA of PDA, Figure 1B; 4% POO of PDA, Supplementary Data Sheet 3: Supplementary Figure S1B). Larvae of chironomids and mosquitoes (Diptera) were ingested while drinking (Figure 1A). All samples that passed quality control (57/57) contained DNA of PDA, suggesting PDA ingestion is frequent.

A seasonal shift in the arthropod community, associated with Mediterranean habitats, was strongly evident in the composition of arthropods detected in the goat feces collected over February–June (ADONIS: RRA, R = 0.256, p = 0.001, Figure 2; POO, R = 0.191, p = 0.001; Supplementary Data Sheet 3: Supplementary Figure S2; Supplementary Table S1). This seasonal shift was also visible in the feces of the six-goat group (i.e., goats with at least three consecutive samples from April–June; ADONIS: RRA, R = 0.251, p = 0.001; Figure 2). Interestingly, Diversity and richness estimates (Shannon H′ and Fisher alpha; calculated based on the six-goat group), were similar among sampling months for all individual goats in the six-goat group (Friedman’s test: Shannon H′, χ²₂ = 1.333, p = 0.513; Fisher’s alpha, χ²₂ = 0.333, p = 0.846; Supplementary Table S8).

The majority of PDA ingested by goats in all sampling months were herbivores (Figure 3A). While the overall relative abundance of herbivores remained steady over the season (Friedman’s test: χ² = 8.898, p = 0.064), endophages were significantly more abundant in May–June (the beginning of summer, Friedman’s test: χ² = 3.234, p = 0.001, Figure 3C). Parasites were more abundant in April (Friedman’s test: χ² = 16.288, p = 0.003), whereas predators, like herbivores, remained relatively steady over time (Friedman’s test: χ² = 8.916, p = 0.063). Although there seemed to be a seasonal trend of ingestion at the genus level for all functional feeding groups (Figures 3B, D), this was non-significant (Chi-square: herbivores, χ²₁ = 3.64, p = 0.457; exophages, χ²₁ = 1.683, p = 0. 794; the remainder had less than 5 observations in some cells thus statistics could not be performed).

Despite feeding together in the same habitat, PDA ingested by individual goats were notably different. The mean number of PDA families ingested by goats varied significantly between individuals in the six-goat group (repeated measures ANOVA: tests of within-subjects’ effects - sampling months, F_2,10 = 0.17, p = 0.846; tests of between-subjects’ effects–individual goats: F_2,10 = 59.827, p = 0.001; Figure 4A). When observing the actual change in the number of families ingested over April–June, some goats (mostly) remained above the median number of families ingested (8) over time, while others remained below this number (Figure 4B). An individual trend was not apparent when comparing the composition of the total arthropod community (ASVs) among all goats or the six-goat group (ADONIS: all goats - R = 0.019, p = 0.279, Supplementary Data Sheet 2: Supplementary Figure S2; six-goat group - R = 0.058, p = 0.445).

Indiviaual (and seasonal) variation especially stood out when examining the ten PDA families with the highest RRA (and POO, Supplementary Data Sheet 3: Supplementary Figure S3; Figure 5; for order abundance for both datasets see Supplementary Data Sheet 3: Supplementary Figure S3). For example, the feces of some goats contained a large proportion (RRA) of Liviidae (plant-parasitic hemipterans), while others contained none at all. Psyllidae (jumping plant lice) and Tephritidae (fruit flies), on the other hand, were present in most of the goats’ feces, while Chrysomelidae (leaf beetles) were present in the feces of only a few goats (Figure 5).

Goats ingested mostly innocuous PDA characterized by reduced mobility. These included three main groups: (1) Endophagous PDA (12% of PDA, Figure 1B) which reside within plant parts from which they are unable to escape, making them vulnerable to ingestion. This group was comprised mainly of leaf miners (Agromyzidae), plant-parasitic mites (Eriophyidae), gall midges (Cecidomyiidae) and boring caterpillars (Lepidoptera). (2) Immature life stages of PDA with mobile adult stages; such as Thrips (Thripidae), hover flies (Syrphidae), herbivorous moths (Pyralidae, Noctuidae) and parasitic wasps (Hymenoptera, the larval stage was likely ingested together with their arthropod hosts). While adult and young PDA stages (i.e., eggs, pupae, larvae) cannot be distinguished from one another using DNA metabarcoding alone, we can assume that the early and immobile life stages were ingested and not highly mobile adults which can escape danger. (3) Small wingless species, such as aphids (Aphididae) and mites (Acari), which are unable to rapidly escape from the plant, like winged or highly mobile species. Nonetheless, certain aphid species are able to escape feeding LMH (Gish et al., 2010; Ben-Ari et al., 2019). The efficiency of the escaping behavior (within and between PDA species) will determine the outcome of the encounter with LMH. Spiders which are common in the study site (Lubin et al., 2011) were also ingested by the goats; namely crab spiders (Thomisidae) and jumping spiders (Salticidae), that are usually found stalking their prey on plant parts such as flowers and fruit. Goats (and LMH in general) might prefer these nutritious plant parts, making PDA present on them more vulnerable to ingestion.

Highly mobile species which are able to rapidly escape the plant when approached by a LMH, such as cicadas, spittlebugs, and grasshoppers, were hardly detected or absent from the fecal samples, despite being common in the study area. Similarly, harmful and noxious PDA [bees and ants, for example, which are also common in the field (Shavit et al., 2009 Finkel et al., 2015)], were also absent from the feces. We recently demonstrated that goats are able to effectively avoid ingesting noxious PDA, such as the spring webworm [Ocnogyna loewii, Lepidoptera (Berman et al., 2017)], while feeding. These caterpillars are covered with long setae which can cause allergic responses in LMH upon contact (Campbell, 2001) and teratogenic diseases if ingested (Webb et al., 2004). While late instar webworms are mobile (moving freely on plants), early instars are not - shortly after hatching the webworms form a common web nest on the vegetation in which they feed and grow. Despite their inability to escape LMH during early instars and their abundance in Mediterranean maquis, especially during the study period (February–June), no webworms (early or late instars) were ingested by the goats. This finding implies that goats are able to sense and avoid caterpillar nests as well as the small webworms within them.

The population dynamics of arthropods in Mediterranean ecosystems follows the seasonality of the Mediterranean climate (Stamou et al., 2004). Furthermore, Mediterranean ecosystems are characterized by strong oscillations in climatic factors (rainfall, temperature and humidity) and resource availability that are impacted by human activity, including grazing. These changes in plant resources greatly influence herbivore assemblages (Li et al., 2022) and the temporal composition of arthropods in the habitat (Stamou et al., 2004; Lazaro et al., 2016; Zhu et al., 2020). We detected a significant seasonal variation in the composition, but not the diversity, of PDA ingested by goats from the end of winter to the beginning of summer (i.e., February to June; Figure 2 Supplementary Data Sheet 3: Supplementary Figures S2A, B; Supplementary Table S8). Meaning, the assemblage structure of PDA ingested by goats changed over the season, yet the diversity and richness indices of these arthropods remained consistent.

Herbivorous and predatory PDA were consistently ingested by the goats through all sampling months (Figure 3A). However, endophagous herbivores were more prevalent in the feces in the end of spring (May–June, Figure 3C), probably since they develop within plant parts providing protection from the hot Mediterranean summer. Parasites were ingested mostly in spring (April, Figure 3A), probably because they depend on the development of their herbivorous hosts which are plentiful during this period due to warm temperatures and increased vegetation growth. Overall, the detection of PDA in the goat feces during February–June was in accordance with the phenology of PDA in Mediterranean habitats during this time period. Future studies may therefore benefit from using seasonal fecal sampling of LMH as a way to infer seasonal changes in arthropod composition in grazed habitats.

By collecting fecal samples directly from marked goats once a month, we were able, for the first time, to observe the extent of individual variation of PDA ingestion in LMH over time. It is important to note that not all goats had a valid sample per month, therefore we focused our analysis on goats with at least three consecutive samples (3 months in a row, i.e., the six-goat group). Despite no significant difference between the six individuals in the overall composition of PDA (ASV’s) detected in the feces (Figure 2), our analysis revealed that the number of PDA families (but not functional feeding groups) ingested by goats (i.e., diversity) varied greatly among individuals (some goats ingested significantly more PDA families than others; Figures 4A, B). In light of our recent studies (Berman et al., 2017, 2019a), it is unlikely that the observed variation is a result of an inability of some goats to detect and avoid PDA ingestion. The fact that the RRA of particular PDA was higher in some goats (such as that of the family Tephritidae which were detected in the feces of all six goats, Figure 5) and that RRA may actually contain some quantitative aspects (Deagle et al., 2018), might suggest, with caution, that different goats ingest different quantities of these PDA. This assumption could be tested when complementing DNA metabarcoding with a quantitative method (Deagle et al., 2005; Guo et al., 2018).

Individual variability in feeding behavior of LMH is influenced by many factors, including availability, distribution and quality of forage or preferred foods, nutritional needs, prior experience with food plants, social interactions (such as competition and learning from conspecifics) and personality (Provenza et al., 2003; Searle et al., 2010; Neave et al., 2018). Since ingestion of PDA by goats is probably incidental, its pattern will likely reflect the dietary choices of each individual, i.e., the plants or plant parts or species consumed by each goat. Individual goats which avoid certain plants based on prior experience, social learning, palatability, or their ability to cope with toxins (Provenza et al., 2003; Neave et al., 2018) will rarely ingest PDA residing on these plants. Personality traits, such as boldness or fearfulness, will influence whether individuals encounter and sample novel food plants (Neave et al., 2018) and the PDA on them. Goats are social animals which feed together under competitive conditions. Thus, intraspecific competition may lead some individuals to exploit resources that are less utilized by conspecifics (Araújo et al., 2011), exposing them to a different variety of PDA. Goat feeding behavior is also influenced by social ranks with dominant individuals having greater access to quality food compared to subordinates (Neave et al., 2018). Such relationships within the herd can affect individual feeding patterns and consequently the variety or amount of PDA ingested.

To date, only one study has examined the extent of PDA ingestion by LMH at the community level using a molecular approach. This study investigated the ingestion of PDA by cattle grazing in Mediterranean grasslands using DNA metabarcoding (Berman and Inbar, 2022). The application of a similar molecular analysis (including the primer set used) to the one used in our current study enabled us to compare between the two.

Goats (i.e., browsers) feed mostly on woody vegetation (fruits, tree and shrub stems and foliage, but also grasses), while cattle (i.e., grazers) feed primarily on grasses and forbs. Beside feed type, these LMH also differ in mouth structure, a fact that may influence bite size and thus feeding rate (Shipley, 1999). Grazers usually have wide muzzles, allowing them to maximize bite size but making it harder for them to select specific plant parts. Browsers, on the other hand, have a narrower muzzle and in some species, such as goats, flexible lips and tongues which allow them to select desired parts but reduce bite size (Lu, 1988; Shipley, 1999). Despite these notable differences, arthropod ingestion was prevalent in both goats and cattle (100 and 76% of feces that passed the quality control contained sequences belonging to PDA, respectively) and their pattern of arthropod ingestion showed high resemblance. Strikingly, goats and cattle ingested similar proportions of herbivores, predators and parasitoids (Figure 6) and a small percentage of aquatic arthropods. The abundance of dung associated arthropods could not be compared between cattle and goats since the time period between defecation and collection was longer in cattle (goat fecal samples were collected at time of defecation), enabling coprophages to arrive and colonize the feces prior to collection. In both studies, PDA characterized by reduced mobility (inability to rapidly escape the plant) were ingested. Harmful or noxious PDA, on the other hand, were largely absent from the fecal samples. This is not surprising as we recently demonstrated that both goats and cattle are able to detect and avoid ingesting noxious PDA (Berman et al., 2017, 2018). Both ingestion patterns also exhibited a similar seasonal shift, with a significant difference in the PDA detected in the feces over time (Berman and Inbar, 2022). Although we did not examine individual variation in cattle, LMH characterized by wide muzzles and large bite size will probably feed in a more homogenous matter than small-muzzled, selective LMH (like goats), reducing the capacity of individual variation in plant consumption (Pansu et al., 2019) and maybe even PDA ingestion. Moreover, goats, which are equipped with flexible lips and tongues, may ingest PDA found in plant parts which are less reachable for other LMH species, including grazers.

The application of DNA metabarcoding on fecal samples collected repeatedly from known individuals enabled us, for the first time, to observe individual variation in the ingestion of PDA by foraging LMH. This finding would be nearly impossible to observe using conventional dietary methods. Despite starting with 18 marked goats, technical issues along the way (notably inability to produce adequate PCR products for part of the samples) prevented us from obtaining a valid sample per month for each individual. Accordingly, we focused our analysis on a sub group of goats with samples from three consecutive months (i.e., the six-goat group). This analysis allowed us to observe individual variation in PDA ingestion at the diversity level (number of different families ingested), but not at the community level (ASV’s). DNA metabarcoding is a strong tool for detecting interactions within food webs (Pringle and Hutchinson, 2020), especially cryptic interactions, such as the presence of endophages and parasitoids which are concealed within their host and would not be detected by trapping or observations. Since the results of DNA metabarcoding studies are affected by recovery biases, contamination of feces by arthropods which were not ingested (i.e., coprophages) may mask actual abundances and interfere with the interpretation of the results. Collecting fecal samples directly from the goats prevented contamination by coprophages, and therefore provides a better picture of the PDA community that was actually ingested.

While DNA metabarcoding enabled us to successfully observe individual variation in the ingestion of PDA by LMH under field conditions, this method has a few disadvantages. For instance, wide range PCR primers may amplify nontarget taxa and may provide less data at the species level (Deagle et al., 2018; Pringle and Hutchinson, 2020; Cuff et al., 2022). We nonetheless chose the Zeale et al. (2011) primers since they enabled us to successfully assign arthropods to functional feeding groups (Berman and Inbar, 2022) and since they are frequently used for arthropod amplification in dietary studies (Alberdi et al., 2018). Future studies should use multiple primer pairs, such as COI and ribosomal markers, to avoid amplification bias (Cuff et al., 2022). Another disadvantage of DNA metabarcoding is that it only provides semi quantitative data (Deagle et al., 2018; Pringle and Hutchinson, 2020), whether using occurrence or RRA data. Both types of analyses used in this study led to a similar outcome in the results and their interpretation (Supplementary Data Sheet 3).