Metabarcoding in two isolated populations of wild roe deer (Capreolus capreolus) reveals variation in gastrointestinal nematode community composition between regions and among age classes

Background Gastrointestinal nematodes are ubiquitous for both domestic and wild ungulates and have varying consequences for health and fitness. They exist as complex communities of multiple co-infecting species, and we have a limited understanding of how these communities vary in different hosts, regions and circumstances or of how this affects their impacts. Methods We have undertaken ITS2 rDNA nemabiome metabarcoding with next-generation sequencing on populations of nematode larvae isolated from 149 fecal samples of roe deer of different sex and age classes in the two isolated populations of Chizé and Trois Fontaines in France not co-grazing with any domestic ungulate species. Results We identified 100 amplified sequence variants (ASVs) that were assigned to 14 gastrointestinal nematode taxa overall at either genus (29%) or species (71%) level. These taxa were dominated by parasites classically found in cervids—e.g. Ostertagia leptospicularis, Spiculopteragia spp. Higher parasite species diversity was present in the Trois Fontaines population than in the Chizé population including the presence of species more typically seen in domestic livestock (Haemonchus contortus, Bunostomum sp., Cooperia punctata, Teladorsagia circumcincta). No differences in parasite species diversity or community composition were seen in the samples collected from three zones of differing habitat quality within the Chizé study area. Young roe deer hosted the highest diversity of gastrointestinal nematodes, with more pronounced effects of age apparent in Trois Fontaines. The effect of host age differed between gastrointestinal nematode species, e.g. there was little effect on O. leptospicularis but a large effect on Trichostrongylus spp. No effect of host sex was detected in either site. Conclusions The presence of some livestock parasite species in the Trois Fontaines roe deer population was unexpected given the isolation of this population away from grazing domestic livestock since decades. Overall, our results illustrate the influence of host traits and the local environment on roe deer nemabiome and demonstrate the power of the nemabiome metabarcoding approach to elucidate the composition of gastrointestinal nematode communities in wildlife. Graphical Abstract Supplementary Information The online version contains supplementary material available at 10.1186/s13071-021-05087-5.

