Introduction

The broad definition of habitat fragmentation covers a decrease in the area and a change in the spatial configuration of a given habitat’s patches, such as a higher number of patches, an increase in their isolation or longer edges (Fahrig 2003). This process is considered to be one of the major causes of biodiversity loss (Palumbi 2001; Haddad et al. 2015). Indeed, accumulating studies rooted in island biogeography theory have shown a decrease of species richness as a result of habitat fragmentation (Collinge 2009; Losos and Ricklefs 2009; Lindenmayer and Fischer 2013). Considering biodiversity at the genetic level, the reduction in habitat patches area and the increase in their isolation are expected to strengthen genetic drift and to lower gene flow between them, and ultimately to lead to a decrease in genetic diversity within populations and an increase in their differentiation (Keyghobadi 2007). A large number of empirical studies validated these predictions (e.g., Keyghobadi 2007; Aguilar et al. 2008; Rivera-Ortíz et al. 2015). Nonetheless, the response to fragmentation appears to be modulated according to landscape characteristics, such as age and level of fragmentation, and ecological traits of the species, such as habitat specialization or dispersal capacities (Driscoll and Weir 2005; Schleicher et al. 2011; Rivera-Ortíz et al. 2015; Khimoun et al. 2016a). Moreover, the expected effects of habitat fragmentation on genetic diversity are theoretically detectable at a spatial scale positively related to dispersal potential of species (Slatkin 1987). However, it is noteworthy that even on mobile organisms such as birds, a few works reported small-scale spatial genetic structure (Arnoux et al. 2014; Bertrand et al. 2014), as a result of habitat fragmentation (Brown et al. 2004; Lindsay et al. 2008; Porlier et al. 2012; Woltmann et al. 2012; Harrisson et al. 2013; Khimoun et al. 2016a).

Genetics is a useful tool to inform conservation management, shedding light on how genetic variation affects individual fitness and extinction risk of populations (Spielman et al. 2004; Sunnucks 2011). It is now widely recognized that preserving a species’ evolutionary potential, through functional genetic variation, is critical to ensuring the long-term persistence of populations, especially in the face of environmental changes. However, studies addressing the genetic consequences of habitat fragmentation are often based on neutral molecular markers. Although neutral genetic variation can be informative for demographic inferences, it is, however, a poor predictor of the adaptive genetic diversity that controls population viability or evolutionary potential (McKay and Latta 2002; Holderegger et al. 2006).

In this context, genes of the immune system have received and still deserve particular attention because infectious diseases are a major driver of the ecology and evolution of natural populations (Ostfeld et al. 2010), and because there is a direct link between the genetic diversity of hosts and their ability to cope with pathogens (Altizer et al. 2003; Spurgin and Richardson 2010). In addition, it is now clear that emerging and reemerging infectious diseases are exposing wildlife and human health to threats at a growing rate (Daszak 2000; Jones et al. 2008). Accordingly, there is a crucial need to understand how habitat fragmentation regulates immunogenetic diversity of hosts beyond academic interest to provide key information for conservation and management decisions. Not only may habitat fragmentation weaken hosts’ resistance through genetic erosion, but it may also expose natural populations to new pathogens in fragmented landscapes in very diverse ways (Patz et al. 2000; Collinge 2009; Gottdenker et al. 2014).

Understanding the relative contribution of drift and selection in shaping patterns of functional diversity in fragmented populations is a central purpose in conservation genetics (Allendorf et al. 2013; Fraser et al. 2014; Frankham et al. 2017). A method commonly used to test for the role and type of selection in shaping adaptive genetic variation consists of contrasting spatial patterns of variation between candidate loci and neutral markers (Bernatchez and Landry 2003; Biedrzycka and Radwan 2008; Strand et al. 2012; Radwan et al. 2014; Santonastaso et al. 2017; Belasen et al. 2019). Indeed, if selection outweighs demographic processes (drift and migration), level of differentiation in candidate loci between populations is expected to be significantly decreased under stabilizing selection compared with neutral markers, whereas it should be higher if populations occurring in heterogeneous environments are exposed to directional selection and are potentially locally adapted. Similar patterns in both candidate and neutral loci reflect the absence of selection or selection pressures that are weaker than demographic processes. In this context, it is now recognized that looking at the structure of adaptive genetic diversity at a small spatial scale is a powerful approach (Richardson et al. 2014). Indeed, micro geographic adaptation processes may result in local adaptations which are not observable when individuals from close populations are clustered to evaluate differences only on larger scales (Gonzalez-Quevedo et al. 2016; Pearson et al. 2018).

