Exploring winter diet, gut microbiota and parasitism in caribou using multi-marker metabarcoding of fecal DNA

doi:10.21203/rs.3.rs-4455656/v1

Download PDF

Article

Exploring winter diet, gut microbiota and parasitism in caribou using multi-marker metabarcoding of fecal DNA

https://doi.org/10.21203/rs.3.rs-4455656/v1

This work is licensed under a CC BY 4.0 License

You are reading this older preprint version

Read the latest preprint version →

In conservation strategies, getting precise and repeatable information on the species’ diet and health without relying on invasive or laborious methods is challenging. Here, we developed an efficient and non-invasive workflow for the sequencing and analysis of four taxonomic markers from fecal DNA to characterize the gut microbiota, parasites, and plants and lichens composing the winter diet of caribou (Rangifer tarandus), Canada's most iconic endangered species. Sequencing of the 18S rRNA gene of eukaryotes from seven locations in Manitoba and Saskatchewan, Canada, allowed for the detection of several parasites in caribou feces but provided limited information about plants and lichens - only algal symbionts were detected. An improved sensitivity and higher taxonomic resolution for plants and lichens was achieved through the sequencing of the ITS2 region, which revealed a rich winter diet in caribou as well as differences among sampling locations. Sequencing of the 16S rRNA gene of prokaryotes highlighted the relationships between the gut microbiota and both the diet and parasites. Overall, our results showed the potential of this multi-marker DNA metabarcoding workflow as an efficient tool to generate relevant information on the diet and health of caribou populations and provide insights into the species biology and ecology.

Biological sciences/Biological techniques/Bioinformatics

Biological sciences/Biological techniques/Sequencing

Biological sciences/Ecology/Boreal ecology

Biological sciences/Ecology/Molecular ecology

Biological sciences/Microbiology/Communities/Metagenomics

Earth and environmental sciences/Ecology

Earth and environmental sciences/Ecology/Boreal ecology

Earth and environmental sciences/Ecology/Conservation

fecal DNA

metabarcoding

gut microbiota

lichens

caribou parasites

Rangifer tarandus

With the ongoing global climate change, many wildlife species are vulnerable or at risk of extinction^1,2. Among them, caribou (Rangifer tarandus), an iconic wildlife species, is declining across its circumpolar distribution³. Boreal caribou, in particular, ranges across the entire boreal ecozone and is currently listed as threatened under Canada’s Species At Risk Act, with only fifteen out of fifty-one local populations being self-sustaining⁴.

The main reasons identified for the decline of boreal caribou are disturbances and modifications to its habitat due to human activities^5,6, in addition to direct and indirect effects of climate change⁷. Modifications to caribou habitat can favor the presence of other ungulates species and predators^8,9 and affect the composition of the ground vegetation. Boreal caribou diet is composed of plants and lichens, the latter being the main source of energy during winter¹⁰. As a result, caribou forage in peatland, open lichen woodland and old growth coniferous forests^11,12, three ecosystems that are experiencing modifications in response to climate change^13,14. Suitable nutritional gain is crucial for caribou survival, particularly during the winter season, when foraging might be difficult, and for pregnant females that rely on lichens, a highly digestible source of energy¹⁵. Deepening our knowledge of the caribou diet could help us better protect its habitat, but so far this has proved difficult: collar camera studies are complex and involve low number of individuals¹⁶ and previous attempts to characterize the diet based on environmental DNA from fecal samples failed to recover lichens^17,18.

The ability of caribou to maximize nutrient absorption is firmly linked to microbial activity, and potential changes in diet could impact its gut microbial community structure and functions¹⁹. The composition of the gut microbiota is also linked to health and metabolism in numerous mammal species^20–22 as well as range size in ungulates²³. Moreover, climate change and modification of habitat could impact caribou-parasite interactions^24–26. A higher presence and abundance of parasites could be linked to a decrease in fitness^27–29 and be particularly deleterious for declining populations. Finally, low quality nutrition can lead to greater susceptibility to parasite infection and create a vicious circle where animals with poor body condition are more likely to be infected that in turn affects their fitness³⁰. In this context, simultaneous knowledge about the animal’s diet, gut microbial community, and parasites would be beneficial.

Metabarcoding of DNA isolated from feces is an effective and non-invasive way to describe the gut microbiota, as exemplified by previous studies on caribou^19,31,32, and can also be applied to the detection of parasites that inhabit mammals’ digestive tract^33–35. It is also commonly being used to describe the diet of various organisms^36–39, but so far with limited success for caribou. Indeed, while the plant diet was successfully characterized by metabarcoding of caribou fecal DNA^17,18,40, previous studies either did not address the lichen diet⁴⁰ or failed to do it because of methodological issues¹⁸ or lack of sequences in databases¹⁷, leading to an important gap that this study is proposing to fill.

While metabarcoding is increasingly being used to describe diet, gut microbiota and parasites, the focus of each study is generally limited to one or two of these aspects and the method is applied to small sample sizes. Moreover, variations in the methodology makes merging of datasets from different studies challenging. For an endangered species with a very large geographic range like caribou, a standardized and high-throughput workflow should be developed to limit over-interpretation of differences among studies that could arise from methodological biases⁴¹. To be comprehensive, the workflow should cover various aspects of caribou ecology and allow for the integration of multiple datasets. Metabarcoding of multiple taxonomic markers, including some allowing to detect lichens, is a promising approach to achieve this goal and could be a great complement to ongoing genetic surveys, for which non-invasive fecal collection campaigns are already ongoing through the EcoGenomics research program.

In this study, we propose a standardized workflow, from sampling to data analysis, to describe caribou diet, gut microbiota and parasites from caribou feces using the ITS2 regions of fungi and plants, the 16S rRNA gene of prokaryotes and the 18S rRNA gene of eukaryotes, as taxonomic markers. The workflow was developed using fecal pellets from 96 individuals among seven locations across central Canada. Fecal pellets were collected during the winter season to ensure optimal conditions for the preservation of DNA and advance knowledge related to the winter diet, which is critical for caribou but remains poorly documented¹⁰.

1. Sample selection

Fecal pellets from 96 individuals sampled in seven locations in Manitoba and Saskatchewan, Canada, were selected from an existing extensive collection held by the EcoGenomics Team and used to generate population demographic parameters⁴². Samples from four to eight males and females from each sampling location were included in the analysis (Fig. 1, Supplementary Table S1). Individuals were identified by genotyping based on variable microsatellite loci amplified from fecal DNA, as previously described^43–46, and caribou-specific primers Zfx/Zfy were used for sex identification⁴⁴. Samples were shipped on dry ice and kept frozen at -80°C.

2. Selection of taxonomic markers and primers

Appropriate taxonomic markers and corresponding primers were selected to ensure the detection of group of organisms of interest from caribou fecal DNA. Broad coverage 16S rRNA gene primers targeting the V4-V5 regions⁴⁷ were selected to characterize the gut microbiota, which is mainly composed of bacteria in ruminants^31,48. Parasites colonizing the caribou digestive tract are members of various taxonomic groups, including the Apicomplexa, Archamoebae, Platyhelminthes, Nematoda, and Arthropoda⁴⁹. To simultaneously detect all these organisms, broad coverage primers targeting the 18S rRNA gene of eukaryotes⁵⁰ were selected. These primers have the potential to also inform on the caribou’s winter diet, which is composed of eukaryotic organisms: lichens ̶ a symbiotic relationship between a fungi and at least one photosynthetic partner, cyanobacteria and/or green algae, and plants¹⁰. However, more specific markers are commonly used to identify lichens and plants, so complementary sets of primers were selected for these groups. The identification of lichens is based on the taxonomy of the fungal symbiont⁵¹ and can be achieved through the sequencing of the nuclear ribosomal internal transcribed spacer (ITS)⁵², and primers targeting the ITS2 region were selected^53,54. For plants, no single widely used marker gene is currently accepted, but sequencing of the plant ITS2 region has been shown suitable for a large variety of plants, allow for an appropriate taxonomic identification at the genus level⁵⁵, and is supported by well-developed gene databases⁵⁶. Previously published primers targeting the plant ITS2 were therefore included in the selection⁵⁷. A list of al primers used in this study is provided in Supplementary Table S2.

