Located in the North of Lebanon, Horsh Ehden Nature Reserve (HENR) (34°18'17.79"N, 35°58'35.33"E) is a highly biodiverse protected area in the Eastern Mediterranean Region covering ca. 1,775 hectares encompassing an elevational gradient from 1100 to 2300 m. To date, 1,058 plants species including 39 tree species, 156 bird species, 4 amphibian species, 19 reptile species and 300 fungi species have been recorded in this reserve (MoE/UNDP/UL, 2004). Additionally, 27 mammal species have been reported, of which 12 were recently identified in the reserve by using an environmental DNA-based approach (Boukhdoud et al. 2021).

One hundred and ninety-five scat samples belonging to 15 mammal species were collected between 2018 and 2020 within HENR (Table 1). All scats found (fresh and dry) were collected randomly from different zones of sparse vegetation, dense forest, and trails at different altitudes ranging from around 1,200 to 2,220 m. Only cliffs were not visited. Samples were collected, geo-located (Suppl. material 2: Appendix 2), and preserved in a sterilized container with silica gel to avoid DNA degradation. Scat samples from five mammal species (four carnivore species and the hare) were collected in all four seasons; however, for the remaining 10 species, we only obtained samples in one, two, or three seasons (Table 1). The internal part of the scat samples was homogenised individually, and approximately 100 mg of the homogenised scats was used in DNA extraction.

DNA extractions were performed after one to 3 days maximum of the sample collection using a QIAamp Fast DNA Stool kit (Qiagen, Valencia, CA, USA) following the manufacturer’s protocol. A negative control was included in each extraction batch. To identify the mammal species that produced each collected scat, we performed DNA barcoding using the 12S rRNA marker as described in Boukhdoud et al. (2021).

To identify the plant species consumed by each mammal species, we amplified the P6 loop region of the chloroplast trnL (UAA) intron (10–143 bp) from the scat-derived DNA extracts. Our first PCR amplified the trnL region of interest using the g and h primers described in Taberlet et al. (2007) combined with overhang adapter sequences (A. Welch, personal communication). The following primers were used:

trnL-g, 5’-ACACTCTTTCCCTACACGACGCTCTTCCGATCT-GGGCAATCCTGAGCCAA-3’ and trnL-h, 5’-GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT-CCATTGAGTCTCTGCACCTATC-3’ (Taberlet et al. 2007). Each DNA sample was amplified in three technical PCR replicates. We included a negative water control for each PCR master mix. The PCR reaction was 25 μl in total volume, including 5 μl Phusion HF Buffer (Thermofisher), 0.2 μl of Phusion enzyme, 0.5 μl of dNTPs (10 mM), 1.25 μl of forward and reverse primers (10 mM), 2.5 μl of DNA template, 1 μl of bovine serum albumin (BSA), and 13.3 μl of sterile water.

The PCR cycling profile was as follows: first denaturation at 98 °C for 30 s, followed by 35 cycles of 98 °C for 10 s, 57 °C for 25 s, and 72 °C for 30 s, and after the last cycle hold at 72 °C for 7 mins.

The PCR products were cleaned using MagnaBind Carboxyl Derivatized Beads (Thermo Scientific) at a 1.7× ratio. After cleaning, products were visualized on a 1% agarose gel.

These triplicate PCR reactions of the same sample were then pooled together. We then performed a second PCR to add iTruSeq dual indices and sequencing adapters (Glenn et al. 2019). This index PCR reaction was performed in 25 μl total volume including 5 μl Phusion HF Buffer (Thermofisher), 0.25 μl of Phusion enzyme, 0.5 μl of dNTPs (10 mM), 2.5 μl of forward and reverse primers (5 mM), 4 μl of DNA template and 10.25 μl of sterile water. Water was used as the negative PCR control.

The cycling profile was as follows: denaturation at 98 °C for 45 s, followed by 13 cycles of 98 °C for 20 s, 61 °C for 30 s, and 72 °C for 30 s, and after the last cycle hold at 72 °C for 7 mins.

Indexed PCR products were visualized on a 0.5% agarose gel and purified with QIAquick Gel Extraction Kit (Qiagen, Hilden, Germany) following manufacturer’s protocol. DNA was quantified using a Qubit dsDNA HS Assay Kit following the manufacturer’s protocol (Thermo Fisher). Samples were pooled together at equimolar concentrations and sequenced on an Illumina MiSeq using a v2, 300 cycle kit, with 150 bp paired-end reads at the Smithsonian Center for Conservation Genomics (CCG) in Washington, DC, USA.

We downloaded demultiplexed sequence data from the BaseSpace Server (Illumina). Sequence reads were then uploaded to mBRAVE, (Multiplex Barcode Research and Visualization Environment) (Ratnasigham 2019), an online platform for analyzing and visualizing DNA metabarcoding data (http://mbrave.net/). This platform is directly linked to the BOLD system database to facilitate species assignment. Parameter values selected on mBRAVE are automatically applied to all runs uploaded in the same project dataset. To merge paired-end reads, a minimum 30 bp overlap between the forward and reverse reads was required, allowing up to 5 nucleotide substitutions. Low quality sequences were removed if the average quality value (QV) was less than 25, allowing maximum 1% of nucleotides with QV value <20 or <10. Sequences shorter than 70 bp were also removed. After filtering, sequences were dereplicated and clustered as Operational Taxonomic Units (OTUs) using a 2% similarity threshold. OTUs were assigned to BINs (Barcode Index Number System) with a 1% ID distance threshold to the reference sequences. Reference sequences were imported from the BOLD database. Reference sequences for the trnL marker were previously generated for the most common tree and shrub species in Lebanon and added to both the BOLD and GenBank database systems (Boukhdoud et al. 2020). The sequence of each OTU not assigned to a species was blasted on the GenBank nonredundant nucleotide (nr/nt) database (accessed January 2021). The taxonomic identification was left to genus level if the similarity was less than 98% and to family level if the similarity was less than 97%. For each run, a summary was generated by the software to record filtering parameters and other criteria including BINs and OTUs counts, QV value distribution graphs, GC composition and sequence length distribution. For each sample, to remove possible sequencing errors or DNA contamination, we inspected the profile of negative controls and compared their read counts with that of dietary samples to exclude samples of possible low-quantity DNA and we excluded the OTUs clustering less than 1% of the sample’s total reads.

We calculated Alpha diversity (Shannon index) by seasons and by animals using PAST software (Hammer et al. 2001). Kruskal-Wallis test (Kruskal and Wallis 1952) was performed using R software to evaluate the statistical difference of plant species richness and diversity between animals and seasons. Plots were created using the package “ggplot2” (Wickham 2009). p values less than 0.05 were considered to be significant.

Prior to DNA extraction, plant seeds, when present, were recovered from the scats and cleaned using diluted ethanol to remove scat residues in order to identify them. The reference seed collection of the Jouzour Loubnan seed bank was used to identify seeds. This seed bank conserves more than 26 million seeds for 100 different taxa (www.lebanon-flora.org; www.jouzourloubnan.org).

DNA extraction and amplification yielded sufficient DNA for further analysis in 177 samples out of 195 (90.76%). The Illumina MiSeq sequencing run generated a total of ~12.1 million paired-end sequence reads; the average number of reads per sample was 15,000 (1,459–86,548). The run’s QV score distribution was 40 and the average sequence length was 80 bp. After filtering, ~3.6 million reads remained. In total, after excluding OTUs representing less than 1% of the sample’s total reads, the analyzed samples produced 133 OTUs of which 44.36% were identifiable to species level, 45.11% to genus level, and 6.77% to the family level.

Five OTU sequences were unidentifiable and did not match any of the species available in public databases.

The obtained 133 OTUs represent at least 133 plant species from 54 different families (Suppl. material 1: Appendix 1). Among the identified species there were angiosperms including Rosaceae species and gymnosperms such as Cedrus libani and Pinus brutia. The OTUs also included different functional types including trees (e.g., Acer tauricolum and Quercus spp.), shrubs (e.g., Berberis libanotica and Sambucus ebulus), forbs and grasses (e.g., Hordeum leporinum and Agropyron panormitanum).

Many identified species were not previously reported in HENR including Sesamum indicum, Ficus carica, Morus alba and Hibiscus trionum. In addition, some identified species are not part of the native flora of Lebanon. We found some ornamental plants such as Lobelia and several other species (e.g., Camellia, Actinidia, Musa, Diospyros, and Arrachis).

The results showed that the red fox (Vulpes vulpes) has the most diverse diet. It consumed at least 84 plant species from 40 different families across the year, unlike the wild cat (Felis silvestris), which consumed only a few plant species. The wild cat consumed Medicago and Prunus species in summer, and Cerastium, Hordeum leporinum and Prunus species in spring. The golden jackal (Canis aureus), grey wolf (Canis lupus), Cape hare (Lepus capensis), beech marten (Martes foina) and the wild boar (Sus scrofa) also have diverse diets; they consumed at least 32, 47, 39, 43, and 45 plant taxa, respectively, from more than 20 different families (Figs 1, 2, Suppl. material 1: Appendix 1).

Figure 1.

Food webs of 15 mammals and the consumed plant families.

Figure 2.

Proportion of sequence reads amplified for each plant family across the seasons. Taxa for which reads number is > 0.5% of the total sequence reads are included.

Six plant family representatives were consumed and shared by the majority of the species: Rosaceae, Poaceae, Apiaceae, Fabaceae, Fagaceae and Berberidaceae. Only members of the Rosaceae family were consumed by all mammals. Species from the Poaceae family and Apiaceae family were consumed by all mammals except the Caucasian squirrel (Sciurus anomalus) and the wild cat, respectively. Some plant families and species were consumed by only one mammal species, including common ivy (Hedera helix, Araliaceae) consumed by wild boar, Actinidia sp. (Actinidiaceae) consumed by the golden jackal, and the commonly called flower-of-an-hour (Hibiscus trionum, Malvaceae) consumed only by the eastern broad-toothed field mouse (Apodemus mystacinus). Many species including Chenopodium album (Amaranthaceae), Corydalis solida (Papaveraceae), Daphne oleoides (Thymelaeaceae), Gagea sp. (Liliaceae), Ixiolinion tataricum (Ixioliriaceae), Onosma sp. (Boraginaceae), and Thesium sp. (Santalaceae) were consumed by red foxes only. The grey wolf was the only mammal to consume the following species: Camellia sp. (Theaceae), Solanum luteum (Solanaceae) and Musa sp. (Musaceae). Cupressaceae species including Juniperus foetidissima were a plant component of the beech marten’s diet only. Pinaceae species such as Cedrus libani and Pinus brutia were, respectively, a component of the red fox and grey wolf’s autumn diets, and they also represent a food resource for the Caucasian squirrel. Some plants were consumed by only two of the 15 mammals including Acer tauricolum, eaten by the red fox and beech marten, Vitis vinifera from the Vitaceae family, consumed by the golden jackal and red fox, and Salix libani, consumed only by the grey wolf and wild boar. In addition, the species Eremurus spectabilis (Xanthorrhoeaceae) was consumed by both grey wolves and Indian porcupines (Hystrix indica). The Lebanese cedar was detected only in Caucasian squirrel and red fox samples (Fig. 1, Suppl. material 1: Appendix 1).

The mean number of plant OTUs obtained across mammal species within seasons was 25 in autumn, 12 in winter, 17 in spring and 17 in summer. For the majority of species, the highest dietary diversity was detected in autumn; low plant dietary diversity was observed in winter compared to the other seasons (Fig. 2).

Rosaceae species sequences were the most numerous in the majority of the species’ scats, and this pattern was especially prominent in autumn and winter. Rosaceae species are the most numerous in the grey wolf’s and beech marten’s diet across all four seasons. In winter, members of the Rosacea family represent 89.1% and 87.2% of the sequences derived from grey wolf (number of species n ≤ 4) and beech marten (n ≤ 6) samples, respectively. They also represent 33% of sequences from the Indian porcupine samples in autumn and 53.2% in winter (n ≤ 4). Rosaceae is also the most speciose plant group in the diet of the field mouse and Caucasian squirrel in the summer and winter seasons, respectively. Rosacea species sequences were most common in the scats of goat (Capra hircus) and cattle (Bos taurus) in spring, and of Eurasian badger (Meles meles) in summer. In spring, Hypericum spp. (Hypericaceae) is consumed by the southern white-breasted hedgehog (Erinaceus concolor) (78.6% of total reads). The least weasel (Mustela nivalis) consumed mainly Cerastium spp. (Caryophyllaceae) in this season (Fig. 2).

For several species, our results showed seasonal shifts in dietary plant components. The diet of the red fox was dominated by Rosaceae species in autumn and winter, but shifted to be dominated by Poaceae species (mainly Hordeum leporinum) in spring and Fabaceae species in summer, especially Onobrychis sp. The diet of the wild boar included mainly Quercus species (Fagaceae) in winter and summer, but shifted to include more Rosaceae species in autumn. Species from the Fabaceae family had higher representation in Cape hare samples during winter and summer compared to spring, when its diet included mostly members of the Poaceae family, and autumn, when Sedum sp. (Crassulaceae) was mainly detected. Concerning the golden jackal, its diverse diet included primarily Ficus carica (Moraceae) in autumn, Rosaceae species in winter, Quercus species in spring and Medicago sp. (Fabaceae) in summer (Fig. 2).

Some plant species were identified in scats in only one season. Anthriscus lamprocarpa, Ribes orientale, Salix libani and Eremurus spectabilis, were only found in autumn. Species including Aethionema coridifolium, Ononis natrix, Corydalis solida, Androsace villosa and Cupressaceae species were identified only in spring. Genista libanotica, Melica angustifolia and Anemone blanda were found only in winter. Satureja cuneifolia, Hibiscus trionum, Morus alba, and Myrtus communis were identified only in samples collected in the summer.

In the Figure 3, we compared the plant species consumption of four carnivora species (golden jackal, grey wolf, beech marten and red fox).

Figure 3.

Plant families consumed across seasons by 4 mammals of the order Carnivora.

Two plant families were consumed by all Carnivora species in all seasons: Poaceae and Rosaceae. Species from the Fabaceae family (e.g., Medicago sp.) were eaten by all four mammals in all seasons, except the golden jackal in the winter. Figs (Ficus carica) from the Moraceae family constituted a component of the diets of golden jackal and grey wolf year-round. In autumn, species from six families were consumed by the four Carnivora species, including Apiaceae, Fabaceae, Fagaceae, Moraceae, Poaceae, and Rosaceae. Lactuca sp. (Asteraceae) and Berberis libanotica (Berberidaceae) species were consumed by grey wolves, beech martens, and red foxes; Hypericum sp. (Hypericaceae) was identified in the scats of all species except the beech marten. In spring, species from eight plant families were consumed by all four carnivore species, including Cerastium sp., Ficus carica, Quercus infectoria, Quercus coccifera, and Agropyron panormitanum. Compared to the three other Carnivora species, the red fox scats contained the most diverse plant taxa (i.e. red fox scats contained the highest number of plant families n = 40) (Fig. 3).

The alpha diversity metrics (Shannon index and observed taxa) revealed the highest diversity and richness of consumed plant species in autumn and the lowest in winter (Fig. 4a and b). Alpha diversity values also showed that plant species richness and diversity were the highest in the red fox’s diet and the lowest in the golden jackal’s diet (Fig. 4c and d). Plant species richness differed significantly between animal species (Kruskal-Wallis test, p<0.05). However, in contrast, no significant difference was observed among seasons (Kruskal-Wallis test, p>0.05) (Fig. 4).

Figure 4.

Comparison of alpha diversity of the consumed plant species in different seasons and between animal species. Shannon index (a) and (c), observed species (b) and (d).

Plant seeds were found in 33% of the collected scat samples. These samples belong to the following six mammal species from three different orders: the golden jackal, grey wolf, Cape hare, beech marten, wild boar, and red fox. Seeds were most present in the scat samples collected during summer-autumn season.

The identified seeds belong to species from several families including Rhamnus cathartica (Rhamnaceae), Ficus carica (Moraceae), Vitis vinifera (Vitaceae), and Rosaceae species (e.g Crataegus spp., Malus trilobata, Rosa canina., Prunus spp., Sorbus spp., and Pyrus syriaca). Many samples contained seeds from more than one plant species. Undigested fruits were also identified in several samples (Fig. 5).

Figure 5.

Plant seeds detected in scat samples.

Syrian pear seeds (Pyrus syriaca) were only identified in the grey wolf samples. Grape seeds (Vitis vinifera) were only found in red fox and the golden jackal samples. Only red fox scat samples contained buckthorn seeds (Rhamnus cathartica). All plant species for which we identified seeds in scat samples were also detected by DNA-barcoding.