Here, we tested whether forest fragmentation affects genetic diversity in a tropical bird species, the wedge-billed woodcreeper Glyphorynchus spirurus, at a small spatial scale. We contrasted patterns of neutral (microsatellite markers) and immunogenetic diversity. Even if studies of immunogenetic diversity in natural populations have been dominated by the use of the major histocompatibility complex (MHC) genes, the complexity of the immune system involves many other genes, which account for a large proportion of the variability in host resistance (Acevedo-Whitehouse and Cunningham 2006; Vinkler and Albrecht 2009; Bollmer et al. 2011). In addition, the frequent duplication of MHC genes in many species, especially passerines, still presents technical challenges to get information at individual loci (Bollmer et al. 2010). Therefore, we estimated immunogenetic diversity from toll-like receptors (TLR) genes, which are part of the innate immune response. TLRs can be amplified in a wild range of non-model species and therefore constitute an interesting marker for immunogenetics studies conducted in natural populations (Alcaide and Edwards 2011; Grueber et al. 2013). They detect invariant features of invading microbes, called “pathogen-associated molecular patterns” (PAMPs) which, through an intracellular signaling cascade, activate an inflammatory response (Alcaide and Edwards 2011). These genes should also evolve in direct response to pathogen-mediated selective pressures (Quéméré et al. 2018), as already suggested by a few studies conducted on several species (Yilmaz et al. 2005; Wlasiuk and Nachman 2010; Alcaide and Edwards 2011; Areal et al. 2011; Grueber et al. 2014).

In order to assess potential parasites pressures and their variation associated with habitat fragmentation, we screened for the presence of ticks in all birds. Indeed, ticks can impact fitness components of hosts, through a reduction of adults’ survival or a reduction of reproductive success due to an increase of nestling mortality (Lehmann 1993). Ticks are also the most important vectors of animal infectious diseases (Jongejan and Uilenberg 2004; Colwell et al. 2011), and tick prevalence or load may be good proxies of the pathogen diversity or load they can transmit. In addition, several putative bacterial pathogens have been recently described in ticks from French Guiana (Binetruy et al. 2020a, b). Finally, several studies suggest interactions between immune gene diversity and tick infestation (Becker et al. 2020). Tick larva burden has been shown to depend on individual TLR4 genotype in a population of water voles (Gavan et al. 2015). Tschirren (2015) reported that allele frequencies of TLR2 in the bank vole were positively associated with the risk of infection by Borrelia burgdorferi, a common tick-transmitted pathogen responsible for the Lyme disease. In addition, bioactive molecules in tick saliva can dampen the TLR4 signaling pathway and increase expression of TLR2, thus allowing ticks to subvert their host immune response (Oliveira et al. 2010).

The specific objectives of this study were: (1) to investigate whether forest fragmentation affects population genetic diversity at small spatial scale; (2) to contrast the patterns of neutral and adaptive genetic diversity to assess the relative effect of demography (genetic drift and migration) and selection pressures; and (3) to assess whether patterns of immunogenetic variability could be associated with variation of parasite pressures.

Materials and methods

Study area and study species

French Guiana is located on the northern Atlantic coast of South America. With a population of 260,000 inhabitants (in 2015) in a territory of 84,000 km2, this is the least densely populated country in tropical America. At the same time, this is the country with the highest forest cover (98% of land area, FAO 2010). The FAO (2001) ranked French Guiana world’s number one in 2000 according to forest area per capita (with 45.6 ha/capita in this country compared with a mean of 2.9 ha/capita in Tropical South America). However, geographical distributions of people and forest resources are highly uneven in this country. The inland forest has always been sheltered from degrading human activities so far, whereas forest habitat in the coastal area has been tremendously modified by the human settlement, especially following arrival of Europeans in the sixteenth century (DIREN Guyane 2007; Calmont 2012). Indeed, 92% of the human population is currently concentrated on the coastal area representing <5% of the country’s area (Galochet and Morel 2015). In addition, human population size has increased dramatically in the last years: there were 157,000 inhabitants in 1999, 245,000 inhabitants in 2013, and 513,000 inhabitants are expected in 2050 (INSEE 2019). As a result, deforestation has occurred in this area mainly because of urbanization, infrastructure development, and land conversion for agriculture, and this ongoing process is expected to increase in the future. For instance, about 7% of forest habitat area in the coastal region was lost in favor of farmland (resulting in an increase of 46% in area) and artificial surfaces (resulting in an increase of 22% in area) in the last 10 years (ONF 2017). Finally, the development of gold mining activities in the last century and a few projects of road links are potentially threatening the integrity of the inland forest in the near future (DIREN Guyane 2007).

The wedge-billed woodcreeper (Glyphorynchus spirurus) is a small passerine (13–16 cm; 10.5–21 g) endemic to the Americas, occurring from southern Mexico to eastern Brazil. This species is found in the understorey of old-growth forests, where it feeds mainly on arthropods and sometimes on vegetable matter, but it also occurs in forest edges and secondary forests (Marantz et al. 2019). This is a territorial species (Blake and Loiselle 2012), which disperses only over short distances (rarely beyond 8 km) when moving within continuous forests (Van Houtan et al. 2007). The wedge-billed woodcreeper is a particularly relevant species for studying consequences of habitat fragmentation, as it is relatively common in both forest fragments and continuous forest.

Field sampling

Field sampling took place during the dry season (between July and November, from 2012 to 2014) in French Guiana. Twelve forest sites were sampled in total: 5 sampling sites were located in the inland continuous forest (subsequently named “control sites”) and 7 sampling sites were located in forest fragments in the coastal area that were isolated from each other and from the inland forest (subsequently named “fragmented sites”; Fig. 1). The area of patches containing fragmented sites was within the same order of magnitude (5–50 km2) compared to the area of the forest containing control sites (>20,000 km2). Then, we did not quantify patch area because it would have resulted in the same dichotomy. Geographic distance between sampling sites ranged between 3.5 and 212 km (mean = 90 km), and both short and large distances were observed between control as well as fragmented sites. We captured birds with mist nets in the morning. Each bird was ringed and about 20 μl of blood was collected from the brachial vein using sterile needles and heparinized tubes. Blood samples were stored in 500 μl of QLB (Queen’s Lysis Buffer; Seutin et al. 1991) until molecular analyses.