3. DNA isolation, library preparation, and sequencing

For each individual sample, three frozen fecal pellets were ground in liquid nitrogen using a mortar and pestle. DNA isolation was performed on 100 mg of ground homogenized samples, using the DNeasy PowerSoil kit (QIAGEN, Netherlands) according to the manufacturer’s protocol (May 2017) for QIACubes instruments (QIAGEN). The Qubit 3.0 fluorometer (Thermo Fisher Scientific, MA, USA) was used for the quantitation of DNA using the dsDNA BR (broad range) assay kit (Life Technologies, CA, USA). DNA extracts were diluted to 5 ng/µL with 10 mM Tris pH 8.0 solution using the QIAgility automated system (QIAGEN).

Library preparation for Illumina sequencing was conducted following the manufacturer’s protocol (16S Metagenomic Sequencing Library Preparation, Part 15044223 Rev. B) with some modifications. The 5’ end of the primers (Supplementary Table S2) was flanked with one or the other of the following Illumina overhang adapter sequences: TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG (forward overhang) and GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG (reverse overhang). PCR reactions were set up by mixing 25.0 µL of 2X HotStarTaq Plus Master Mix (QIAGEN, Germany), 0.5 µL of each 10 µM HPLC-purified (or desalted) primer (Invitrogen, USA), 19 µL of RNase-free water (QIAGEN) and 5.0 µL of DNA diluted at 5 ng/µL, for a total volume of 50 µL. Thermal cycling conditions for the 16S rRNA gene were as follows: initial denaturation at 95°C for 5 min; 35 cycles (40 cycles for the 18S, fungal ITS2 and plant ITS2) at 94°C for 45 s, 50°C (47°C for 18S, 50°C for fungal ITS2 and 49°C for plant ITS2) for 45 s, 72°C for 1 min; and a final elongation step at 72°C for 10 min. PCR products were visualized on GelRed-stained 1% agarose gel, using the ChemiGenius Bioimaging System (Syngene, Cambridge, UK). Target fragment sizes were ~ 479 bp for the 16S rRNA gene, ~ 484 bp for the 18S rRNA gene, ~ 347–487 bp for the fungal ITS2 region and ~ 230–378 bp for the plant ITS2 region.

PCR products were purified using AMPure XP reagent (Agencourt, Beckman Coulter Life Science, CA, USA) according to Illumina’s protocol, with gentle pipetting as the preferred option for mixing. Combinatorial dual-index barcodes were added by amplifying 5 µL of the purified PCR product with 25 µL of 2X KAPA HIFI HotStart Ready mix (Roche, South Africa), 5 µL of each Nextera XT index primer (Nextera XT Index kits v2, Illumina Inc., CA, USA) and 10 µL of RNase-free water, for a total volume of 50 µL. Thermal cycling conditions were as follows: 3 min at 95°C; 10 cycles of 30 sec at 95°C, 30 sec at 55°C, 30 sec at 72°C; and a final elongation step of 5 min at 72°C. The resulting dual-indexed libraries were purified with the AMPure XP reagent following the procedure described above, quantified using the SYNERGY Mx system (BioTek Instruments, VT, USA) and pooled at equimolar ratio. To compensate for low base diversity in libraries, 15% of PhiX Control V3 (Illumina) was spiked in the pool. Sequencing was performed using a Miseq v3 600 cycle kit on an Illumina MiSeq system at the Next-Generation Sequencing Platform, Genomics Centre, CHU de Québec-Université Laval Research Centre, Québec, Canada.

4. Pipeline development and database optimization

The pipeline Q2Pipe was developed to ensure a user-friendly and standardized treatment of amplicon sequencing data with QIIME2⁵⁸. Briefly, QIIME2 and its dependencies were loaded to a Singularity image (https://zenodo.org/records/4667718), which allows users to work in a batch environment by installing Singularity, downloading the Singularity image and cloning the Q2Pipe Github repository (https://github.com/NRCan/Q2Pipe). Users can also build the Singularity image by using the image recipe provided in the Q2Pipe repository. The newly built pipeline allows for the selection and modification of all options available in QIIME2 in a single text file which can easily be shared among users and published, while also offering the option to run secondary external programs like FUNGuild⁵⁹, directly from the pipeline. For compatibility purposes, Q2Pipe can also be used with the QIIME2 Anaconda installation method if the user installs the external dependencies in the Anaconda virtual environment.

As several errors have been detected in the genus name of Eukaryota within the SILVA SSU database 138–99⁶⁰, a custom python script was used to screen and correct taxonomical disparities between the name at level 6 (i.e. Genus) and binomial name at level 7 (i.e. species). This correction affected 56% of the Eukaryota genera (not applied to Bacteria and Archaea). Furthermore, a low number of entries for lichens were observed in the UNITE database⁶¹, which is intended for the molecular identification of fungi based on sequencing of the ITS region. To improve the taxonomic assignment for this group, UNITE was enriched with sequences from the OLICH database⁶² as well as from twenty lichen specimens collected in Newfoundland, Canada, and formally identified based on anatomical/morphological features⁶³. Protocols for DNA isolation, amplification, sequencing and bioinformatic treatment of sequences from lichen samples was as described above for fecal samples. Only ITS2 region sequences corresponding to the genus identified by anatomical/morphological features and representing more than 10% of sequences from a sample were conserved. This resulted in twenty-two lichen sequences from nineteen species that have been added to the customized database: Alectoria sarmentosa, Cetraria aculeata, Cladonia amaurocraea, Cladonia arbuscula, Cladonia gracilis var gracilis, Cladonia maxima, Cladonia rangiferina (three sequences from two specimens), Cladonia squamosa, Cladonia stellaris, Cladonia uncialis, Lasalia papulosa, Lobaria pulmonaria, Lobaria quercizans, Platismatia norvegica, Sphaerophorus globosus, Stereocaulon tomentosum, Tuckermannopsis americana and Usnea filipendula. Lichen sequences are available on NCBI under OQ792092 to OQ792112 accession numbers. Lichen sample metadata are available on NCBI under BIOProject PRJNA904780, BIOSamples numbers SAMN34138716 to SAMN34138735.

5. Bioinformatics analyses on sequences from caribou fecal samples

Bioinformatics analyses for each sequence dataset obtained from caribou feces were performed using QIIME2, version 2020.8 within Q2Pipe. Option files for each gene analyzed can be found as Supplementary Material 1. Although parameters were slightly different for each taxonomic marker, the overall procedure was similar. Briefly, demultiplexed, per-sample, gzipped fastq files were imported and sequencing reads were truncated at their 5’ and 3’ ends based on 1) the length of the primer (16S and 18S) or specific primer sequence (ITS2) and 2) per base sequence quality score. The QIIME2 plugin DADA2⁶⁴ was selected as the denoising method, which allowed for filtering, dereplication, merging paired-end reads and chimera identification. This resulted in inferences of amplicon sequence variants (ASVs). The least abundant ASVs were filtered out if their frequency was less than 0.05% of the mean ASV frequency. Naive Bayes trained classifier were generated from the corrected and amended SILVA SSU database 138–99⁶⁰, the amended UNITE v8.3 dynamic fungal database⁶¹, and the PLANiTS v2020.3 database⁵⁶. Taxonomic assignment of 16S and 18S rRNA gene ASVs (SILVA database), fungal ITS2 ASVs (UNITE database) and plant ITS2 ASVs (PLANiTS) was then performed using the QIIME2 plugin feature-classifier classify-sklearn, which is a machine-learning-based classification method that requires trained classifiers. For the 16S rRNA gene, ASVs assigned to Eukaryota, mitochondria, chloroplast and unassigned ASVs at the kingdom level were filtered out, while ASVs from Bacteria, Mammalia and unassigned ASVs at the kingdom level were filtered out for the 18S rRNA gene. For plant and fungal ITS2 regions, taxonomic assignment followed the same process as for the 16S and 18S rRNA genes but using respectively the and the customized, respectively, whereas taxonomic filtering was performed to keep only either Streptophyta or Fungi. ASV tables of the 18S rRNA gene and fungal ITS2 region were further analyzed using FUNGuild⁵⁹ to infer the ecological guilds of fungi based on their taxonomy. Raw sequence data are available in the Sequence Read Archive (SRA) of the NCBI under accession numbers SRR23356982 to SRR23357367. Sample metadata are also available on NCBI under BIOSamples numbers SAMN32246636 to SAMN32246732, BIOProject PRJNA904780.