Beaumelle et al. Parasites & Vectors (2021) 14:594 Background Gastrointestinal nematode parasites of wildlife and livestock have a worldwide distribution. Multiple species usually co-infect a single host, and they have major ecological, economic and animal welfare impacts [1][2][3]. Gastrointestinal nematodes of wildlife can have negative impacts on the health and individual fitness of hosts [4] as well as on their population dynamics [5], which, in extreme cases, can lead to population extinction [6]. Consequently, understanding host-parasite relationships and determining parasite diversity and abundance are crucial for sustainable wildlife management [7,8].
The immunological response can help the host to control parasite infestations but it incurs a significant cost to its resources and must be balanced with other energy demanding functions [4,9], especially when the quality of habitat is poor [10]. Therefore, the physiological and behavioral attributes of individual hosts influence infection level, parasite community composition and their effects [11,12]. The age, sex and reproductive status of hosts are important factors in this context and determine the parasite community structure, i.e. the abundance and diversity of parasites. As their immune system is still immature, the highest abundance and diversity of parasites are generally observed in young individuals (e.g. roe deer, Capreolus capreolus [13]; Dall's sheep, Ovis dalli dalli [14]; feral sheep, Ovis aries [3,15]). Similarly, because of the deterioration of their immune system, senescent individuals usually host a higher abundance and diversity of parasites than younger adults [16,17]. Variation in parasite abundance and diversity is also observed between sexes in various species, including ungulates (e.g. roe deer [13], African buffalo, Syncerus caffer [18], mouflon, Ovis gmelini musimon × Ovis sp. [19], feral sheep [3]), and is generally higher in males [20]. The higher susceptibility of males to infection has been ascribed to the immunosuppressive effects of sex steroid hormones in males [21], sexual size dimorphism [22] and also increased likelihood of pathogen exposure [23]. In addition to the influences of life history traits, host populations inhabiting different areas can be infected by various parasites because of the direct effect of environment on parasite communities (e.g. temperature, precipitations) [24] or via the effect of environment on hosts [25].
No relation (e.g. moose, Alces alces) [26] or negative relationship (e.g. chamois, Dall's sheep) [14,27] on the health and individual fitness of host with increased parasite richness were reported. This might be explained by the nematode community composition with different levels of pathogenicity among parasites, even among closely phylogenetically related species, and may have synergistic or antagonistic interactions [18,[28][29][30]. Thus, refining the identification of the parasite communities would help to improve our knowledge on host-parasite relationships. However, investigating parasite diversity in wildlife can be challenging [31]. Recently, next-generation sequencing metabarcoding approaches have been developed to investigate gastrointestinal nematode communities [31,32], although rarely applied to wild ungulates so far [33,34]. Metabarcoding has a number of advantages including the study of endoparasite communities, without prior assumption of which species are likely to be present. It also allows large numbers of samples to be processed simultaneously and can be highly sensitive to detect rare species [35]. Metabarcoding does not quantify the absolute number of parasites in a sample, but provides information on the proportions of each parasite species.
The aim of this study was to investigate the gastrointestinal nematode communities present in two contrasting natural populations of roe deer (Capreolus capreolus), a common and widely distributed ungulate species in Europe. The gastrointestinal nematodes of roe deer have been previously investigated in wild populations [36][37][38][39] and have been succinctly described in one of the populations studied here [13]. The two populations of roe deer examined in this study have been spatially isolated in fenced areas for more than 60 years, preventing ungulate movement. Therefore, we expected the diversity of nematodes to be dominated by roe deer-specific parasites (H1). We also expected to observe a higher parasitic diversity in young individuals of a few months of age than in adults (H2) and in males than in females (H3). Because differences in immune responses and differences of body conditions were reported between Chizé and Trois Fontaines [10], we predicted differences of nemabiome with a higher diversity of gastrointestinal nematodes in roe deer living in poor (i.e. Chizé) compared with good environmental conditions (i.e. Trois Fontaines). We predicted similar results at local scale, among the different sectors of Chizé [33], with a higher diversity in the sector with the poorest environmental conditions (H4).

Study areas
The two study sites, Chizé and Trois Fontaines, are located in western (46°05' N, 0°25' W) and northeastern France (48°43' N, 2°61' W), respectively. Both sites are fenced forests of 2614 ha in Chizé (since 1978) and 1360 ha in Trois Fontaines (since 1976), preventing roe deer and other ungulates to move into or out of the study sites. No ungulates other than roe deer are present within the two sites, except wild boar. Although some rare hare Lepus europaeus inhabit the site of Chizé, there are no lagomorphs in Trois Fontaines (Office Français de la Biodiversité, personal communication). The two sites have different climate, soil and forest attributes. The oceanic climate in Chizé is characterized by mild winters and dry summers in contrast to Trois Fontaines where the continental climate results in rainy warm summers and cold winters. In both sites, forests are dominated by oak (Quercus sp.), beech (Fagus sylvatica) and hornbeam (Carpinus betulus), but with a more heterogeneous spatial distribution of species and coppice abundance and quality in Chizé [40] compared to Trois Fontaines. In Chizé, the variation of coppice type is categorized into three sectors characterized by a rich, medium or poor environment for roe deer (Fig. 1). With its low productive forest, Chizé is still considered a poor habitat for roe deer, in contrast to the highly productive forest of Trois Fontaines [41]. In fact, roe deer in Trois Fontaines have access to high quantity and quality of food in the whole study area, while roe deer in Chizé experience a trade-off between high quantity or quality of forage [42]. Consequently, in Trois Fontaines, the population of roe deer is increasing, whereas the population growth rate is close to 0 in Chizé [10]. Both populations are kept stable by yearly removals (Office Français de la Biodiversité, personal communication).

Roe deer data
Roe deer of known age are captured annually as part of a long-term monitoring program. Briefly, roe deer were caught using long nets surrounding forest patches and beaters and maintained in wooden boxes until the end of the drive. Details of the capture procedure are available in Gaillard et al. [43]. Fecal samples were collected rectally from immobilized roe deer captured between January For each individual, we recorded the sex and age using tooth eruption pattern [43] as well as the approximate location of capture. We also verified that individuals of different age and sex are equally distributed throughout the capture period and spatial localization (data not shown).

Parasitological analyses and culture
Feces collected from roe deer per rectum were placed in plastic bags and air removed before sealing. In the field, bags were stored in an insulated box to avoid exposure to low temperatures, which could influence survival and development of some parasite species [44]. Samples were then either directly transported or mailed to the parasitology laboratory of VetAgro Sup (Lyon, France) arriving within 24-48 h of collection. On receipt, we set up coprocultures with 5.4 ± 3.7 g (minimum = 1 g, maximum = 23 g) of feces mixed with vermiculite. After an incubation period of 10-12 days at 24 °C with regular mixing and moisture, the larvae present were collected with a Baermann apparatus after 24 h of sedimentation in water. After determining the total number of L3 collected by counting an aliquot, we stored them in 70% ethanol until DNA extraction.
We also counted eggs of gastrointestinal nematodes following a modified McMaster protocol [45], with a solution of ZnSO 4 (density = 1.36) and a McMaster slide, providing a theoretical sensitivity of 15 eggs per gram (epg) of feces. We also prepared a "control slide" to improve sensitivity of the technique by fully filling a 14-ml tube with the remaining solution and covering it with a coverslip. After centrifugation (5 min at 1200 rpm), the coverslip was transferred on a microscope slide before microscopical observation. We attributed the value 7.5 epg for parasite species detected on the control slide but not on the McMaster slide. Coproscopic data were reported only for samples containing enough L3 collected from coproculture and analyzed by metabarcoding.

DNA extraction and metabarcoding
DNA extraction was performed using the Qiagen DNeasy ® powerSoil kit (Qiagen, Hilden, Germany) following the manufacturer's instructions with an elution volume of 50-150 µl. ITS2 rDNA metabarcoding was performed using the protocol described by Avramenko et al. [32]. Briefly, the rDNA ITS2 target was PCR-amplified from a 1:10 dilution of genomic DNA template. Negative experimental controls were included consisting of water added to the PCR reaction instead of template DNA. These controls were handled in exactly the same way as the actual samples. All samples (n = 149) and negative controls (n = 22) were tagged with unique barcode identifiers to allow the pooling (or normalization) into a single 100 ng amplicon library. The final concentration of the pooled library was assessed with the KAPA qPCR Library Quantification Kit (KAPA Biosystems, USA) following the manufacturer's recommended protocol. The prepared pooled library was run on an Illumina MiSeq Desktop Sequencer using a 600-cycle pair-end reagent kit (MiSeq Reagent Kits v3) at a concentration of 15 nM with the addition 25% PhiX Control v3 (Illumina). A standard demultiplexing protocol was employed that generated FASTQ files for all barcoded samples in the pooled library.

Bioinformatics
All analyses were carried out using R 3.6 (R Core team, 2020). We curated DNA sequence data in two steps (see Additional file 1: Figure S1). First, the sequence reads were passed through a DADA2-based pipeline [46]. Specific details are available at www. nemab iome. ca. Primers were removed with Cutadapt [47] after the removal of ambiguous bases. Next, reads were discarded based on the number of errors (> 2 and > 5 in forward and reverse reads, respectively), length (< 50 bp) and quality (truncG = 2 and phiX genome). Forward and reverse reads were merged only if they overlapped by > 12 bp with low mismatches (i.e. maximum 3 for > 100 bp, maximum 2 for > 50 bp and 1 for < 50 bp reads). Putative chimeras were automatically removed but we verified the effect of this filtering on the data set post-analysis, as these programs are known to produce many false positives [48]. A conservative approach was used for taxonomic assignment of amplicon sequence variants (ASV) in which three different assignment methods were used (Additional file 1: Table S1): (1) IDTaxa using the default conservative threshold of 60, which corresponds to the confidence at which to truncate the output taxonomic classification [49]; (2) assignTaxonomy [46,50] and (3) BLASTn [51]. The nematode ITS2 rDNA database 1.1.0 [52] containing 9811 complete and non-redundant sequences, with a few additional manual corrections of species assignation corresponding to morphotypes, was used for both IDTaxa and assignTaxonomy. For BLASTn, the three best hits were considered only if they reached a threshold of 90% similarity against NCBI database sequences (2021/04/13). All species-level assignments were ascribed a 'confidence level' (Additional file 1: Table S1): high confidence indicates that all three methods assigned to the same species and moderate confidence indicates that two out of the three methods agree at the species level and also the third method was in agreement but could only be assigned to the genus level. Taxonomic assignment was limited to the genus level when there was disagreement across any of the three methods at the species level (but consistency at the genus level) or if there was no similar sequence on GenBank with > 98% of similarity and a query cover > 85% with the corresponding ASV.
The raw ASV data set was curated following a procedure developed by Taberlet et al. (2018). ASVs observed only once in the whole data set were considered spurious and were removed [53]. To delete potential contaminants [54], we discarded any ASVs more abundant in the negative controls (maximum 55 reads) than in other samples and reads with < 0.13% (maximizing the removal of sequences with low abundance in controls) of the total ASV abundance in the entire data set. We only kept ASVs identified as gastrointestinal nematodes at least to the genus level. Finally, any samples with < 1000 reads were removed.

Phylogeny
A phylogenetic analysis of ASV was conducted using a maximum likelihood tree with the HKY + G [55] model according to "ModelTest" function [56,57]. We used FigTree 1.4.4 [58] to visualize ASV phylogeny. Although the generation of ASVs enables genetic diversity to be characterized to the highest resolution, the biological relevance of much of this diversity remains unknown. A limitation of using ASVs to characterize genetic diversity is that large sample sizes are required to detect significant statistical differences among samples (sensitive to type II error). To overcome this hurdle, we also clustered ASVs into "taxa" and OTUs (operational taxonomic units) using the "IdCluster" function of the package DECIPHER [59] and a posteriori cutoff of 0.03.

Statistical analyses on measures of parasitic community
We removed the lungworm Varestrongylus from the data set because the focus of the present study was on gastrointestinal nematodes. To control the representativeness of the data, we drew the accumulation curves for the ASV and taxa with the "exact" method of the function "specaccum" of the package vegan [60]. Instead of the exact read count per sample and items, we used the relative frequencies of reads. The read relative frequencies in samples were plotted with ggtree [61]. After controlling for the year, the capture date, the epg and the number of L3 used for DNA extraction, we assessed the effects of age (i.e. young [< 1 year] or adults) and sex (i.e. males or females) on two proxies of diversity of gastrointestinal nematodes community: alpha (measure of community diversity within individual samples) and beta (measure of comparison of dissimilarity between each pair of samples). We measured the beta diversity by using the Bray-Curtis dissimilarity matrices and the alpha diversity by computing three indices: richness (the number of different items), Shannon-Weaver index [62] and Simpson index (1-D) [63] indices.
To test for differences in nemabiome among individual hosts of different age, sex and site/sector of Chizé, we tested for the influence of these factors, and their interactions, on the alpha and beta diversity. We ensured that all samples come from a unique individual to avoid individual effect in statistical analyses. Doing so, we removed 19 samples (17 from Chizé and 2 from Trois Fontaines; Additional file 1: Figure S2), selecting for each repeated measure the samples with the higher mass of feces used for the coprocultures. Age and sex effects were first assessed at broad scale with the two sites (Chizé and Trois Fontaines) and, second, at the local scale, comparing the three sectors of Chizé (rich coppices, medium coppices and poor beech stand). We used generalized linear models in a model selection approach to test the effects of site, sex and age of individuals on the alpha diversity, with a Poisson family for richness, a Gaussian family for Shannon and a Gaussian (with logit link because Simpson index ranges from 0 to 1) family for Simpson. For the beta diversity, we used a custom function to compute Akaike's information criterion corrected for small sample size (AICc) based on residual sums of squares [64]. Following the recommendations of Burnham and Anderson [65], we identified models that are biologically meaningful and considering limits due to the sample size. It resulted in 18 competitive models. We then ranked the models for each set of candidate models using the AICc and calculated ΔAICc and AICc weights. We selected the model with the lowest AICc value. Models with ΔAICc ≤ 2 were considered equivalent [65], and in this case, we considered the most parsimonious one, i.e. the model with the lowest degree of freedom.
Differences in nemabiome among samples (beta diversity) were visualized on two dimensions with a nonmetric multidimensional scaling (NMDS). We used specifically the NMDS because Bray-Curtis dissimilarities are non-Euclidean embeddable distances [66]. We verified the consistency of patterns between samples resulting from the three taxonomic levels with Procrustes analyses.

Parasite sampling
ITS2 rDNA nemabiome analyses were performed on 149 samples (n = 36 and 113 samples in Trois Fontaines and Chizé, respectively; Table 1, Fig. 1). The year of sampling had no influence on the nemabiome (Adonis test and Kruskal-Wallis rank sum test; see Additional file 1: Table S2), and this factor was therefore not considered in the following analyses.
The median number eggs per gram of feces was low and heterogeneous (7.5 eggs/g [0; 114] 95%IQR ) as well as the number of L3 used for DNA extraction (80 [12; 993] 95%IQR ; Table 1, Fig. 2). Both had no significant effects on alpha and beta diversity (Pearson correlation test and Adonis test; see Additional file 1: Table S2).

Taxonomic and community diversity
Considering the whole data set, we detected 100 ASVs resolved across 14 taxa of gastrointestinal nematodes and 2 taxa of lungworms (Table 1 and Additional file 1: Table S1). ASVs were identified either at the genus (43%) or species (57%) level (Additional file 1: Table S1). Ten genera, including ten species, were identified, and their phylogeny was consistent with previous knowledge on nematode taxonomy (Additional file 1: Figure S3). Sampling captured most of the taxa as the taxa accumulation  Fig. 2 Egg count, number of L3 and read relative frequencies of gastrointestinal nematodes at taxa level. Each stacked bar chart represents the species composition of 149 roe deer samples within which each taxon is defined by one color. The numbers of L3 and fecal egg counts are displayed above each bar chart. The data are split based on site location (Chizé or Trois Fontaines), age (young or adult roe deer) and sex (F: females and M: males). In addition, samples from Chizé are subdivided according to the three different sectors: sector 1: "rich coppices;" sector 2: "medium coppices;" sector 3: "poor beech stand" curves rapidly reached a plateau (Additional file 1: Figure  S4). Forty-seven ASVs were found only in one sample for which we were not able to disentangle whether their origin was intra-specific genetic diversity or PCR/sequencing error, but following Callahan et al. [46], we considered those ASVs in analyses whatever their frequencies in the data set were. Procrustes analyses performed on ASV with taxa have tight correlations (P = 0.001). Thus, we only show results at the ASV taxonomic level, but results from taxa and OTUs are available in Additional file 1: Table S3-S4, Figure S5 and Additional file 2, respectively.
According to the best models on the data subset without repeated measures ( Table 2; Additional file 1: Table S3), the richness and Simpson diversity index were higher in Trois Fontaines than in Chizé (estimate ± standard error, richness: β = 0.20 ± 0.06, P = 0.002; Simpson: β = 0.24 ± 0.12, P = 0.041). The Shannon index, which considers the relative abundances of ASV, included significant interactions between age and site, indicating a higher alpha diversity in younger roe deer in Trois Fontaines (β = 0.32 ± 0.13, P = 0.017) ( Fig. 3; Table 2). At the local scale of Chizé, no differences among the three sectors were detected for all alpha diversity indexes (Fig. 2).
According to the most parsimonious model selected for the Bray-Curtis dissimilarity, the differences of nemabiome between samples were mostly explained by the site (F 1,127 = 35.92, P = 0.001) and age (F 1,127 = 12.68, P = 0.001) variables. Locally in Chizé, the selected model for Bray-Curtis dissimilarity depended only on age (F 1,94 = 6.36, P = 0.002). Similarly, we observed distinct groups in the NMDS, with partially separated clouds of points for Chizé and Trois Fontaines and a marked separation between young and adult roe deer in Trois Fontaines (Fig. 4).

Table 2 Best generalized linear models selected for ASVs and each diversity index (richness, Simpson, Shannon)
Gaussian families were used for Simpson and Shannon regression and Poisson family for richness. The effect of site (Chizé as reference), age (adult as reference) and sex (females as reference) and the interaction between them are reported when included in the selected models. Parameter estimates with SE are reported with the corresponding z-value (Poisson family) or t-value (Gaussian family) and P-value. Statistical significance is represented by *P < 0.05, **P < 0.01 and ***P < 0.  (Table 2) and the significative difference of nemabiome between the two sites ( Table 3). We also observed variation in the prevalence and relative abundance of five Trichostrongylus spp. ASVs between the two sites and between the host age groups, with greater differences of abundance between age in Trois Fontaines compared to Chizé (Fig. 5, A11), as indicated by the significant interaction between site and age on Shannon diversity index.

Discussion
ITS2 rDNA nemabiome metabarcoding provided fine resolution data on the gastrointestinal nematode parasite community composition of two populations of roe deer in France. The gastrointestinal nematode species detected in Trois Fontaines and Chizé correspond largely to species previously detected on necropsies of roe deer in the two study sites [13,67]. We also detected species not previously described in these two populations [13,67], but previously found in other roe deer populations (e.g. in Spain [37] and in the Netherlands [36]). This study provided new insights into siteand age-specific gastrointestinal nematode diversity because we showed that young roe deer have different nemabiomes compared to adults, with more pronounced effects of age in one site (i.e. Trois Fontaines).

Diversity of gastrointestinal nematodes in roe deer
In accordance with previous studies on helminths of roe deer (e.g. [37,38]), the main species detected were O. leptospicularis, Trichostrongylus spp. and Spiculopteragia spp. These gastrointestinal nematodes are commonly found in roe deer, confirming our hypothesis H1 that the dominant nematode species were roe deer-specific parasites given that these populations are enclosed in fenced forest areas. However, we found some differences between the two study sites. We have detected generalist species (i.e. gastrointestinal nematodes species found in different host species) at a relatively high frequency in Trois Fontaines but not in Chizé. In addition to C. ovina previously detected in Trois Fontaines [13], we have detected for the first time the presence of the highly pathogenic species H. contortus and B. trigonocephalum in this population. In Chizé and Trois Fontaines, we noted the presence of T. axei and O. venulosum, at low relative abundance, also considered as generalist species, and T. colubriformis (only in Chizé) and T. circumcincta, which generally infect sheep and goats [30,68]. The identifications of H. contortus and T. axei are at a "high confidence level" because these species were well represented in the nemabiome database and the three different methods of assignment provided the same species assignment. For Bunostomum sp., ASVs were only assigned at the genus level as a result of a relatively low level of identity (92% over 283 bp [query cover of 100%, e-value = 9e −110 ]) to the best BLAST hit in GenBank (Additional file 1: Table S1). To our knowledge, B. trigonocephalum of sheep and B. phlebotomum of calves are the only Bunostomum species which have been identified in roe deer in Europe (e.g. [36]). These two species are highly divergent, and considering the low level of intra-specific variation observed within B. trigonocephalum [69], it is possible that the Bunostomum sp. identified in Trois Fontaines is a previously unidentified, or cryptic, species closely related to B. trigonocephalum. For some ASVs, the taxonomic identification did not reach the species level, e.g. Trichostrongylus spp. (67% among this genus), Spiculopteragia spp. (6%) or Ostertagia spp. (67%). This may be largely due to the lack of reference sequence available for these species in the nemabiome database and ultimately in GenBank. This issue is accentuated when applying nemabiome metabarcoding to wildlife hosts compared with domestic livestock where parasite communities are better defined. For example, there is just a single ITS2 reference sequence for Trichostrongylus capricola compared to 299 ITS2 reference sequences for Haemonchus contortus in the nemabiome database (see also Aivelo and Medlar 2018 for the SILVA database). The potential presence of cryptic species and/ or hybrids of closely related species [70] may further complicate accurate species assignation. Consequently, the generation of reference sequences from morphologically validated parasites should continue to be a research priority for the scientific community as a whole.
Ashworthius sidemi was not detected in Trois Fontaines, although it has previously been isolated on roe deer carcasses during winter in this study area [13,71]. This may be due to the hypobiosis in this species [72] or to the fact that during winter females contained only non-embryonated eggs [73]. Further investigations on the seasonal nemabiome will be required to determine the complete gastrointestinal nematode species richness and variations during the year and species.

Effects of site, age and sex class on the nemabiome
In Trois Fontaines, the age of individuals had a major influence on the gastrointestinal nematode species community found in roe deer, with a higher alpha diversity (β = 0.32 ± 0.13, P = 0.017) observed in young compared to adult animals in accordance with our prediction (H2). This age effect on beta diversity is mostly driven by some specific variants of Trichostrongylus spp. in both sites and of H. contortus in Trois Fontaines that are more prevalent, and of higher relative abundance, in young roe deer compared to adults. The limited influence of age on prevalence and relative abundance observed for some of the other parasite species (see also e.g. [74]) may result from survival mechanisms of these parasite species, such as the dysregulation of the immune response of the host [75], from an already developed and efficient immune response of the host when young roe deer were captured against these parasite species (i.e. at 8-10 months of age in our study) or from an immune response that is inefficient. The higher prevalence and relative abundance of some gastrointestinal nematode species found in adults in Chizé but not in adults in Trois Fontaines suggest a higher ability to resist to certain parasite species (e.g. some Trichostrongylus variants) in adult roe deer in Trois Fontaines. There is probably a synergistic interaction among the body condition of roe deer, their immune response and the richness of parasites to which they have been exposed [11]. Individuals in lower body condition are more susceptible to parasite infection or/ and high levels of parasitism are suspected to negatively impact body condition [19]. Such negative correlation between physiological performance and parasitic burden has previously been observed to be more pronounced in Chizé than in Trois Fontaines [76]. In addition, in Trois Fontaines, the higher body condition of roe deer [10] and strongyle species diversity may have allowed roe deer to invest in innate immunity [77], which could prevent the establishment or facilitate the expulsion of strongyles [78]. In contrast, roe deer in Chizé, exposed to a lower species richness of strongyles and being in lower body condition, invest more in adaptive response, which is considered less costly than innate response [79].
Contrary to our expectation (H3), the alpha diversity of the parasite communities did not differ between males and females in this study. Indeed, the probability of exposure to gastrointestinal species is equal between males and females as they do not spatially segregate [80], and few sexual differences in their feeding behavior or activity have been observed between males and females [81,82]. In addition, no sex-specific physiological status occurred during the study period, such as nursing young [4] or rutting [83], that could have diverted or limited the allocation of resources to parasite defenses. Sampling during such sex-specific energetically costly periods could be pertinent to consider the mechanisms underlying sex and age differences in gastrointestinal nematodes communities.
Despite the differences in habitat quality, the nemabiome of roe deer was not different among the three sectors of Chizé (H4). This suggests that local ecological variation and the resulting differences in body conditions among roe deer inhabiting these different habitats [40] do not affect the nemabiome. This may also be explained by the movement of roe deer among sectors that can spread the parasites among the feeding areas [40,84].

Detection of parasite species generally associated with domestic ungulates
The occurrence in roe deer of Trois Fontaines of parasites commonly detected in small domestic ruminants, including highly pathogenic species (e.g. H. contortus, B. trigonocephalum), confirmed previous studies showing that wildlife can share some parasites found in livestock (e.g. [33,38]). To our knowledge, only four to six domestic sheep have been bred for several years on the other side of a road surrounding the fenced study site of Trois Fontaines (C. Warnant, personal communication). Even if the fence of the study site is regularly checked for damages, the fence does not totally prevent a few roe deer from leaving and/or entering the reserve (e.g. after the Lothar hurricane in 1999). Such movements may have facilitated the introduction of generalist parasites in the population of Trois Fontaines. The consequences of these pathogenic species on roe deer fitness remain however to be determined.

Conclusion
In the present study we investigated the community of gastrointestinal nematode species in two isolated populations of roe deer using ITS2 rDNA nemabiome metabarcoding. Most gastrointestinal nematodes detected in the two sites were parasites classically found in roe deer. Even though the two populations inhabit fenced areas, we also observed gastrointestinal nematodes usually found in domestic ungulates in one of the populations (Trois Fontaines). Despite a fence restraining movements of animals from and to the study site, migrations of a few individuals may have contributed to the movements of parasites among populations, including generalist parasites from domestic ungulates. This could have important implications for the management of wild populations and the introduction of non-native species to wildlife communities. Therefore, it appears highly relevant to study roe deer and domestic ungulates living in sympatry in anthropogenic landscapes dominated by farmland to better understand the dynamic of gastrointestinal nematodes among wild and domestic populations and the risk of cross-transmission. This work also illustrated the power and limitations of ITS2 nemabiome metabarcoding for wildlife studies. It allows the screening of large numbers of samples for a large number of parasites, with identification to the species level in many cases. However, some parasites could still only be identified to the genus level emphasizing the need to continue to enrich the DNA sequence databases for more comprehensive identification.