Fig. 1: Location of sampled sites of wedge-billed woodcreeper in French Guiana and population genetic structure detected by STRUCTURE algorithm.
figure 1

The green color on maps represents the forest habitat. Individuals were sampled in fragmented forest fragments (named Fxx in white label) and in the continuous inland forest (i.e., control sites named Cxx in black label). For each population, pie chart represents the membership proportion in the three genetic clusters, based on microsatellite genotypes. The size of each chart is proportional to sample size. Figures are based on the run with the highest posterior probability of the data for each K-number of genetic clusters. Values of increase in the posterior probability of the data (IPPD) are provided for each model in the box and the best model is indicated in red.

A total of 254 wedge-billed woodcreeper were captured, corresponding to an average of 21 birds per sampled site (min = 15, max = 32; Supplementary Material S4).

Lab work

Total DNA extraction was performed on 240 μl of QLB-stored blood, according to a standard phenol–chloroform protocol optimized from Hillis et al. (1996), initiated by a first step of digestion with 10 μl of proteinase K (overnight incubation at 55 °C).

Microsatellite genotyping

Following the approach described in Khimoun et al. (2016b), we optimized a set of 14 microsatellite loci for individuals’ genotyping (Supplementary Material S1). Based on preliminary tests (amplification success, quality of peak profile, and Hardy–Weinberg equilibrium (HWE)), we retained 2 loci previously isolated in the wedge-billed woodcreeper (Milá and Bardeleben 2005) and 12 cross-species amplifying loci developed from expressed sequence tags (EST-SSR; Dawson et al. 2010, 2013). Although EST-SSRs are highly conserved, they have been proved valuable molecular markers to reveal weak genetic structure (Khimoun et al. 2017). A full description of amplification protocols is provided in Supplementary Material S2. Loci were amplified in simplex using a fluorescent labeling method of PCR products with M13-tail (Schuelke 2000). PCR products were multiloaded and analyzed in an automated sequencer (ABI3730), and then allele scoring was performed using GENEIOUS 8 (Kearse et al. 2012).

TLR genotyping

We used primers developed by Alcaide and Edwards (2011) and Grueber et al. (2012) to amplify nine of the ten TLR loci present in birds (TLR7 was not retained because of evidence for duplication events; Cormican et al. 2009; Grueber et al. 2015). These primers target exons encoding the leucine-rich repeat (LRR) regions of the extracellular domain of each TLR, involved in pathogens recognition (Alcaide and Edwards 2011; Grueber et al. 2012). Preliminary tests were performed using different PCR conditions to assess TLRs amplification success on agarose gel and sequence quality (i.e., unambiguous sequence editing), using a subsample of 25 individuals. From these preliminary results, we selected four TLRs for amplification and sequencing on the whole sample set: TLR3 and TLR21 binding viral ligands and TLR4 and TLR5 binding bacterial ligands (Keestra et al. 2013). The remaining five TLRs (TLR1LA, TLR1LB, TLR2A, TLR2B, and TLR15) were discarded, as amplification success was not consistent among individuals. A full description of PCR protocols is provided in Supplementary Material S2.

PCR products were sequenced in an automated sequencer (ABI3730) and sequences were edited and aligned with GENEIOUS 8 (Kearse et al. 2012). Single nucleotide polymorphisms (SNPs) were confirmed from examination of both forward and reverse chromatograms. Individual haplotypes were resolved using the program PHASE (Stephens et al. 2001; Stephens and Donnelly 2003) implemented in DNAsp v5 (Librado and Rozas 2009) with the following parameters, using three independent runs with different seed values: 1000 iterations, a thinning interval equal to 1, and a burn-in value of 100. We used results from the run with the highest goodness of fit and we kept all haplotype pairs with a probability above 0.6. This threshold reduces the number of unresolved genotypes with little or no increase of false positives (Garrick et al. 2010). Unique observations of putative haplotypes were confirmed by reamplification and resequencing.

Parasites screening

To characterize potential drivers of local selection acting on bird populations, we quantified the prevalence of ticks in each sampling site. We screened the head of each bird for ticks, and we collected a few ticks at several sites (from 15 individual birds from different sites) for morphological and genetic identification following Binetruy et al. (2019).

Statistical analyses

We assessed both within- and between-population genetic diversity through a multilocus approach when using microsatellite loci. On the contrary, TLR loci were analyzed individually, because they may have evolved independently in response to different selection pressures (Grueber et al. 2015).

Genetic diversity within populations