6. Characterization of the diet and detection of potential parasites

Potential components of the diet were filtered from non-rarefied ASV tables among plants (plant ITS2 and 18S datasets), lichen-forming fungi (fungal ITS2 and 18S datasets) and lichen-forming algae (18S dataset). Lichen-forming fungal ASVs were identified based on FUNGuild trophic mode assignment, while ASVs from lichen-forming algae were identified based on the list published by Sanders & Masumoto⁶⁵. Following steps were carried out on filtered datasets containing only ASVs belonging to these potential components of the diet. To increase the reliability of our data and discard potential contaminants, ASVs were considered as part of the diet when accounting for more than 1% of the sequences in a sample⁶⁶. Then, relative abundance data was converted to presence/absence, as numerous factors may lead to bias in sequence representation of food items in DNA isolated from fecal samples⁶⁷. The frequency of occurrence (i.e. the number of samples that contain a given food item, referred to here as FOO), as proposed by Ando et al.⁶⁸, was used as a semi-quantitative approach to underline diet differences between locations. Plant and lichen ASVs were classified according to their morphology (hereby called classes: herbaceous, non-vascular or lignified for plants; crustose, epiphyte foliose, epiphyte fruticose, foliose or fruticose for lichens). The number of ASVs for the plant and lichen components of the diet as well as for each class were calculated and used as a proxy for the richness of the diet.

To identify parasites, ASVs belonging to the Apicomplexa, Archamoebae, Platyhelminthe, Nematoda and Arthropoda were extracted from the non-rarefied 18S rRNA ASV table; ASVs from these groups were then listed and parasites were identified based on the list of caribou parasites from Tryland and Kutz⁴⁹. ASVs assigned to the Entamoeba genus, a common parasite in mammals, were also included in the parasite list even if, to the best of our knowledge, it has never been reported as a caribou parasite. Moreover, ASVs assigned to the Blastocystis genus, a ubiquitous anaerobic protist which develops in the gastrointestinal tract of a wide variety of hosts such as mammals (including some ungulates), birds and reptiles^69,70 were also included, although it belongs to the Stramenopila lineage⁷¹ and was not identified through our initial taxa selection. As errors in the taxonomic assignment of Eukaryotes was previously observed with the SILVA database, each taxonomic assignment was verified by conducting a BLAST analysis of ASV sequences against the NCBI Genbank database. When the BLAST results were not concordant with the taxonomic assignment obtained at the genus level by our bioinformatic pipeline, the sequences were discarded from the list. As a result, all sequences assigned to the Cryptosporidium genus were discarded as they matched to uncultured eukaryote clones, the first Cryptosporidium hit in NCBI sharing only around 91% similarity with sequences from our dataset. Overall, five genera of potential parasites were detected and included in further analyses. As for the diet, parasite detection was described by the frequency of occurrence and the number of ASVs assigned to each genus.

7. Statistical analyses

All statistical analyses were performed in R (4.2.2 Patched) using 0.05 as a significance threshold value.

For alpha and beta diversity analyses, datasets were normalized by rarefaction in QIIME2 to control for uneven sequencing across samples⁷² (see Q2Pipe Option files in Supplementary Material 1 for details). Optimal library size for each taxonomic marker was selected considering both the rarefaction curves and the number of reads of the smallest libraries. Sampling depth was set at 18,430 sequences for the 16S rRNA gene, resulting in 1,769,280 sequences clustered into 7,192 ASVs observed among 96 samples (100%); 2,750 sequences for the 18S rRNA gene dataset, resulting in 264,000 sequences clustered into 5,174 ASVs observed among 96 samples (100%); 3,700 sequences for the plant ITS2 region, resulting in 351,500 sequences clustered into 1,585 ASVs observed among 95 samples (99%); and 7,100 sequences for the fungal ITS2 region, resulting in 674,500 sequences clustered into 7,326 ASVs observed among 95 samples (99%). The ASV richness (Chao1 index) and diversity (Shannon index) for all markers were estimated with the “vegan” package. To evaluate differences in these indices as well as in richness (i.e., number of ASVs) of different functional groups (parasites, plant diet and lichen diet) among locations and between sex, Kruskal-Wallis tests (kruskal.test function) were performed, followed by Dunn post hoc test (dunn.test function) when significant. The same approach was used to assess differences in the relative abundance of various taxonomic groups (i.e. eukaryotic, fungal and bacterial phyla, plant families, and the twenty most abundant genera of each dataset) among locations and between sex. Correlations between the abundance of the twenty most abundant bacterial genera and diet and parasite richness were also tested using the Spearman's Rank correlation (cor.mtest function).

The effects of sampling locations and sex on beta diversity for each marker gene were tested by permutational multivariate analysis of variance (PERMANOVA) using the “adonis” function of the “vegan” package, based on Bray-Curtis dissimilarity matrices obtained from the “vegdist” function. When significant, differences between sampling locations were tested by pairwise PERMANOVA using the pairwise.adonis function. Non-metric multidimensional scaling (NMDS) plots were generated for each taxonomic marker using the “metaMDS” function, based on Bray-Curtis dissimilarity to visualize differences among locations and between sex. The “envfit” function of the “vegan” package was used to test whether the parasite, plant and lichen richness were correlated to the bacterial NMDS space. Only significant correlations were mapped on the NMDS plots.

18S rRNA gene

As expected, sequencing of the 18S rRNA gene allowed for the detection of all three Eukaryote kingdoms. A total of 44 phyla were detected across the seven locations. The Phragmoplastophyta phylum had the highest relative abundance (21.3%), followed by the Archamoebae (17.4%), Chlorophyta (11.2%) and Nematoda (10.9%) phyla (Fig. 2a). Significant differences in the relative abundance of several phyla among locations and between sex were detected (Supplementary Table S3A). For example, the relative abundance of Neocallimastigomycota was lowest for Central Saskatchewan and Wabowden and highest for Charron Lake and La Ronge. The relative abundance of Nematoda was higher in males, while the relative abundance of Archamoebae was higher in females.

Eukaryotic richness (Chao1 index) and diversity (Shannon index) were generally higher in La Ronge when compared to the other locations (Supplementary Table S3B). The composition of the eukaryotic community also differed among locations (Fig. 3a) and the variance in beta diversity was mainly explained by this factor (R² = 0.205, p = 0.001) and, to a lesser extent, by sex (R² = 0.034, p = 0.001). All Adonis pairwise comparisons between locations were significant, except for Flin Flon, which was only different from Central Saskatchewan, and Naosap-Reed and Wabowden, which were not significantly different (Supplementary Table S3C).

Five genera that host potential caribou parasites were detected in the 18S rRNA gene dataset: Eimeria (Apicomplexa), Nematodirella and Parelaphostrongylus (Nematoda), Entamoeba (Archamoebae), and Blastocystis (Stramenopila lineage). Entamoeba was the most frequently detected and was found in all locations, followed by Parelaphostrongylus and Blastocystis (Fig. 4). Eimeria and Nematodirella were only detected in two and three locations, respectively (Fig. 4). Differences in the number of ASVs assigned to parasitic genera among locations were also detected: Naosap-Reed had more Entomoeba ASVs and, overall, more ASVs assigned to parasitic genera than other locations, Charron Lake had more Nematodirellla ASVs, while Flin Flon and Naosap-Reed had more Blastocystis ASVs (Fig. 5a, Supplementary Table S3D). Interestingly, Nematodes (Parelaphostrongylus and Nematodirella) were more frequently found in males than in females (Fig. 4) and the number of ASVs assigned to these genera was also higher in males than in females (Supplementary Tables S3D). However, no difference was found between males and females for the total parasite richness.

As mentioned above, the sequencing of the 18S rRNA gene also allowed for the detection of Neocallimastigomycota fungi at relatively high abundances. The number of Neocallimastigomycota ASVs was higher at Naosap-Reed and The Bog and lower at Central Saskatchewan (Supplementary Table S3E).

Although the 18S rRNA gene is not commonly used to detect and describe plant and fungal communities, as more specific markers are available and provide better taxonomic resolution for these groups, a total of ten plant genera and nine lichen-forming algae genera were detected, providing some information about caribou’s winter diet (Supplementary Fig. S1). However, only one ASV assigned to a lichen-forming fungus was detected, while ASVs assigned to Asterochloris, the algal symbiont of Cladonia, were found at low relative abundances and were therefore not identified as part of the diet by our workflow.

2. Plant ITS2 region

A total of 45 plant families were detected based on the sequencing of the ITS2 region, with Ericaceae as the most abundant (39.28%), followed by Cyperaceae (16.84%), Sphagnaceae (7.24%), Hylocomiaceae (5.13%) and Betulaceae (4.74%) (Supplementary Fig. S2). Another 8.57% of the sequences were not assigned at the family level. The majority of the twenty most abundant plant families showed differences in their relative abundances among locations (Supplementary Table S4A). As an example, the relative abundance of Ericaceae was much higher at The Bog and Central Saskatchewan than in other locations.

No difference in plant richness was found among locations (Supplementary Table S4B). However, the Shannon diversity index was significantly higher in Naosap-Reed than in Central Saskatchewan. No difference in plant alpha diversity was detected between sex.

As for the 18S rRNA gene, beta diversity based on ITS2 sequencing was significantly influenced by the sampling location (Fig. 3b, R² = 0.237, p = 0.001). All Adonis pairwise comparisons between locations were significant, except for Naosap-Reed and Flin Flon which were similar (Supplementary Table S4C).

To describe the plant part of the diet, the non-rarefied plant community matrix was used. Of the 1,599 plant ASVs detected, 148 ASVs assigned to 30 genera passed the 1% threshold and were kept for diet analysis. Given the high number of ASVs, the FOO was calculated for each plant genus instead of each ASV, which resulted in three plant groups (Fig. 6). The first group consisted of the most frequently detected plant genera: Andromeda, Kalmia, Carex, Pleurozium and Sphagnum, which were detected in almost all, if not all, locations. The second group was composed of Rhododendron, Betula and Salix and could be defined as a “secondary diet” for which the FOO is similar for all locations (Fig. 6). The third group was composed of plant genera that were not frequently detected in caribou fecal DNA (Alnus, Arctostaphylos, Calliergon, Chamaedaphne, Comarum, Eriophorum, Gaylussacia, Juncus, Juniperus, Larix, Menyanthes, Picea, Pinus, Polytrichum, Populus, Potentilla, Ptilidum, Rubus, Sarmentypnum, Sarracenia, Triglochin, Vaccinium) or that were frequently detected in a few locations, such as Populus and Polytrichum in La Ronge, Sarracenia in The Bog and Potentilla in The Bog and Naosap-Reed (Fig. 6).

While the total number of plant ASVs composing the diet was not significantly different among locations, plant diet was composed of more ASVs assigned to lignified plants in Central Saskatchewan and The Bog, non-vascular plants in La Ronge and Charron Lake, and herbaceous plants in Naosap-Reed and Flin Flon (Fig. 5b, Supplementary Table S4D). No difference was found between sex.

3. Fungal ITS2 region

A total of nine fungal phyla were detected, with Ascomycota as the most abundant (65.24%), followed by Basidiomycota (16.91%) and Neocallimastigomycota (7.57%) (Supplementary Fig. S3). Unidentified fungi represented 9% of the total number of ASVs. Some differences in relative abundances among locations were detected (Supplementary Table S5A). As for the 18S rRNA gene, sequencing of the fungal ITS2 region showed a lower relative abundance of Neocallimastigomycota in Central Saskatchewan than in other locations.

Fungal richness (Chao1) was higher in Central Saskatchewan than in Naosap-Reed, The Bog and Wabowden (Supplementary Table S5B) and higher for females than for males. Shannon diversity index was not significantly different among locations but was higher for females than for males (Supplementary Table S5B).

The fungal beta diversity was primarily driven by location (R² = 0.198, p = 0.001), with Central Saskatchewan samples clustering separately from the others (Fig. 3c). Adonis pairwise comparisons showed that all locations were significantly different from each other, except for Flin Flon which was similar to Wabowden, Naosap-Reed and The Bog (Supplementary Table S5C), a trend that can be visualized in Fig. 3c.

As for plants and parasites, the non-rarefied fungal community matrix, which contained a total of 603 lichen ASVs, was used to describe the diet. A total of 382 lichen ASVs, assigned to 24 genera, passed the 1% threshold and were kept for diet analysis. Given the high number of ASVs, the FOO was calculated for each genus instead of ASV, which resulted in three lichen groups (Fig. 7). The first group consisted of three genera that were frequently detected in all locations: Cladonia, Bryoria and Usnea. The second group was composed of Evernia, Hypogymnia and Parmelia, epiphytic lichens more frequently detected in Central Saskatchewan (Fig. 7). The third group was composed of lichen genera detected at a low frequency in most locations (Buellia, Flavopunctelia, Imshaugia, Lecanora, Micarea, Parmeliopsis, Pectenia, Phylliscum, Porpida, Pseudevernia, Pycnora, Raesaeneniana, Ramalina, Sarea, Scoliciosporum, Stereaucolon, Violella and Xanthoparmelia), except for Lecanora in Central Saskatchewan and Raesaeneniana and Stereaucolon in La Ronge (Fig. 7).

The total number of ASVs assigned to lichens was lower in The Bog and higher in Central Saskatchewan and La Ronge. The number of ASVs for crustose, epiphyte-foliose and epiphyte-fruticose lichens was also higher in Central Saskatchewan whereas the number of ASVs for fruticose lichens was lower (Fig. 5c, Supplementary Table S5D).

4. 16S rRNA gene

A total of 39 bacterial phyla were detected, with Firmicutes as the most abundant (72.38%), followed by Bacteroidota (23.59%), Proteobacteria (0.94%), Verrucomicrobiota (0.94%) and Cyanobacteria (0.58%) (Supplementary Fig. S4). Some differences in the relative abundance of the twenty most abundant phyla were detected among locations (Supplementary Table S6A). At the genus level, the most abundant taxa were Ruminococcaceae UCG-005 (19.95%), Unassigned Firmicutes (15.47%), Bacteroides (11.8%), Ruminococcaceae UCG-010 (9.26%) and Christensenellaceae R-7 group (5.26%) (Fig. 8). The Central Saskatchewan and La Ronge locations were the more distinct, with significant differences in relative abundance detected for 8 out of the 20 top bacterial genera (Supplementary Table S6B). Of note, the relative abundance of Prevotella was much lower in Central Saskatchewan compared to other locations, but higher for Christensenellaceae R7 (Fig. 2b, Supplementary Table S6B). The richness (Chao1) and the Shannon diversity index were slightly lower in La Ronge and Central Saskatchewan than in other locations (Supplementary Table S6C).

As for other markers, the variance of beta diversity of the 16S rRNA gene was explained by location (R²=0.218, p = 0.001) while sex was nearly significant (R²=0.012, p = 0.076) (Fig. 3d). All Adonis pairwise comparisons between locations were significant, except for Flin Flon which was similar to The Bog and Naosap-Reed (Supplementary Table S6D). Analysis of correlations between the NMDS space of the bacterial community and diet and parasite richness variables showed that the distinct gut bacterial community observed in Central Saskatchewan was linked to the dissimilar diet observed in this location (higher diversity of epiphyte lichens). The richness of Nematodirella and Entamoeba was also linked to the structure of the gut microbial community (Fig. 3d, Supplementary Table S6E).

ASV richness of some plant and lichen groups was significantly correlated to differences in the abundance of some bacterial genera (Fig. 8). Of note, the diversity of epiphytic lichens was positively correlated to the relative abundance of the UCG 005 and Christensenellaceae R7 group but negatively correlated to the relative abundance of Prevotella and UCG 010 (Firmicutes, Oscillospirales class). Some significant correlations were also detected between the richness of parasites and the relative abundance of some bacterial genera (Fig. 8).