Linkage disequilibrium (LD) between all pairs of microsatellite loci and deviation from HWE for each microsatellite locus were tested in each population using GENEPOP 4 (Rousset 2008). Significance of tests was assessed after applying Holm–Bonferroni correction accounting for the number of populations for LD (to test for nonindependence between each locus pair across populations) and HWE (to test for consistent disequilibrium due to potential amplification issues across populations) tests. Expected heterozygosity (HE) and allelic richness (A) are two common measures of within-population genetic diversity. Several studies have reported no or only a weak congruence between HE estimated from microsatellites (SSR-HE) and from genome-wide SNPs, suggesting that SSR-HE may not be a useful proxy of genome-wide diversity. On the contrary, as significant correlations between allelic richness based on microsatellites (SSR-A) and genome-wide SNPs have been described, SSR-A is assumed to be a good surrogate for genome-wide diversity (Väli et al. 2008; Vilas et al. 2015; Fischer et al. 2017). Thus, within-population genetic diversity was estimated from allelic richness calculated for each TLR and averaged across microsatellite loci. TLRs allelic richness was calculated from (1) nucleotide sequences (i.e., considering all SNPs of the sequence) and (2) amino acid sequences (i.e., considering only SNPs causing non-synonymous substitutions). In order to be compared among sampled sites, allelic richness based on microsatellites and on each TLR locus was estimated from a standardized sample size (i.e., standardized allelic richness, As) of ten individuals, using a rarefaction procedure implemented in ADZE 1 (Szpiech et al. 2008). We compared standardized allelic richness between control and fragmented sites using a Wilcoxon–Mann–Whitney test, using R 3.6.3 (R Core Team 2020). We also tested the correlation in allelic richness between microsatellite loci and each TLR locus with a Spearman’s rank-order correlation, in R 3.6.3 (R Core Team 2020). A lack of significant correlation would suggest that TLR and microsatellite genetic variations have been structured by different processes in wedge-billed woodcreeper populations.

Genetic differentiation between populations

Population genetic structure was assessed following two complementary approaches. First, it was evaluated without a priori information about geographic location of sampled individuals. We looked for the number of genetic clusters (K) represented in the data set, under the assumption of HWE and linkage equilibrium between the loci within clusters. This analysis was performed from microsatellite genotypes only, using a Bayesian clustering approach implemented in STRUCTURE 2.3.3 (Pritchard et al. 2000). Several models were fitted to the data with values of K ranging from 1 to 12, and the optimal value of K was determined from the successive increase of the posterior probability of the data (IPPD) at each increment of K (Garnier et al. 2004). For each value of K, ten independent runs of one million iterations were performed after a burn-in period of 100,000 iterations for a model with admixture (Falush et al. 2003), correlated allelic frequencies, and with the LocPrior option (Hubisz et al. 2009). In the second approach, we considered each sampling site as a distinct population. We calculated pairwise FST (Weir and Cockerham 1984) over all populations for microsatellite loci and for each TLR gene (using either nucleotide or amino acid sequences). We also calculated pairwise GST because this metric is independent of the level of genetic variation (Hedrick 2005). We tested genotypic differentiation between pairs of populations using 10,000 permutations and the Holm–Bonferroni correction (accounting for the number of population pairs), using both microsatellite loci and each TLR gene. Estimations and tests of differentiation were performed using GenAlEx 6.5 (Peakall and Smouse 2012). We also tested the correlation between patterns of genetic differentiation assessed from each TLR and microsatellites, using Mantel tests performed using the R package vegan (Oksanen et al. 2019).