This study showed that multi-marker metabarcoding of DNA isolated from non-invasively collected fecal pellets has the potential to generate new insight into the biology and ecology of at-risk species such as caribou. Our optimized workflow is based on four taxonomic markers with a broad coverage: the 16S rRNA gene for prokaryotes, 18S rRNA gene for eukaryotes and ITS2 regions of plants and fungi. This resulted in a broad characterization of the winter diet (including both plants and lichens) and of the gut microbiota, as well as the detection of genera that include caribou parasite species.

Sequencing of the 18S rRNA gene for the identification of various functional groups (i.e., parasites, plants, lichens) from diverse phyla present in ruminant feces had, to our knowledge, never been attempted in previous studies. Indeed, for most animal species, the diet is restricted to a few groups of organisms and specific marker genes are generally used⁶⁸. Here, by targeting the 18S rRNA gene, plants and lichens as well as parasites were simultaneously detected in an herbivorous species. Five parasite genera belonging to the phyla Apicomplexa, Archamoebae, Nematoda and to the Stramenopila lineage were detected in caribou feces using this approach. Nematodes belonging to the genera Nematodirella and Parelaphostrongylus were more frequently observed and more diverse in males than in females. The trend was not tested independently for each location, given the limited number of individuals per location, but the pattern seems to be present in at least three locations (Central Saskatchewan, Charron Lake and Flin Flon) and would deserve further investigation. In mammals, it is recognized that males are typically more parasitized than females due to their behavior and biology, a fact that could be linked to male-biased mortality⁷³. Unfortunately, species-level identification of parasites was not possible based on the sequencing of the V4 region of the 18S rRNA gene but would have been of great interest for some genera. For example, the Parelaphostrongylus genus includes at least three species identified as caribou parasites (P. odocoilei, P. tenuis, and P. andersoni) but their frequency of occurrence and impact on caribou health is quite different: P. andersoni, a muscle worm infecting caribou and other Cervidae in Canada, is common in boreal caribou, whereas P. tenuis, the meningeal worm, is rare but lethal²⁵. Other complementary approaches, such as the sequencing of the ITS2 region using primers specifically designed to target nematodes³³ or the use of qPCR detection assays allowing to differentiate species, as proposed for other nematodes (see Reslova⁷⁴ et al. as an example), could be considered in this case. Several common parasitic nematodes such as Varestrongylus eleguniensis⁷⁵, Ostertagia gruehneri²⁸ or Platyhelminte species such as Moniezia benedeni⁷⁶ were expected, given the relatively high prevalence reported in above-mentioned studies, but went undetected in our fecal samples. This could be linked to the time of sampling, as parasitic infections are known to be seasonal in many species including caribou⁷⁶. Here, sample collection was exclusively done during the winter to ensure integrity of the fecal DNA in the context of caribou genetic studies.

Although this seasonal bias could be a limitation for the detection of some parasites by our workflow, the fact that sequencing of the 18S rRNA gene allowed to detect some known genera of caribou parasites, but also genera that were not previously reported, indicate that this approach could contribute to the detection of emerging parasites and be useful in the context of a biosurveillance strategy. Of note, the genus Entamoeba, which include parasite species that could be pathogenic, such as Entamoeba hystolytica in human⁷⁷, or with a level of pathogenicity which is still unresolved, such as Entamoeba bovis in ruminants⁷⁸, was observed at high frequencies in caribou feces from all locations. While no study to date has reported Entamoeba infections in caribou, new hosts are frequently reported⁷⁹. Another new potential parasite detected in caribou feces is Blastocystis, found in six of the seven locations. Reindeer (Rangifer tarandus) in northeastern China as well as other wild ungulates, including white-tailed deer (Odocoileus virginianus) in the USA, have also been identified as hosts of this ubiquitous protist^80–83. Further studies will be needed to identify the Entamoeba species colonizing caribou and assess the impact of both Entamoeba and Blastocystis, if any, on caribou health.

While useful to provide information about potential parasites at the individual and location levels, sequencing of the 18S rRNA gene was found to be less informative than more specific taxonomic markers for the characterization of the caribou winter diet. Among the fungal ASVs detected, only one was associated with a lichen-forming fungus. Lichens (and probably more so the fungal symbiont) are highly digestible for caribou¹⁵ and it is therefore possible that the fungal lichen DNA remaining in feces was not abundant enough (when compared to DNA from other eukaryotes) to be detected at the sequencing depth applied in this study. Alternatively, primers bias, a common issue in metabarcoding studies⁸⁴, may have favored the amplification of the 18S rRNA gene of other groups of eukaryotes. Nevertheless, several genera of lichen-forming algae were detected, but the information provided was less precise, as the same algal species can form various lichens, and the lichen taxonomy is therefore based on the taxonomy of the fungal symbiont⁵¹. Moreover, Asterochloris, the algal symbiont of Cladonia, was not identified as part of the diet, which is not consistent with results from previous studies^10,16 and from our fungal ITS2 analysis. At last, description of the plant diet based on the 18S rRNA gene was also not as informative as the plant ITS2 region, with a much lower number of plant genera identified. Overall, while the reasons behind these discrepancies among taxonomic markers still need to be identified, this observation further highlights the importance of using stable and consistent methods to assess the diet and infer ecological links based on metabarcoding data.

The ITS2 region is more commonly used as a plant taxonomic marker than the 18S rRNA gene, particularly to describe the plant diet, and was also found more suitable in this study, allowing for the detection of a rich winter plant diet in boreal caribou. The five most frequently detected plant genera are all growing as groundcover (Carex, Pleurozium and Sphagnum) or in the understory forest layer (Andromeda and Kalmia). The presence of mosses and Ericaceae in the diet have already been described using metabarcoding of fecal DNA in migratory caribou¹⁷ and in several other studies assessing the caribou diet using histological methods, as reviewed by Webber et al¹⁰. However, Thompson et al.¹⁶ observed, using collar cameras, that the proportion of mosses and Ericaceae in the boreal caribou’s winter diet was much lower than that of lichens, an information that can’t be obtained from metabarcoding data. The secondary food group was composed of lignified plant genera from the understory and midstory layers, known to be more foraged during the summer¹⁰. The sequencing of the plant ITS2 region made it possible to distinguish samples from different locations, especially La Ronge, where ericaceous shrubs (Kalmia and Andromeda) were less frequently detected as part of the diet while the diversity of non-vascular plants composing the diet was higher. This shift in the plant diet may be associated with the transition between boreal plain and boreal shield forest at this location⁸⁵, indicating that boreal caribou adapts its plant diet to what is present in its environment¹⁰.

Interestingly, while La Ronge was the most distinct location based on the sequencing of the plant ITS2 region, Central Saskatchewan was the most distinct when considering the fungal ITS2. The large difference observed in fungal community composition at this location was driven by a lower relative abundance and diversity of Neocallimastigomycota as well as a distinct lichen diet. Neocallimastigomycota fungi have been linked to fibrous plant digestibility and lignin degradation⁸⁶, and diet was shown to be a key factor influencing these fungal populations⁸⁷. Their low abundance within the gut microbial community of caribou in Central Saskatchewan is interesting as these samples had the most diverse lignified plant diet and a higher relative abundance of Ericaceae. On the one hand, this could be linked to a digestive bias, as less digestible food items have more chances to be detected in feces⁶⁸, while on the other hand, the lack of Neocallimastigomycota fungi in the gut microbiota can lead to a decrease in fiber digestion potential, as demonstrated in sheep⁸⁸, which could have further improved the detection of these less digestible food items in feces from Central Saskatchewan.