Finally, we tested whether fragmentation of the forest cover may lead to a decrease of landscape functional connectivity, by impeding gene flow among bird populations. To do so, we compared the relative support of two alternative models to predict observed patterns of genetic structure. First, the isolation-by-distance IBD model, corresponding to a null model, assumes that gene flow among populations is only constrained by geographic distance and that genetic differentiation between populations is expected to increase with Euclidian distance among populations. The alternative model, a landscape connectivity model, accounts for landscape habitat heterogeneity and its potential consequences in promoting/impeding gene flow, by relying on ecological distances. Here, ecological distances among populations integrated the potential cost associated with bird dispersal through a non-forest habitat. Then bird dispersal and gene flow between populations are assumed to occur following the path associated with the least cost, as described in the least cost LCP model (Adriaensen et al. 2003). In order to compute least-cost distances, we constructed a raster surface of dispersal cost values from a forest cover raster from 2010 with a 30 × 30 m resolution (Hansen et al. 2013). Each pixel received a cost value that was inversely related to the proportion of forest in this pixel (ranging from a cost value of 0 when the whole pixel is forested to a cost value of 1 when there is no forest in the pixel). Then we calculated the cumulative cost along the least-cost path between each population pair, using the R package gdistance (Van Etten 2017). This measure of ecological distance (hereafter named NF-distance) reflects the minimum distance that a bird has to cross in non-forest habitats between each pair of populations. We tested the relative performances of geographic (IBD model) and NF (LCP model) distances to explain population pairwise genetic differentiation (using either FST or GST) assessed from microsatellites and each TLR locus, using Mantel tests performed with the R package vegan (Oksanen et al. 2019). Finally, we compared the regression slopes between microsatellites and TLR loci considering either geographic distances or NF-distances, using the MantelPiece R-script version 1.0 (http://www.erikpostma.net/resources.html) that provides 95% confidence intervals for the difference in slopes using jackknifing.

Tests of selection on TLR loci

We used two different selection tests to characterize the evolutionary patterns of TLR polymorphism. First, we ran the codon-based Z-test of selection, which detects selection at the locus level (Pond and Frost 2005). This test compares the frequency of non-synonymous mutations per non-synonymous site and the frequency of synonymous mutations per synonymous site. This test was performed for each TLR locus in MEGA version X (Kumar et al. 2018) with default settings. Because this test can only detect strong signal of selection (Pond and Frost 2005), we also used the fast unbiased Bayesian approximation (FUBAR) (Murrell et al. 2013) method available through the Datamonkey server (http://datamonkey.org), which detects selection at the codon level. This test allows determining whether specific codons, instead of the entire sequence, evolved under positive or purifying selection. A codon under purifying selection would play an important structural or functional role, while a codon under positive selection may suggest, in the case of immune genes, an adaptive response to pathogen pressures (Bernatchez and Landry 2003).

Parasite prevalence variation among populations

We tested for variations of tick prevalence between control and fragmented sites using a generalized linear mixed model (“glmer” in R package “lme4”) with binomial errors. This model contained sampling site as a random effect and the type of site (control/fragmented) as a fixed effect. Finally, we tested for the correlation between the prevalence of ticks (i.e., the percentage of infected birds per population) and within-population allelic richness estimated from microsatellites loci and each TLR locus using a Spearman’s rank-order correlation test, in R 3.6.3 (R Core Team 2020).

Results

Genetic diversity within populations

All 14 microsatellite loci were polymorphic with a number of alleles ranging from 3 to 25 (mean = 12, SD = 7). Three microsatellite loci were in Hardy–Weinberg disequilibrium, each one in a different population, meaning that we detected no amplification issues in a particular locus. Regarding LD, seven tests were significant after Holm–Bonferroni correction, but there was no pattern regarding a particular locus pair or population, meaning that the 14 microsatellite loci are statistically independent. Standardized allelic richness averaged over loci ranged from 4.53 to 6.34 (Supplementary Material S4).

Length of TLR sequences ranged from 610 bp for TLR21 to 1105 bp for TLR3 (Supplementary Material S3) and no stop codon or frameshift mutations were detected. Levels of polymorphism were rather heterogeneous among the four TLR (Supplementary Materials S3 and S4). When considering nucleotide sequences, TLR3 and TLR21 were much less variable than TLR4 and TLR5, with standardized allelic richness ranging from 3.75 to 9 alleles per population for TLR3 and TLR21 against 7.96 to 15.40 alleles for TLR4 and TLR5. Such a dichotomy was not observed when considering amino acid sequences. TLR loci rather ranged along an ascending gradient of standardized allelic richness, from TLR21 (with the lowest population allelic richness: 1.0–3.9) to TLR5 (with the highest population allelic richness: 6.30–12.0), with TLR3 (2.88–5.90) and TLR4 (3.58–8.96) presenting intermediate levels of population allelic richness (Supplementary Material S4).

The effect of forest fragmentation on within-population genetic variability depended on the type of genetic markers considered (Table 1). We found that fragmented sites had lower neutral standardized allelic richness compared with control sites (U = 30, p = 0.048). We also found a negative effect of fragmentation on TLR5 standardized allelic richness, when considering both nucleotide (U = 30.5, p = 0.042) and amino acid sequences (U = 31, p = 0.030). On the opposite, within-population variability at the TLR21 locus was higher in fragmented sites than in control sites when considering amino acid sequences (U = 5, p = 0.049), whereas no effect was detected when considering nucleotide sequences (Table 1). Finally, no effect of fragmentation was detected on TLR3 and TLR4 variability (Table 1).

Table 1 Effect of forest fragmentation on within-population allelic richness.

In agreement with previous results, levels of within-population variability assessed from microsatellites and TLR5 were positively correlated, when considering both nucleotide and amino acid sequences (ρ = 0.69, p = 0.014 and ρ = 0.62, p = 0.035, respectively). No significant relationship was detected between variability assessed from neutral markers and the other three functional loci, except for standardized allelic richness of TLR4 nucleotide sequences that was positively correlated with microsatellites allelic richness (Supplementary Material S5).

Genetic differentiation between populations

Neutral genetic differentiation

STRUCTURE clustering analyses based on microsatellite genotypes detected three different gene pools among the sampled individuals (Fig. 1). The main genetic structure in two different genetic pools segregated individuals from three fragmented sites (F05–F07) located close to the main city (i.e., Cayenne) from all other sampled individuals. A finer genetic structure was detected, with a third genetic pool encompassing individuals from the fragmented site F04.

Over 66 exact pairwise tests for genotypic differentiation at microsatellite loci, 37 (56%) were significant after Holm–Bonferroni correction. No genetic differentiation was detected among control sites, whereas almost every population pairs (67%), including at least one fragmented population, were genetically differentiated (except for fragmented sites F01 and F02 that were never differentiated from any control sites; Supplementary Material S6). Multilocus FST estimates ranged from 0 (between several population pairs) to 0.103 (between sites F04 and F06). Pairwise FST values tend to increase from control site pairs to fragmented site pairs, with intermediate levels of differentiation for pairs involving one control and one fragmented site (Supplementary Material S6). GST estimates exhibited the same trend, with values reaching 0.213 between fragmented sites F04 and F06.

Immunogenetic differentiation

Based on TLR genotypes, the number of pairwise tests for genotypic differentiation that remained significant after Holm–Bonferroni correction ranged from 2 for TLR3 to 14 for TLR5 when considering nucleotide sequences and from 0 for TLR21 to 10 for TLR4 when considering amino acid sequences (Supplementary Material S6). All these significant tests involved a fragmented site (except one test that was significant between control sites C03 and C05 for TLR21 nucleotide sequences). Pairwise genetic differentiation, whether estimated from nucleotide or amino acid sequences, and whether using FST or GST, was highly variable. FST estimates ranged from about 0.07 for TLR4 and TLR5 to 0.112 for TLR3 and 0.245 for TLR21 when considering nucleotide sequences, whereas GST estimates were much higher for TLR4 and TLR5 (0.778 and 0.983, respectively) than for the two other loci (0.282 for TLR3 and 0.518 for TLR21; Supplementary Material S6). This illustrates that some population pairs shared hardly any TLR5 nucleotide sequences. We observed the same trend when considering amino acid sequences, with maximum pairwise GST estimates ranging from 0.124 for TLR21 to 0.714 for TLR4.

Correlations between genetic differentiation at microsatellite loci and each TLR locus

Although absolute levels of differentiation were often highly contrasted between FST and GST estimates for TLR loci, correlations between neutral and functional differentiation levels were quite similar when tested using FST or GST estimates (Supplementary Material S7). When considering TLR nucleotide sequences, there was a significant positive correlation between FST values estimated from microsatellites and TLR4 and TLR5 (ρ = 0.61, p < 0.001 and ρ = 0.69, p < 0.001, respectively), whereas this correlation was lower and marginally significant for TLR3 (ρ = 0.3, p = 0.06) and TLR21 (ρ = 0.22, p = 0.06). When considering TLR amino acid sequences, a positive correlation with neutral differentiation levels was observed for TLR4 (ρ = 0.53, p = 0.002) and TLR5 (ρ = 0.37, p = 0.037), but not for TLR3 and TLR21 (Supplementary Material S7).

Test of isolation-by-distance and least-cost-path models

We did not detect any IBD pattern, neither for neutral nor for adaptive loci (Table 2). However, neutral differentiation levels were positively correlated with NF-distance between population pairs (FST: ρ = 0.66, p = 0.003; GST: ρ = 0.63, p = 0.005). Using both nucleotide and amino acid sequences, genetic differentiation assessed from TLR4 and TL5 was positively correlated with NF-distance with correlation coefficients ranging from 0.28 to 0.42 for FST values and from 0.19 to 0.37 for GST values, though p values ranged between 0.06 and 0.11 for FST values and between 0.03 and 0.19 for GST values. In agreement, we detected no difference in the regression slope of pairwise genetic differentiation (estimated from microsatellites, TLR4, or TLR5) and NF-distances. Alternatively, genetic differentiation levels estimated from TLR3 and TLR21 were not explained by NF-distance (p > 0.25 and p > 0.34 when considering FST and GST values, respectively). When considering FST values, the slope of the regression between differentiation based on TLR3 and NF-distance was lower than the slope of the regression considering neutral differentiation. We observed a similar trend for amino acid sequences of TLR21. When considering GST values, the slope of LCP models based on amino acid sequences of TLR3 and TLR21 was lower than the slope of LCP model based on microsatellite markers (Table 2).

Table 2 Tests of IBD model (correlation between pairwise genetic and geographic distances) and LCP model (correlation between pairwise genetic and NF-distances) for microsatellites and TLR loci (Spearman correlation coefficient and probability).

Tests of selection on TLR loci

Overall, we found a lower frequency of non-synonymous substitution per non-synonymous site compared to the frequency of synonymous substitutions per synonymous site at every TLR locus (Z-test = −2.48, p = 0.01 for TLR3; Z-test = −2.46, p = 0.02 for TLR4; Z-test = −3.05, p = 0.003 for TLR5; Z-test = −1.7, p = 0.10 for TLR21). In agreement, the FUBAR method detected many codons under purifying selection across the four TLRs. More precisely, we found 9, 8, 19, and 17 codons under purifying selection for TLR3, TLR4, TLR5, and TLR21, respectively. However, the FUBAR method also reported evidence of positive selection (i.e., an excess of non-synonymous substitutions) on ten codons, across the four TLRs, all located within the extracellular region of the gene involved in ligand binding (LRR region). We observed these codons in TLR3 (1 codon), TLR4 (3 codons), and TLR5 (6 codons).

Parasite prevalence variation among populations

The global prevalence of ticks in wedge-billed woodcreeper in French Guiana was 18.4%. The sample of ticks taken from wedge-billed woodcreeper belonged almost exclusively to one species, Amblyomma longirostre (Acari: Ixodidae; 33 ticks out of the 35 collected). We collected two other species represented by a single individual, Amblyomma varium and Amblyomma geayi. There are large differences in prevalence between sites. Birds were completely free of ticks in some sites (F01, C01, C03, or C04), but local prevalence reached 67% in site F05 (Supplementary Material S4). Tick prevalence was significantly higher in fragmented sites than in control sites (likelihood ratio test: χ2 = 6.05, df = 1, p = 0.014; Fig. 2).

Fig. 2: Correlation between tick prevalence and within-population allelic richness at microsatellite loci.
figure 2

Black and white dots represent control and fragmented sites, respectively.

Finally, we found a negative correlation between tick prevalence and within-population genetic variability assessed from neutral markers (ρ = −0.70, p = 0.017; Fig. 2) and TLR5 (ρ = −0.60, p = 0.037 and ρ = -0.9, p < 0.001 based on nucleotide and amino acid sequences, respectively). This correlation was nonsignificant when considering the three other TLR (Supplementary Material S8).

Discussion

The present study highlights the negative consequences of forest fragmentation in an understorey tropical bird species. Our results revealed an erosion of neutral genetic diversity and a substantial genetic differentiation among fragmented populations resulting from a decrease in landscape connectivity. Such a decrease in population connectivity leads to the divergence of three distinct genetic pools, with two of them encompassing only fragmented sites. Patterns of neutral and functional genetic diversity and structure were not always congruent, depending on TLR loci. This suggests a role of selective pressures potentially related to variations in parasite prevalence such as those we observed between control and fragmented sites for ticks.

Regarding microsatellite genetic diversity, we detected a lower allelic richness within sites sampled in forest fragments, compared to sites sampled within the continuous forest (i.e., control sites). As deforestation in French Guiana occurs in the coastal region, forest fragments are thus located at the margin of forest dweller distribution, whereas control sites are located more in the inland part of their distribution. Such an erosion of neutral genetic diversity observed in forest fragments might result from this spatial distribution of fragmented/control sampled sites. According to the central-margin hypothesis (Eckert et al. 2008), the effective population size and gene flow among populations are expected to decline toward geographical range margins (Lesica and Allendorf 1995), leading to (1) a decreased within-population genetic diversity and (2) a higher genetic differentiation among populations, from the center to the periphery of species geographic ranges. However, according to the wide geographic range of the wedge-billed woodcreeper, extending from southern Mexico to eastern Brazil (covering around 4500 km from north to south and 5000 km from west to east), control sites sampled at 20 km from the coast must also be considered at the margin of the species geographic range. Thus, the central-margin hypothesis, alone, cannot explain the erosion of neutral genetic diversity observed in forest fragments, putting forward the hypothesis of a decreased landscape functional connectivity driven by forest fragmentation. This hypothesis is supported by patterns of neutral genetic structure observed among sampled sites.

Neutral genetic differentiation among sampled sites was not explained by an isolation-by-distance model, highlighting that gene flow among wedge-billed woodcreeper populations is not constrained by geographic distance. Indeed, there was no differentiation between sites C1 and C5 separated by 200 km of continuous forest, whereas sites F05 and F06 located in different forest patches were significantly differentiated although separated by only 3.5 km. In contrast, genetic distances were significantly explained by ecological distances that integrate the nature of habitat (forest vs. non-forest) between each pair of populations into the calculation of least-cost distances. Indeed, we observed a decrease of population genetic connectivity along with the increase of NF-distance between populations. This result reveals that forest fragmentation induces a decrease in landscape functional connectivity by impeding gene flow and increasing genetic differentiation of fragmented populations. Individuals from the most fragmented sites (F04–F07) constitute particular gene pools, differentiated from the other genetic cluster encompassing the remaining sampled individuals. We thus reported a striking case of substantial genetic differentiation at a very restricted spatial scale, in the absence of physical barriers and for bird species with no particular reproductive system or social structure. All these results underline that the negative effects of habitat fragmentation on intraspecific genetic diversity (i.e., depletion of within-population genetic diversity and divergence of genetic pools) may be substantial also for highly mobile species such as birds and for widely distributed and continental species. Based on a binary vision of landscape heterogeneity (forest vs. non-forest habitat), our results pinpoint the need for further investigation of the relative dispersal cost of each non-forest habitat. However, such landscape genetic analyses would require a higher number of sampled sites. This would provide insight about how matrix composition and configuration may mitigate the consequences of forest fragmentation on population connectivity and which areas deserve particular interest in terms of conservation/restoration.

When considering between-population genetic differentiation, TLR4 and TLR5 paralleled the neutral pattern observed from microsatellites, whereas no such relationship was found for the two remaining TLRs (TLR3 and TLR21). These results show that the dominant evolutionary force shaping immunogenetic diversity may be different depending on loci. Indeed, the congruence between neutral, TLR4, and TLR5 patterns of genetic structure suggests that genetic diversity and structure observed for TLR4 and TLR5 loci are mainly shaped by neutral demographic processes (i.e., genetic drift and gene flow) that we found to be substantially altered by forest fragmentation. In agreement, we found an erosion of TLR5 within-population genetic diversity in fragmented sites. No such pattern was detected for TLR4 and the time lag between gene flow reduction and depletion of genetic diversity may explain this result (Frankham et al. 2010). Such a main effect of genetic drift in shaping TLR variation was previously reported in oceanic island populations of Berthelot’s pipit, as patterns of TLR diversity were explained by the severity of the bottlenecks that occurred during island colonization (Gonzalez-Quevedo et al. 2015). Here, we show that drift may also have pronounced effects on contemporary evolution of functional diversity as a response to habitat fragmentation, which might impair the wedge-billed woodcreeper ability to respond rapidly to emerging infectious disease. Conversely, when considering TLR3 and TL21, the observed patterns of genetic structure were discordant from a neutral pattern. This result suggests that neutral evolutionary processes are not the main drivers of TLR3 and TLR21 current patterns of genetic diversity, suggesting a role of selection counteracting the effect of demographic processes (Yilmaz et al. 2005; Wlasiuk and Nachman 2010; Alcaide and Edwards 2011; Areal et al. 2011; Grueber et al. 2014). Indeed, the lower genetic differentiation observed based on TLR3 and TLR21 compared to levels of neutral genetic differentiation suggest stabilizing selection. Such a contrast in the dominant evolutionary processes observed across TLR loci, and more generally observed for genes across the immune system (Evans et al. 2010; Downing et al. 2010; Bateson et al. 2015; Knafler et al. 2017; Minias et al. 2018; Sagonas et al. 2019), is commonly reported and could be explained by the different classes of recognized pathogens. Indeed, TLR3 and TLR21 are located in the endolysosomal compartment and are proposed to detect DNA or RNA molecules of viruses, whereas TLR4 and TLR5 target bacteria as they are located on the surface of eukaryotic cells and detect a lipid molecule located on the surface of bacteria (Lipopolysaccharide) and a preserved protein (flagellin) conforming the bacteria flagellum, respectively (Areal et al. 2011; Keestra et al. 2013). It is worth mentioning that our results provide an additional instance of lower diversity of viral-binding TLRs (here, TLR3 and TLR21) compared to bacterial-binding TLRs (here, TLR4 and TLR5), that might be related to the far lower diversity of viral PAMPs compared to bacterial PAMPs (Cavaillon 2017), as viral nucleic acids are highly conserved among many viruses (Kawai and Akira 2006).

Several studies reveal that immune gene diversity may be driven by parasite-mediated selection, though determining the precise mechanism is still challenging (Spurgin and Richardson 2010). Our results regarding the molecular evolution of TLR in the wedge-billed woodcreeper were in agreement with previous results reported in two other bird species, the Lesser Kestrel and the House Finch (Alcaide and Edwards 2011). Indeed, our results showed that these markers are mainly evolving under purifying selection, and thus corroborate the classical view claiming that polymorphism of innate immune genes has been strongly optimized by natural selection and evolve under purifying selection (Parham 2003). However, several loci showed evidence of positive selection in the ligand-binding region, thus suggesting that innate immune genes involved in pathogen recognition could also evolve in direct response to pathogen-mediated selective pressures.

Interestingly, we detected a higher tick prevalence in fragmented sites that might be explained by differences in abundance of main hosts of adults’ ticks. Indeed, while immature stages of A. longirostre almost exclusively feed on passerine birds, adults feed instead on arboreal mammals such as sloths but also rodents of the Erethizontidae family (porcupines) (i.e., Guglielmone et al. 2014). In addition, rodents are known to proliferate in fragmented landscapes (Keesing et al. 2009; Young et al. 2015; Galetti et al. 2015), potentially resulting from the absence of predators or competitors, and from an increase in resource availability due to the positive effect of forest edges (Dirzo et al. 2014; Young et al. 2015).

Finally, the higher tick prevalence in fragmented sites could contribute to maintaining or even generating a higher immunogenetic diversity in these sites, as we observed for TLR21, through balancing selection (Bernatchez and Landry 2003; Eizaguirre et al. 2012). Indeed, ticks are the most important vectors of animal infectious diseases: they transmit more pathogens than any other arthropod (Jongejan and Uilenberg 2004; Colwell et al. 2011), and tick prevalence or load may be a good proxy of the pathogen diversity or load they can transmit. Recent surveys in French Guiana revealed the presence of several putative bacterial pathogens in bird-associated ticks, and notably in A. longirostre, the most common tick we found on Glyphorynchus spirurus (Binetruy et al. 2020a, b). This includes an intracellular bacteria, Rickettsia amblyommatis (Binetruy et al. 2020a), and a novel Borrelia species, Candidatus Borrelia mahuryensis, that is closely related to Lyme disease agents and infects >5% of A. longirostre specimens (Binetruy et al. 2020b). A previous study also reported a positive association between tick load and MHC diversity (Radwan et al. 2014). Thus, regarding TLR3 and TLR21, (1) the lack of genetic structuration related to forest fragmentation and (2) the lowest levels of population genetic differentiation compared to neutral pattern might be explained by a pool of a few viruses shared between controlled and fragmented sites, leading to the maintenance of similar allelic frequency by stabilizing selection. Indeed, viruses can propagate between different ecosystems and lead to homogeneous communities over large scale (Sano et al. 2004; Breitbart and Rohwer 2005). These results call for characterizing the local microbiota, both in ticks and birds to test the effect of forest fragmentation on these communities.

In conclusion, our results clearly showed that habitat fragmentation negatively impacts population genetic connectivity, even at a local scale for a priori mobile species such as birds. They also show that the NF-distance, as a good proxy of gene flow, may provide practitioners with an easy-to-measure metric to assess landscape functional connectivity. Nevertheless, consequences of habitat fragmentation on genetic diversity can be heterogeneous depending on genes that are considered, and pathogen selective pressures may be higher in fragmented environments, thus potentially leading to the maintenance of an adaptive genetic diversity despite an erosion of neutral genetic diversity. Finally, the contrasted results obtained across the four TLR loci pinpoints the relevance of these immune genes in a candidate gene approach in immunogenetics. As they suffer less technical drawbacks, they may also offer an interesting alternative to the widely used MHC genes in conservation genetics.