Sequencing of the ITS2 region provided great insight into the caribou lichen diet based on fecal DNA, which had not been previously described for caribou using a metabarcoding approach due to methodology failures^17,18. Surprisingly, despite the use of a conservative threshold of 1% to include ASVs into the diet, 24 lichen genera were detected across all samples, which is higher than the 13 lichen genera reported by Webber et al.¹⁰ across 30 different studies based on the histological analysis of caribou feces or rumen content, indicating that DNA-based approaches are likely more sensitive. As expected, Cladonia, Usnea and Bryoria were the most frequently detected lichens as they are commonly reported as important food items for caribou^10,16. The FOO for lichen genera was similar between locations, except for Central Saskatchewan, where epiphytic lichens were detected more frequently and with a higher diversity than in other locations, whereas soil fruticose lichens were less frequently detected and less diverse. As the diversity of the fruticose class is mostly driven by Cladonia in our dataset, it can be assumed that caribou from Central Saskatchewan had limited access to this main food source during the winter¹⁶. Two hypotheses could be put forward to explain why caribou from this location would forage more than others on epiphytic lichens and plants composing the midstory vegetation layer. First, caribou dig less on terrestrial lichens as the snow thickness and density increase during the winter⁸⁹ or in the presence of ice crust that prevent caribou from grazing on the ground, as recently demonstrated in Peary caribou⁹⁰. Second, the region is under intensive forest management, an activity known to affect the caribou’s winter diet^16,91. Although further studies will be needed to identify drivers of caribou lichen diet in this area, our findings clearly demonstrate the importance of considering lichens when studying the winter diet of caribou, something that our workflow allows to successfully do for the first time using a DNA-based approach.

Sequencing of the 16S rRNA gene is frequently used to describe the gut microbiota of animals. As previously described in reindeer^31,92,93, the gut microbiota of caribou was mainly composed of Firmicutes and Bacteroidota in our study, and several of the most abundant taxa detected were also found as dominant in other ruminants worldwide⁹⁴, further supporting that a core microbiota exists in ruminants. Interesting correlations between the gut microbiota and the diet have been detected, confirming trends observed in other ruminants^19,94,95. First, the beta diversity of the bacterial community was correlated with the diversity of the lichen diet and, to a lesser extent, to the herbaceous plant diversity. Alpha diversity of the lichen diet was also better correlated than the plant diet to the relative abundance of several bacterial taxa known for their functional role in the metabolism and fermentation process in ruminants^96,97, with both positive (e.g., Christensenellaceae R7 Group) and negative (e.g., Prevotella) correlations detected, especially for the epiphytic lichens. These important bacterial taxa of the caribou digestive tract could potentially be proposed as indicators of the quality of the diet, and eventually of the habitat, but more research is needed on this front. To a lesser extent, the gut microbial community was also linked to the diversity of parasites. As indicated by studies performed on livestock (including ruminants) in the context of helminth infections, reviewed by William et al.⁹⁸, while both diet and parasite infections can alter the gut microbiota, diet-induced modifications can also affect the host susceptibility to parasites. The workflow proposed here, which allows simultaneous characterization of the diet, parasites, and gut microbiota from fecal DNA, could greatly contribute to the study of the complex interactions between these various factors, and eventually lead to the identification of robust individual and population-level health indicators.

Overall, through the optimization of the various steps of the workflow, including the development of a tool allowing to simplify and standardize the bioinformatic analyses of sequences and the enrichment of the fungal ITS2 database with sequences from lichenous fungi, this work has allowed for the implementation of an efficient, non-invasive, and standardized multi-marker metabarcoding workflow for the characterization of caribou fecal DNA. The application of this workflow to samples from seven locations in Manitoba and Saskatchewan highlighted its capacity to reveal similarities and differences between locations, especially for the winter diet, as well as to establish associations between datasets, such as those detected between the lichen diet and the gut microbiota. These results showed that the development of diet and health indicators derived from the metabarcoding of caribou fecal DNA, in parallel to host genetic data, is a promising approach to gather relevant information in the context of population dynamics studies, including niche partitioning and effects of natural and anthropogenic disturbances and habitat loss assessment, therefore providing critical information to adequately protect boreal caribou and its habitat.

COMPETING INTERESTS

The authors declare no competing interests.

Author Contribution

C.M. and M.M. conceived the study. M.M, P.W., C.M. and A.D. acquired funding. G.P., M.M. and P.W. provided samples and metadata. P.G. developed a data analysis tool. M.N., MJ.B. and P.G. collected data and conducted analyses under C.M. supervision. All co-authors contributed to data interpretation. M.N. wrote the manuscript. All other co-authors reviewed the manuscript.

Acknowledgement

The authors thank Bridget Redquest for preparing and shipping the samples, as well as Caroline Bourdon, Julianne Mathon-Dufour, Camille Bernier, and Marie-Josée Morency for technical support in the laboratory. Funding from this research was provided by the Genomic Applications Partnership Program of Genome Canada and the Cumulative Effect program of the Canadian Forest Service.

Data Availability

Raw sequence reads as well as sample metadata are deposited in the Sequence Read Archive (BioProject PRJNA904780, accession numbers SRR23356982 to SRR23357367 and OQ792092 to OQ792112).

Arneth, A. et al. Post-2020 biodiversity targets need to embrace climate change. Proc. Natl. Acad. Sci. 117, 30882–30891 (2020).
Urban, M. C. Accelerating extinction risk from climate change. Science 348, 571–573 (2015).
Vors, L. S. & Boyce, M. S. Global declines of caribou and reindeer. Glob. Change Biol. 15, 2626–2633 (2009).
Environment and Climate Change Canada. Amended Recovery Strategy for the Woodland Caribou (Rangifer Tarandus Caribou), Boreal Population, in Canada. Species at Risk Act Recovery Strategy Series. 143 (2020).
Festa-Bianchet, M., Ray, J. C., Boutin, S., Côté, S. D. & Gunn, A. Conservation of caribou (Rangifer tarandus) in Canada: an uncertain future. Can. J. Zool. 89, 419–434 (2011).
Stewart, F. E. C. et al. Boreal caribou can coexist with natural but not industrial disturbances. J. Wildl. Manag. 84, 1435–1444 (2020).
Neilson, E. W. et al. The direct and habitat-mediated influence of climate on the biogeography of boreal caribou in Canada. Clim. Change Ecol. 100052 (2022) doi:10.1016/j.ecochg.2022.100052.
Serrouya, R., McLellan, B. N., Boutin, S., Seip, D. R. & Nielsen, S. E. Developing a population target for an overabundant ungulate for ecosystem restoration. J. Appl. Ecol. 48, 935–942 (2011).
Wittmer, H. U., Mclellan, B. N., Serrouya, R. & Apps, C. D. Changes in landscape composition influence the decline of a threatened woodland caribou population. J. Anim. Ecol. 76, 568–579 (2007).
Webber, Q. M. R., Ferraro, K. M., Hendrix, J. G. & Vander Wal, E. What do caribou eat? A review of the literature on caribou diet. Can. J. Zool. 100, 197–207 (2022).
Leblond, M., Dussault, C. & St-Laurent, M.-H. Development and validation of an expert-based habitat suitability model to support boreal caribou conservation. Biol. Conserv. 177, 100–108 (2014).
Rettie, W. J. & Messier, F. Hierarchical habitat selection by woodland caribou: its relationship to limiting factors. Ecography 23, 466–478 (2008).
Payette, S. & Delwaide, A. Tamm review: The North-American lichen woodland. For. Ecol. Manag. 417, 167–183 (2018).
Price, D. T. et al. Anticipating the consequences of climate change for Canada’s boreal forest ecosystems. Environ. Rev. 21, 322–365 (2013).
Parker, K. L., Barboza, P. S. & Stephenson, T. R. Protein conservation in female caribou (Rangifer tarandus): effects of decreasing diet quality during winter. J. Mammal. 86, 610–622 (2005).
Thompson, I. D. et al. Factors influencing the seasonal diet selection by woodland caribou (Rangifer tarandus tarandus) in boreal forests in Ontario. Can. J. Zool. 93, 87–98 (2015).
Brodeur, A., Leblond, M., Brodeur, V., Taillon, J. & Côté, S. D. Investigating potential for competition between migratory caribou and introduced muskoxen. J. Wildl. Manag. 87, e22366 (2023).
Mitchell, G. et al. DNA metabarcoding of faecal pellets reveals high consumption of yew (Taxus spp.) by caribou (Rangifer tarandus) in a lichen-poor environment. FACETS 7, 701–717 (2022).
Salgado-Flores, A. et al. Rumen and cecum microbiomes in reindeer (Rangifer tarandus tarandus) are changed in response to a lichen diet and may affect enteric methane emissions. PLOS ONE 11, e0155213 (2016).
Desselberger, U. The mammalian intestinal microbiome: composition, interaction with the immune system, significance for vaccine efficacy, and potential for disease therapy. Pathogens 7, 57 (2018).
Kinross, J. M., Darzi, A. W. & Nicholson, J. K. Gut microbiome-host interactions in health and disease. Genome Med. 3, 14 (2011).
Moeller, A. H. & Sanders, J. G. Roles of the gut microbiota in the adaptive evolution of mammalian species. Philos. Trans. R. Soc. B Biol. Sci. 375, 20190597 (2020).
Wolf, J. F. et al. Gut microbiome composition predicts summer core range size in two divergent ungulates. FEMS Microbiol. Ecol. 97, fiab048 (2021).
Laaksonen, S. et al. Climate change promotes the emergence of serious disease outbreaks of filarioid nematodes. EcoHealth 7, 7–13 (2010).
Kutz, S. J. et al. Chapter 2 - Parasites in ungulates of arctic North America and Greenland: A view of contemporary diversity, ecology, and impact in a world under change. in Advances in Parasitology (eds. Rollinson, D. & Hay, S. I.) vol. 79 99–252 (Academic Press, 2012).
Morales-Castilla, I. et al. Forecasting parasite sharing under climate change. Philos. Trans. R. Soc. B Biol. Sci. 376, 20200360 (2021).
Albon, S. D. et al. The role of parasites in the dynamics of a reindeer population. Proc. R. Soc. Lond. B Biol. Sci. 269, 1625–1632 (2002).
Dickinson, E. R., Orsel, K., Cuyler, C. & Kutz, S. J. Life history matters: Differential effects of abomasal parasites on caribou fitness. Int. J. Parasitol. 53, 221–231 (2023).
Irvine, R. J., Corbishley, H., Pilkington, J. G. & Albon, S. D. Low-level parasitic worm burdens may reduce body condition in free-ranging red deer (Cervus elaphus). Parasitology 133, 465–475 (2006).
Beldomenico, P. M. & Begon, M. Disease spread, susceptibility and infection intensity: vicious circles? Trends Ecol. Evol. 25, 21–27 (2010).
Ilina, L. A. et al. The comparative analysis of the ruminal bacterial population in reindeer (Rangifer tarandus L.) from the Russian Arctic zone: regional and seasonal effects. Animals 11, 911 (2021).
Pope, P. B. et al. Metagenomics of the Svalbard reindeer rumen microbiome reveals abundance of polysaccharide utilization loci. PLOS ONE 7, e38571 (2012).
Avramenko, R. W. et al. Exploring the gastrointestinal “nemabiome”: deep amplicon sequencing to quantify the species composition of parasitic nematode communities. PLOS ONE 10, e0143559 (2015).
Gasser, R. B., Bott, N. J., Chilton, N. B., Hunt, P. & Beveridge, I. Toward practical, DNA-based diagnostic methods for parasitic nematodes of livestock — Bionomic and biotechnological implications. Biotechnol. Adv. 26, 325–334 (2008).
Vlčková, K. et al. Diversity of Entamoeba spp. in African great apes and humans: an insight from Illumina MiSeq high-throughput sequencing. Int. J. Parasitol. 48, 519–530 (2018).
Erickson, D. L. et al. Reconstructing a herbivore’s diet using a novel rbcL DNA mini-barcode for plants. AoB PLANTS 9, plx015 (2017).
Li, J. et al. Divergence of fecal microbiota and their associations with host phylogeny in Cervinae. Front. Microbiol. 9, (2018).
McShea, W. J. et al. Metabarcoding reveals diet diversity in an ungulate community in Thailand. Biotropica 51, 923–937 (2019).
Srivathsan, A., Ang, A., Vogler, A. P. & Meier, R. Fecal metagenomics for the simultaneous assessment of diet, parasites, and population genetics of an understudied primate. Front. Zool. 13, 17 (2016).
Christopherson, V., Tremblay, J.-P., Gagné, P. N., Bérubé, J. & St-Laurent, M.-H. Meeting caribou in the alpine: Do moose compete with caribou for food? Glob. Ecol. Conserv. 20, e00733 (2019).
O’Rourke, D. R., Bokulich, N. A., Jusino, M. A., MacManes, M. D. & Foster, J. T. A total crapshoot? Evaluating bioinformatic decisions in animal diet metabarcoding analyses. Ecol. Evol. 10, 9721–9739 (2020).
Hettinga, P. N. et al. Estimating size and trend of the North Interlake woodland caribou population using fecal-DNA and capture–recapture models. J. Wildl. Manag. 76, 1153–1164 (2012).
Ball, M. C., Finnegan, L., Manseau, M. & Wilson, P. Integrating multiple analytical approaches to spatially delineate and characterize genetic population structure: an application to boreal caribou (Rangifer tarandus caribou) in central Canada. Conserv. Genet. 11, 2131–2143 (2010).
Ball, M. C. et al. Characterization of target nuclear DNA from faeces reduces technical issues associated with the assumptions of low-quality and quantity template. Conserv. Genet. 8, 577–586 (2007).
McFarlane, S., Manseau, M. & Wilson, P. J. Spatial familial networks to infer demographic structure of wild populations. Ecol. Evol. 11, 4507–4519 (2021).
Priadka, P. et al. Partitioning drivers of spatial genetic variation for a continuously distributed population of boreal caribou: Implications for management unit delineation. Ecol. Evol. 9, 141–153 (2019).
Parada, A. E., Needham, D. M. & Fuhrman, J. A. Every base matters: Assessing small subunit rRNA primers for marine microbiomes with mock communities, time series and global field samples. Environ. Microbiol. 18, 1403–1414 (2016).
Sundset, M. A., Praesteng, K. E., Cann, I. K. O., Mathiesen, S. D. & Mackie, R. I. Novel rumen bacterial diversity in two geographically separated sub-species of reindeer. Microb. Ecol. 54, 424–438 (2007).
Tryland, M. & Kutz, S. J. Reindeer and Caribou: Health and Disease. (CRC Press, 2018).
Stoeck, T. et al. Multiple marker parallel tag environmental DNA sequencing reveals a highly complex eukaryotic community in marine anoxic water. Mol. Ecol. 19, 21–31 (2010).
Lücking, R., Hodkinson, B. P. & Leavitt, S. D. The 2016 classification of lichenized fungi in the Ascomycota and Basidiomycota – Approaching one thousand genera. The Bryologist 119, 361–416 (2017).
Kelly, L. J. et al. DNA barcoding of lichenized fungi demonstrates high identification success in a floristic context. New Phytol. 191, 288–300 (2011).
White, T. J., Bruns, T., Lee, S. & Taylor, J. Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics. in PCR protocols: a guide to methods and applications 315–322 (Academic Press, 1990). doi:10.1016/B978-0-12-372180-8.50042-1.
Menkis, A. et al. Occurrence and impact of the root-rot biocontrol agent Phlebiopsis gigantea on soil fungal communities in Picea abies forests of northern Europe. FEMS Microbiol. Ecol. 81, 438–445 (2012).
Hollingsworth, P., Graham, S. & Little, D. Choosing and using a plant DNA Barcode. PloS One 6, e19254 (2011).
Banchi, E. et al. PLANiTS: a curated sequence reference dataset for plant ITS DNA metabarcoding. Database 2020, baz155 (2020).
Chen, S. et al. Validation of the ITS2 region as a novel DNA barcode for identifying medicinal plant species. PLOS ONE 5, e8613 (2010).
Bolyen, E. et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat. Biotechnol. 37, 852–857 (2019).
Nguyen, N. H. et al. FUNGuild: An open annotation tool for parsing fungal community datasets by ecological guild. Fungal Ecol. 20, 241–248 (2016).
Quast, C. et al. The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucleic Acids Res. 41, D590–D596 (2012).
Abarenkov, K. et al. The UNITE database for molecular identification of fungi – recent updates and future perspectives. New Phytol. 186, 281–285 (2010).
Marthinsen, G., Rui, S. & Timdal, E. OLICH: A reference library of DNA barcodes for Nordic lichens. Biodivers. Data J. 7, (2019).
Hill, D. J. Keys to lichens of North America: revised and expanded. The Lichenologist 48, 255–255 (2016).
Callahan, B. J. et al. DADA2: High-resolution sample inference from Illumina amplicon data. Nat. Methods 13, 581–583 (2016).
Sanders, W. B. & Masumoto, H. Lichen algae: the photosynthetic partners in lichen symbioses. The Lichenologist 53, 347–393 (2021).
Ando, H. et al. Evaluation of plant contamination in metabarcoding diet analysis of a herbivore. Sci. Rep. 8, 15563 (2018).
Pompanon, F. et al. Who is eating what: diet assessment using next generation sequencing. Mol. Ecol. 21, 1931–1950 (2012).
Ando, H. et al. Methodological trends and perspectives of animal dietary studies by noninvasive fecal DNA metabarcoding. Environ. DNA 2, 391–406 (2020).
Jiménez, P., Muñoz, M. & Ramírez, J. D. An update on the distribution of Blastocystis subtypes in the Americas. Heliyon 8, (2022).
Anamalai, S. et al. Blastocystis in captivated and free-ranging wild animals worldwide: a revie. Trop. Biomed. 39, 338–372 (2022).
Grattepanche, J.-D. et al. Microbial diversity in the Eukaryotic SAR Clade: illuminating the darkness between morphology and molecular data. BioEssays 40, 1700198 (2018).
Schloss Patrick D. Rarefaction is currently the best approach to control for uneven sequencing effort in amplicon sequence analyses. mSphere 9, e00354-23 (2024).
Moore, S. L. & Wilson, K. Parasites as a viability cost of sexual selection in natural populations of mammals. Science 297, 2015–2018 (2002).
Reslova, N. et al. The identification and semi-quantitative assessment of gastrointestinal nematodes in faecal samples using multiplex real-time PCR assays. Parasit. Vectors 14, 391 (2021).
Verocai, G. G. et al. The biogeography of the caribou lungworm, Varestrongylus eleguneniensis (Nematoda: Protostrongylidae) across northern North America. Int. J. Parasitol. Parasites Wildl. 11, 93–102 (2020).
Turgeon, G., Kutz, S. J., Lejeune, M., St-Laurent, M.-H. & Pelletier, F. Parasite prevalence, infection intensity and richness in an endangered population, the Atlantic-Gaspésie caribou. Int. J. Parasitol. Parasites Wildl. 7, 90 (2018).
Stauffer, W. & Ravdin, J. I. Entamoeba histolytica: an update. Curr. Opin. Infect. Dis. 16, 479–485 (2003).
Matsubayashi, M. et al. First detection and molecular identification of Entamoeba bovis from Japanese cattle. Parasitol. Res. 117, 339–342 (2018).
Jacob, A. S., Busby, E. J., Levy, A. D., Komm, N. & Clark, C. G. Expanding the Entamoeba universe: new hosts yield novel ribosomal lineages. J. Eukaryot. Microbiol. 63, 69–78 (2016).
Betts, E. L. et al. Genetic diversity of Blastocystis in non-primate animals. Parasitology 145, 1228–1234 (2018).
Maloney, J. G., Jang, Y., Molokin, A., George, N. S. & Santin, M. Wide genetic diversity of Blastocystis in white-tailed deer (Odocoileus virginianus) from Maryland, USA. Microorganisms 9, (2021).
Kaczmarek, A. et al. Blastocystis occurrence and subtype diversity in wild European terrestrial mammals – The case of Białowieża Primeval Forest (NE Poland). Int. J. Parasitol. Parasites Wildl. 16, 120–125 (2021).
Wang, J. et al. Distribution and genetic diversity of Blastocystis subtypes in various mammal and bird species in northeastern China. Parasit. Vectors 11, 522 (2018).
Kelly, R. P., Shelton, A. O. & Gallego, R. Understanding PCR processes to draw meaningful conclusions from environmental DNA studies. Sci. Rep. 9, 12133 (2019).
Padbury, G. A., Acton, D. F., Stushnoff, C. T. & Center, U. of R. C. P. R. The Ecoregions of Saskatchewan. (University of Regina Press, 1998).
Grenet, E., Breton, A., Barry, P. & Fonty, G. Rumen anaerobic fungi and plant substrate colonization as affected by diet composition. Anim. Feed Sci. Technol. 26, 55–70 (1989).
Gruninger, R. J. et al. Anaerobic fungi (phylum Neocallimastigomycota): advances in understanding their taxonomy, life cycle, ecology, role and biotechnological potential. FEMS Microbiol. Ecol. 90, 1–17 (2014).
Gao, A. W., Wang, H. R., Yang, J. L. & Shi, C. X. The Effects of elimination of fungi on microbial population and fiber degradation in sheep rumen. Appl. Mech. Mater. 295–298, 224–231 (2013).
Thomas, D. C., Edmonds, E. J. & Brown, W. K. The diet of woodland caribou populations in west-central Alberta. Rangifer 337–342 (1996) doi:10.7557/2.16.4.1275.
Gautier, C., Langlois, A., Sasseville, V., Neave, E. & Johnson, C. A. Remote sensing, snow modelling, survey data and Indigenous Knowledge show how snow and sea-ice conditions affect Peary caribou (Rangifer tarandus pearyi) distribution and inter-island and island–mainland movements. Polar Res. 41, (2022).
Kivinen, S., Moen, J., Berg, A. & Eriksson, Å. Effects of modern forest management on winter grazing resources for reindeer in Sweden. AMBIO 39, 269–278 (2010).
Glendinning, L., Genç, B., Wallace, R. J. & Watson, M. Metagenomic analysis of the cow, sheep, reindeer and red deer rumen. Sci. Rep. 11, 1990 (2021).
Zielińska, S., Kidawa, D., Stempniewicz, L., Łoś, M. & Łoś, J. M. New insights into the microbiota of the Svalbard reindeer Rangifer tarandus platyrhynchus. Front. Microbiol. 7, (2016).
Henderson, G. et al. Rumen microbial community composition varies with diet and host, but a core microbiome is found across a wide geographical range. Sci. Rep. 5, 14567 (2015).
Wu, X. et al. Characterizing the alteration in rumen microbiome and carbohydrate-active enzymes profile with forage of muskoxen rumen through comparative metatranscriptomics. Microorganisms 10, 71 (2022).
Perea, K. et al. Feed efficiency phenotypes in lambs involve changes in ruminal, colonic, and small-intestine-located microbiota1. J. Anim. Sci. 95, 2585–2592 (2017).
Patra, A. K. & Yu, Z. Genomic insights into the distribution of peptidases and proteolytic capacity among Prevotella and Paraprevotella species. Microbiol. Spectr. 10, e02185-21 (2022).
Williams, A. R. et al. Emerging interactions between diet, gastrointestinal helminth infection, and the gut microbiota in livestock. BMC Vet. Res. 17, 62 (2021).

No competing interests reported.

Download PDF

Editorial decision: Revision requested
02 Jul, 2024
Reviews received at journal
28 Jun, 2024
Reviews received at journal
19 Jun, 2024
Reviewers agreed at journal
19 Jun, 2024
Reviewers agreed at journal
18 Jun, 2024
Reviewers invited by journal
18 Jun, 2024
Editor assigned by journal
10 Jun, 2024
Editor invited by journal
27 May, 2024
Submission checks completed at journal
23 May, 2024
First submitted to journal
21 May, 2024

You are reading this older preprint version

Read the latest preprint version →

Exploring winter diet, gut microbiota and parasitism in caribou using multi-marker metabarcoding of fecal DNA

Status:

Version 1

Abstract

Figures

INTRODUCTION

MATERIALS AND METHODS

RESULTS

DISCUSSION

Declarations

COMPETING INTERESTS

Author Contribution

Acknowledgement

Data Availability

References

Additional Declarations

Supplementary Files

Status:

Version 1