<figure><img src="image_1.png" /></figure>

HAL
open science

# Genomic analysis of geographical structure and diversity in the capercaillie (Tetrao urogallus)

Lídia Escoda, Josep Piqué, Ladislav Paule, Kévin Foulché, Emmanuel Menoni, Jose Castresana

<figure><img src="image_2.png" /></figure>

## To cite this version:

Lídia Escoda, Josep Piqué, Ladislav Paule, Kévin Foulché, Emmanuel Menoni, et al.. Genomic analysis of geographical structure and diversity in the capercaillie (Tetrao urogallus). Conservation Genetics, 2023, 25, pp.277 - 290. 10.1007/s10592-023-01567-6 . hal-05352032

HAL Id: hal-05352032

https://ofb.hal.science/hal-05352032v1

Submitted on 6 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Genomic analysis of geographical structure and diversity in the capercaillie (Tetrao urogallus)

Lídia Escoda<fnref n="1" /> · Josep Piqué<fnref n="2" /> · Ladislav Paule<fnref n="3" /> · Kévin Foulché<fnref n="4" /> · Emmanuel Menoni<fnref n="4" /> · Jose Castresana<fnref n="1" />

Received: 14 November 2022 / Accepted: 29 August 2023

© The Author(s) 2023

## Abstract

The capercaillie is widespread throughout the boreal forests of northern Eurasia, but faces serious conservation challenges in the southernmost mountain ranges, where the populations are fragmented and some are critically endangered. To develop effective conservation strategies for these populations, it is essential to have information on both their genetic diversity and the genetic structure of the species. In this work, we used a reduced representation (ddRAD) genomic sequencing technique to analyze the genetic structure of the capercaillie across its European range and to assess the inbreeding levels in some of the most threatened populations. Our population structure analysis suggested the existence of two evolutionarily significant units, one formed by the two populations from the Iberian Peninsula and the other by the populations from the rest of Europe. Genetic diversity showed a significant decrease in the Iberian populations with respect to the other European populations. An isolation-with-migration model supported these results and allowed us to estimate the parameters of the population tree. The inbreeding coefficients estimated for the Iberian capercaillies showed relatively low levels in the Pyrenees. However, some individuals with very high inbreeding values were detected in the Cantabrian Mountains, suggesting that some of its subpopulations are substantially isolated. The population structure results and the genomic monitoring method we used to assess inbreeding levels may be crucial for the conservation and recovery of the most endangered capercaillie populations.

**Keywords** Capercaillie · *Tetrao urogallus* · ddRAD · Genetic structure · Inbreeding · Isolation-with-migration

## Introduction

The capercaillie (*Tetrao urogallus*) is a Eurasian species that inhabits the extensive boreal forests of the north as well as the large mountain ranges of southern Europe, including the Balkans, the Carpathians, the Alps, the Jura Mountains, the Pyrenees, and the Cantabrian Mountains (Leclercq and Ménoni 2018). Most capercaillie populations live in mature coniferous forests that are generally open and structurally diverse, although in some parts of the Pyrenees they inhabit

mixed coniferous and deciduous forests and, in some areas, such as the Cantabrian Mountains, the species can be found in deciduous beech, oak, and birch forests (Leclercq and Ménoni 2018; de Juana and Kirwan 2019).

Thirteen subspecies have been described, but it has been suggested that a revision of this taxonomy is necessary (Leclercq and Ménoni 2018; de Juana and Kirwan 2019). Genetic studies on the capercaillie at the global species level have so far been conducted either using mitochondrial genes (Duriez et al. 2007; Rodríguez-Muñoz et al. 2007; Segelbacher and Piertney 2007; Bajc et al. 2011), which provide a partial view of the genetic characteristics of populations, or with microsatellites (Rodríguez-Muñoz et al. 2007), which have a reduced resolution compared to genome-wide data. Most mitochondrial studies detected two major lineages in the capercaillie unrelated to the described subspecies: a southern lineage comprising the Iberian populations (Pyrenees and Cantabrian Mountains) and a boreal lineage encompassing the rest of the Eurasian populations. Moreover, due to the existence of Pyrenean-like mitochondrial haplotypes in the Balkan populations, specifically in

<figure><img src="image_4.png" /></figure>

<footnote>Jose Castresana
jose.castresana@csic.es</footnote>

<footnote n="1">¹ Institute of Evolutionary Biology (CSIC-Universitat Pompeu Fabra), Passeig Marítim de La Barceloneta 37, Barcelona, Spain</footnote>

<footnote n="2">² TRAGSATEC, Madrid, Spain</footnote>

<footnote n="3">³ Faculty of Forestry, Technical University in Zvolen, Zvolen, Slovakia</footnote>

<footnote n="4">⁴ DRAS, Office Français de La Biodiversité, Villeneuve de Rivière, France</footnote>
the Bulgarian populations, it was proposed that there was a recent connection between capercaillie populations located in the southern part of the range (Duriez et al. 2007; Bajc et al. 2011; Klinga et al. 2015). In terms of microsatellites, while some studies supported the singularity of some sub-species (Storch et al. 2006; Rodríguez-Muñoz et al. 2007), others did not find a good correlation between current capercaillie taxonomy and genetics (Klinga et al. 2015).

The capercaillie is considered a species of Least Concern (LC) by the IUCN at the global level (BirdLife 2016). However, many populations are locally endangered. For this reason, the species is legally protected in several western and central European countries, which have mostly small and declining populations (Storch 2007). The decrease in these capercaillie populations has been attributed to multiple factors, some of the main threats being destruction of its forest habitat, overhunting, human disturbance, collisions with cables and fences, ungulate pressure on its habitat, predation pressure, and climate change (Quevedo et al. 2006; Moreno-Opo et al. 2015; Rodríguez-Muñoz et al. 2015; BirdLife 2016; Jahren et al. 2016; González et al. 2017; de Juana and Kirwan 2019; Kämmerle et al. 2020). Of the declining populations, those located at the southwestern edge of the range, in the Iberian Peninsula, have suffered a drastic regression in recent years. In particular, the latest estimates in the Cantabrian Mountains point to less than 200 specimens remaining in this population, which would represent a 90% reduction in population size over the last four decades (Jimenez et al. 2022). The Pyrenean population has also suffered a significant decline since the middle of the last century, particularly on the southern slopes. Thus, a population reduction of 33% was estimated for the southeastern Pyrenees in a single decade (Canut et al. 2021), with a 58% decline in the south-central area over two decades (Gil et al. 2020). Although not as strong, the northern slopes of the Pyrenees have also suffered a slow but significant population erosion over the past two decades (Bal et al. 2021).

Conservation genomics can be useful in several ways when analyzing and monitoring populations in strong regression, and thus can greatly aid in conservation and recovery plans (Steiner et al. 2013; Leroy et al. 2018; Supple and Shapiro 2018; Barbosa et al. 2020; Hohenlohe et al. 2021). First, it is necessary to assess the global genetic structure of the species and determine its evolutionarily significant units (Ryder 1986; Moritz 1994; Funk et al. 2012). This is important because population management may require population reinforcement with specimens from other areas. The use of source populations with the greatest genetic similarity to the affected populations should help ensure that possible evolutionary adaptations are not lost, thereby increasing the probability of success of these programs (Supple and Shapiro 2018; Hohenlohe

et al. 2021). However, it is also important to characterize the genetic health of the populations targeted by conservation plans. This is necessary to identify not only the populations with the worst genetic parameters and greatest conservation needs, but also those that can be used as donors in captive breeding programs or for genetic rescue (Supple and Shapiro 2018). In addition, each individual should be characterized so that the most suitable individuals from potential source populations can be selected. The first level of genetic health information is individual heterozygosity, which reflects the genetic diversity of each individual. This value allows easy comparison of individuals, populations, and species and is highly important as a measure of population viability (Scott et al. 2020). However, it is an average of the overall genetic diversity of the genome and does not reflect the existence of genome regions with lower diversity. Another parameter that provides more complete information on the genetic health of an individual is the inbreeding coefficient (Kardos et al. 2018), which measures the proportion of an individual's genome in homozygosis due to mating between close relatives. High inbreeding values can lead to reduced biological fitness of specimens and inbreeding depression (Hedrick and Garcia-Dorado 2016; Kardos et al. 2016). Specifically, the inbreeding coefficient can have values ranging from 0, for individuals whose parents are unrelated, to 1, for individuals whose parents have a maximum degree of relatedness and are completely inbred. A favorable value of this coefficient can be considered when it is less than 0.1 (Frankham et al. 2017). Finally, it has become clear in recent years that genetic load should also be considered when planning reinforcement programs (Mathur and DeWoody 2021; Bertorelle et al. 2022). For example, large populations tend to have a higher genetic load, so their use as source populations could be detrimental (van der Valk et al. 2019; Kyriazis et al. 2021). It is therefore essential to consider different genomic parameters in conservation plans.

In this work, we used a reduced representation (ddRAD) genomic sequencing technique (Peterson et al. 2012) to assess the genetic structure, diversity, and inbreeding of the capercaillie. This technique enables hundreds to thousands of single nucleotide polymorphisms (SNPs) to be obtained along the whole genome, making it possible to analyze the genetic structure at a higher resolution than with mitochondrial genes or microsatellites, as well as to characterize inbreeding levels of individuals, as demonstrated in other species (Querejeta et al. 2016; Escoda et al. 2017; Lavretsky et al. 2019). Samples for this work were obtained throughout the mountains of southern Europe and the western part of the boreal region to determine the evolutionarily significant units (Moritz 1994; Funk et al. 2012) of the species in Europe. In addition, we included a larger representation of the Pyrenean and Cantabrian populations to estimate

<figure><img src="image_5.png" /></figure>
inbreeding coefficients in these populations, which are
among the most endangered of the species.

**Materials and methods**

**Samples, DNA extraction and PCR amplification**

We used samples from 66 capercaillie specimens collected between 1987 and 2020 (Table S1 in the Supporting Information and Fig. 1). All the samples used were either tissue or blood and had been well preserved in ethanol or frozen in order to build a high quality genomic reference for the capercaillie. The samples were available from different institutions and scientists, and had been collected for previous studies (see Acknowledgements section).

The DNA was purified using the DNEasy Blood and Tissue Kit (QIAGEN), following manufacturer's instructions and previously described details (Escoda et al. 2019). The amount of DNA in each sample was measured by quantitative PCR (qPCR) using primers that amplify a conserved chromosomal fragment in vertebrates, as previously described (Escoda et al. 2017).

From the DNA extract, a 385 bp fragment of the
mitochondrial D-loop gene was amplified and sequenced
using capercaillie-specific primers (Duriez et al. 2007),
including extraction and PCR controls to rule out possible
contamination, and with PCR conditions as in previous
studies (Escoda et al. 2019). Both DNA extractions and

DNA amplification reactions of the D-loop and genomic
libraries were performed in separate rooms under
maximum sterile conditions, including irradiation with
ultraviolet light, to avoid contamination.

**Analysis of D-loop sequences**

Mitochondrial D-loop sequences were aligned using MAFFT v7.130b (Katoh and Standley 2013) and a maximum-likelihood phylogenetic tree was reconstructed with RAxML v8.0.19 (Stamatakis 2014) using a GTR model of nucleotide substitution and a gamma distribution of evolutionary rates. Mid-point rooting was used to represent the tree.

**Preparation and analysis of ddRAD libraries**

A ddRAD library (Peterson et al. 2012) was constructed
from the DNA, as detailed in previous studies (Escoda et al.
2017, 2019). Briefly, a minimum of 200 ng of DNA from
each sample was digested with two restriction enzymes
(EcoRI and MspI). Adapters were then ligated with a
specific 5 bp sequence (barcode) for each sample to allow
its subsequent identification. The DNA fraction between
300 and 400 base pairs was then selected on an agarose
gel. The DNA thus selected for each sample was amplified
in three PCR reactions in order to obtain a homogeneous
representation of the amplified genomic fragments. The
three PCRs from all the samples were pooled into a single
genomic library, which was sent for sequencing to the

<figure><img src="image_6.png" /><figcaption>**Fig. 1** Map of capercaillie samples analyzed by ddRAD sequencing, with the species range according to different sources shown in gray (Robles et al. 2006; Coppes et al. 2015), and indicating the main geographical origin of the samples. The distribution of the species con-</figcaption></figure>

tinues into the eastern part of the Asian continent. Some of the points
represent more than one sample. The shapefile of the countries was
obtained from Eurostat (https://ec.europa.eu/eurostat)
Genomic Sequencing Service at Pompeu Fabra University, using an Illumina NextSeq system, and generating single reads of 150 base pairs.

The reads obtained were processed with the *Stacks* 1.48 package (Catchen et al. 2011) using the same methods as in previous studies (Escoda et al. 2017, 2019). After removing the barcode, 145 bp sequences remained. Only loci assembled in at least 90% of the individuals were considered (*r*=0.9) for exporting sequence loci and SNPs. A minimum allele frequency (MAF) filter of 0 was used, except where indicated.

## Population structure analysis

Population structure analyses were performed with the SNPs dataset of a homogeneous group of capercaillies. All the European samples not coming from the Iberian Peninsula (n=28) were included in this group but, in order to maintain sampling homogeneity across different populations, 8 samples were selected from the Pyrenees together with 8 from the Cantabrian Mountains, so that this analysis included 44 samples in total, as indicated in Table S1. The dataset contained 8,627 SNPs. The *Structure* 2.3.4 program (Pritchard et al. 2000) was used with the correlated model and 300,000 generations of calculation after a burn-in of 30,000 generations. The number of populations (*K*) considered ranged from 2 to 6. A total of 10 runs per *K* were performed (only one run did not converge properly, at *K*=5, and this was removed as all the others converged to the same result). The Evanno method (Evanno et al. 2005) was used to identify the most supported number of populations from the *ΔK* graph. The uncorrelated model of *Structure* was also used for comparison. In addition, the program was tested with a 6,740 SNPs dataset from which singletons were removed (with MAF=0.0125), as suggested in other works (Linck and Battey 2019). Population structure was also analyzed with the program PopCluster (Wang 2022) using the 8627 SNPs dataset and *K* ranging from 1 to 6.

Principal component analysis (PCA) was performed on the 8627 SNPs dataset of the homogeneous group of capercaillies using the *KING* 2.2.5 program (Manichaikul et al. 2010).

## Genome diversity analysis

Genomic heterozygosity was calculated from the sequences for all the samples as the number of heterozygous positions divided by the total length of the genomic sequence analyzed, which consisted of 2,320,800 base pairs of genomic DNA sequence obtained from 16,040 loci. The significance of heterozygosity differences between populations was calculated using the Tukey's test (*p* < 0.05) in the R environment. Nucleotide diversity (π) was

calculated with the program *populations* from the *Stacks* package (Catchen et al. 2011) using the variant positions in the six capercaillie populations obtained from the *Structure* analysis.

## Application of an isolation-with-migration model

To estimate population split times, effective population sizes and migration rates in an isolation-with-migration (IM) framework, we used the program IMa3 version 1.12 (Hey et al. 2018), with similar conditions as in a previous work (Balmori-de la Puente et al. 2022). We considered the six geographical populations of Fig. 1, which were also supported by the *Structure* analysis: Pyrenees, Cantabrian Mountains, Alps, Balkans, Carpathians, and Western boreal zone. Only individuals of the homogeneous group (Table S1) were used. In addition, three individuals that were genetically very distant from their respective populations in the PCA analysis were excluded, as their presence in these populations may be due to recent migration events. The population topology used, based on a TreeMix analysis without migration (Pickrell and Pritchard 2012), was: ((Pyrenees, Cantabrian), (Alps, Balkans, (Carpathians, Boreal)))). Alternative topologies were also tested but did not converge properly. The final analysis was performed with 100 randomly selected loci, as analysis with a higher number of loci did not converge, but we checked that different sets of randomly chosen alignments produced similar results. The infinite sites model was used for all loci. We chose the priors based on their convergence in preliminary analyses: maximum split time *t*<sub>μ</sub>=0.05; maximum population size mutation rate 4N<sub>μ</sub> (q)=0.5; and maximum migration rate *m*/<sub>μ</sub>=2. Runs were performed with 480 chains on 60 processors. The burn-in was set to 500,000 steps and 15,000 genealogies were saved, sampling every 100 steps. The parameter estimates reported were the histogram bins with the highest value and the confidence intervals (C.I.) were the 95% highest posterior density intervals. The population tree figure was generated with the IMfig program (Hey 2010). The significance of population migration rates (2Nm, effective number of migrant genes per generation) was based on a log-likelihood-ratio test (Nielsen and Wakeley 2001). To scale parameters to demographic units, we used a mutation rate of 1.91 × 10<sup>-9</sup> mutations/site/year, determined for the chicken genome using four-fold degenerate sites (Nam et al. 2010). The chicken belongs to the same family, Phasianidae, as the capercaillie, so it is likely that they have similar mutation rates. In addition, we have shown in a previous work that four-fold degenerate sites have similar mutation rates to ddRAD loci (Balmori-de la Puente et al. 2022). The per-base mutation rate was multiplied by 145 to obtain the mutation rate of the loci. The generation time used to calculate
the effective population size was 3.8 years (Leclercq and Ménoni 2018).

## Relatedness and inbreeding coefficient

Individual inbreeding coefficients and pairwise relatedness coefficients were estimated from the SNPs dataset with the R library *related* 1.0 (Pew et al. 2015), using the full nine-states identity-by-descent model. Confidence intervals (95%) of the relatedness coefficients were calculated using 100 bootstraps over loci and these coefficients were considered significant if the lower limit of the confidence interval was greater than 0. Estimates of inbreeding and relatedness coefficients require information on the background allele frequencies of the population, so these were only calculated for the Cantabrian and Pyrenean populations, where the higher number of samples made this estimation possible.

## Results

### Mitochondrial phylogeography

A 385 bp fragment of the mitochondrial D-loop gene was sequenced for the 66 capercaillie individuals, which resulted in 15 different haplotypes. The phylogenetic tree constructed with these sequences showed the presence of two major mitochondrial clades (Fig. 2), one containing the majority of Pyrenean and Cantabrian individuals (Iberian group) and another corresponding to the rest of the European individuals (European group). As observed in previous studies (Duriez et al. 2007; Bajc et al. 2011; Klinga et al. 2015), several haplotypes from Balkan samples, specifically three from Bulgaria, were found to be associated with the Iberian clade, while some European haplotypes were found in the Pyrenees.

### Genetic structure of capercaillie populations

The *Structure* analysis was performed using 8627 SNPs from the homogeneous set of 44 samples (Table S1). The most supported number of populations (*K*), according to the Evanno method (Evanno et al. 2005), was 2 (Fig. 3a). This analysis showed that the capercaillie samples were firstly divided into an Iberian group and a group comprising the rest of the individuals from the European continent (Fig. 3b and c). Within the European group, all the specimens from the Alps presented a substantial level of admixture, since their genomic composition consists of approximately 75% of the European group with the rest corresponding to the Iberian group. The Pyrenean capercaillies also showed some admixture, with around 5% of the European component. Higher *K* values divided the European group into subgroups

<figure><img src="image_7.png" /><figcaption>**Fig. 2** D-loop-based mitochondrial phylogeny of capercaillie individuals, with different colors indicating the Iberian specimens and the rest of the European specimens. The root of the tree was placed at the midpoint. The scale is in substitutions/position</figcaption></figure>

associated with the large European mountain ranges. With *K*=3, the Alpine and Balkan populations formed one group while the Carpathian and western boreal populations formed another. With *K*=5, the European groups were separated into the following subgroups: Alps, Balkans, Carpathians and western boreal zone. Finally, with *K*=6 the Iberian group was divided into its two geographical populations.

Different conditions known to affect the *Structure* analysis in some circumstances were tested. When we used an uncorrelated frequency model, the most supported number of populations with the Evanno analysis was again 2 and the results with this *K* value were similar to those obtained with the correlated model. However, the components that appeared with values of *K*>2 were residual, giving rise to a structure similar to that of *K*=2 (Figure S1a). The analysis after filtering out singletons again supported *K*=2 and gave similar results to those obtained with the full SNPs dataset (Figure S1b). These results indicated that the Evanno result is robust under different conditions.

The population structure analysis performed with PopCluster was highly consistent with that of *Structure* with the correlated model from *K*=2 to *K*=5 (Figure S2), although at *K*=6 the two Iberian groups were not split,
<figure><img src="image_8.png" /><figcaption>**Fig. 3** a Estimated log likelihood of the data for the 10 runs performed for each number of populations (*K*) in the *Structure* analysis and the Δ*K* graph for estimating the best supported *K* according to the Evanno method. b *Structure* analysis of the homogeneous group</figcaption></figure>

but rather the Alpine population was split into two groups. In addition, the most supported number of populations according to this method was *K*=5, which differs from the *Structure* result.

of capercaillie samples with a correlated frequency model and *K* from 2 to 6. Each bar represents an individual and the colors indicate each of the inferred genomic groups. c Results of the *Structure* analysis with *K*=2 represented on the map

A PCA based on the SNP allele frequencies of the homo-
geneous set of capercaillies, in which the first two compo-
nents explain 15.1% and 11.3% of the variance, respectively,
arranged most of the samples in a similar way to their geo-
graphical distribution (Fig. 4). This analysis supported the
<figure><img src="image_9.png" /><figcaption>**Fig. 4** Principal component analysis of the capercaillie specimens. The main populations and some samples that appear separated from their respective populations are indicated. The percentage of variance captured by each principal component is indicated on the corresponding axis</figcaption></figure>

existence of two large groups separated along the principal axis (PC1), corresponding to the Iberian and European groups. Within each group, there was a separation into subgroups that approximately corresponded to the geographical separation between them, as is the case of the separation observed between the Iberian populations (Pyrenees and Cantabrian Mountains) on the one hand, and between the European populations (Alps, Balkans, Carpathians and western boreal zone), on the other. Among the capercaillie individuals that showed the greatest geographical discordance, the one from Bosnia was very far from the other Balkan specimens, and those from Ukraine and Romania were far from the other specimens from the Carpathians.

## Genomic diversity

The mean individual heterozygosity for the 66 capercaillies was 549 SNPs/Mb (Fig. 5a). The box plot of heterozygosity grouped by population and the statistical test performed showed that the differences between the European populations were not significant, whereas all differences between the European and Iberian populations and between the Cantabrian and Pyrenean populations were significant (Fig. 5b). A mean heterozygosity of 670 SNPs/Mb was calculated for the European group, 527 for the Pyrenean population and 366 for the Cantabrian population. A similar trend was observed for π: 0.16, 0.12 and 0.10 for the European, Pyrenean, and Cantabrian groups, respectively.

## Isolation-with-migration model

The population tree resulting from the application of the IM analysis is shown in Fig. 6. The main population divergence in the tree, that between the Iberian populations and the rest of the European populations, was estimated to have occurred 36 Kyr ago (C.I.: 26–58 Kyr), while the minimum divergence estimated was that between the two Iberian populations, 4 Kyr ago (C.I.: 1–15 Kyr) (Table S2). Estimates of effective population sizes were much smaller for the two Iberian populations, with 2 thousand individuals each (C.I. of 0–5 and 0–8 thousand individuals for the Cantabrian and Pyrenean populations, respectively), than for the European populations, with the Alpine population reaching 47 thousand individuals (C.I.: 25–92 thousand individuals; Table S2). Migration rates between different populations were detected. Although none was significant, 2Nm values higher than 0.1 were found between different populations and the Alpine population (Fig. 6).

## Inbreeding and relatedness coefficients in the Pyrenees and the Cantabrian Mountains

For the 14 capercaillies analyzed from the Cantabrian Mountains, the mean inbreeding coefficient was 0.10, with values ranging from 0 to 0.34 (Fig. 7a and b and Table S3). Inbreeding coefficients were lower for the 24 capercaillies from the Pyrenees, with a mean of 0.04 and values ranging from 0 to 0.11 (Fig. 7c and d and Table S3). In both populations, a slight increase in the inbreeding coefficient was observed over the years (Figs. 7b and 7d).

With respect to the relatedness coefficients, very few significant relationships were found and, in general, of a low degree, both in the Pyrenees and in the Cantabrian Mountains, so these could not be used to build connectivity maps. It is likely that the low number of kinship relationships detected is due to the wide range of years over which the capercaillies were sampled (Fig. 7b and d).

## Discussion

### Population structure and evolutionarily significant units in the capercaillie

The genetic structure analyses performed here were based on 8627 SNPs and methods that allow the genetic variability to be partitioned into its main groups, such as *Structure* and PCA, commonly used to delineate conservation units (Funk et al. 2012; Barbosa et al. 2018). Most results from these analyses indicated that the European capercaillie specimens fall into two major groups, in basic agreement with previous mitochondrial work (Duriez et al. 2007; Segelbacher and
<figure><img src="image_10.png" /><figcaption>**Fig. 5** a Individual genomic heterozygosity (in SNPs/Mb) of the capercaillies analyzed. The individuals have been ordered according to their populations of origin. b Box plot of heterozygosity grouped by population</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 6** Population tree derived from the application of an isolation-with-migration model. The five split times between populations are depicted as solid horizontal lines, with estimated values and the confidence intervals in Kyr on the left. Arrows represent migration rates (2Nm) higher than 0.1. The width of the boxes is proportional to the estimated effective population sizes, whose values and confidence intervals are shown in Table S2</figcaption></figure>

Piertney 2007; Bajc et al. 2011). Neither the most supported *Structure* analysis, which favored two main groups, nor further subdivisions, which essentially corresponded to geographical areas, are equivalent to the much more complex

map of the currently recognized subspecies, mostly dissociated from these geographical areas (Leclercq and Ménoni 2018; de Juana and Kirwan 2019). The two main genetic groups found in the *Structure* analysis can be considered

<figure><img src="image_12.png" /></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 7** Inbreeding coefficients of capercaillie from the Cantabrian Mountains represented **a** on a map and **b** as a function of year; and from the Pyrenees represented **c** on a map and **d** as a function of year. Sample numbers can be found in Table S1 in the Supporting Information</figcaption></figure>

as the two evolutionarily significant units (Moritz 1994; Funk et al. 2012) into which the species is divided: the Iberian unit on the one hand, and the rest of the populations of the species on the other hand, which we call the European unit. With respect to the Iberian unit, the *Structure* analyses showed that it was subdivided into distinct populations only at *K* = 6. In contrast, the European unit was divided into different subgroups with smaller *K* values (Alps, Balkans, Carpathians and western boreal zone). However, these subdivisions are not as strongly supported as the first one, according to the Evanno method (Fig. 3a). In addition, there is a high level of admixture between the European subgroups, indicating gene flow between them. The PCA analysis corroborated the existence of these two main genetic groups. The separation that can be observed between the subgroups in this plot could be explained by a relatively recent separation between them and genetic drift.

We also performed an IM analysis, which allowed us to obtain preliminary estimates of basic population parameters, such as population divergence times, effective population sizes, and migration rates or gene flow between populations, and visualize them in a population tree (Fig. 6). The results of this model were consistent with and supported the various analyses of population genetic structure and diversity that we performed. The high divergence observed in the IM analysis between the Iberian populations and the others explains their early separation in the *Structure* analysis (Fig. 3b) and the optimal *K* = 2 determined by the Evanno method

(Fig. 3a), as well as the primary separation into these two groups observed in the PCA (Fig. 4). Using a mutational rate inferred from the chicken genome, the IM model estimated that this divergence happened 36 Kyr ago (C.I.: 26–58 Kyr). This value has a wide confidence interval associated with it, but it suggests that this split would have occurred during the last glacial period. On the opposite site of the tree, the most recent divergence involved the Cantabrian and Pyrenean populations and took place 4 Kyr ago (C.I.: 1–15 Kyr). The uncertainty is also large for this estimate, but it allows us to infer that this split would have occurred during the Holocene. The relatively low divergence between the Iberian populations explains why they only split at *K* = 6 in the *Structure* analysis (Fig. 3b) and are very close in the PCA (Fig. 4). The gene flow inferred by the IM model from the Iberian and Carpathian populations to the Alps explains the admixture observed in this population at different *K* values (Fig. 3b). Finally, the smaller population sizes estimated by IM for the Iberian populations, especially for the Cantabrian one, agree with the measures of individual heterozygosity (Fig. 5) and genetic diversity (π) carried out. However, it should also be noted that, due to the computational complexity required to estimate the parameters of the IM model, only 100 loci of 145 bp could be used in the analyses, which contributes to the high variance of some parameters, particularly the most recent divergence times, and the lack of statistical support for the migration rates. Using longer fragments derived from complete genomes would probably improve the accuracy of
the results. Nevertheless, the IM analysis performed here
with ddRAD data shows the great potential of using large-
scale genomic data to estimate critical population parameters
in conservation genomics.

It can be hypothesized that both evolutionary units were
isolated in different refugia during the glaciations. Following
Leclercq and Ménoni (2018), one of these would have been
an Iberian refugium while the other populations could have
been located in various southeastern refugia, probably with
great interconnection. The split time of 36 Kyr (C.I.: 26–58
Kyr) estimated by IM for the two evolutionary units, during
the last glacial period, agrees with this view. The isolation
of these two groups for long periods of time could have
originated specific adaptations in each of them, meaning
that the integrity of the two groups should be protected in
order to maintain their evolutionary potential (Moritz 1994;
Funk et al. 2012). However, specific analyses should be
performed in future work, particularly with whole genomes,
to understand whether there are important adaptations in
these two groups (Barbosa et al. 2020; Hohenlohe et al.
2021).

As in previous mitochondrial works (Duriez et al. 2007;
Bajc et al. 2011; Klinga et al. 2015), we found two main
haplogroups, mostly corresponding to an Iberian group
and a European group (Fig. 2). As also reported in these
works, we found capercaillie samples from Bulgaria whose
mitochondrial haplotypes clustered with those from the
Pyrenees (Fig. 2). However, both *Structure* and IM analyses
(Figs. 3 and 6) showed no genomic affinity between the
Pyrenean and Balkan populations. The presence of these
Pyrenean-type haplotypes in the Balkans could be explained
by an evolutionary model involving some kind of past
connectivity between the two populations (Leclercq and
Ménoni 2018), perhaps quite ancient; however, it did not
leave a significant footprint in the nuclear genome of the
Balkan capercaillies. The discordances between nuclear
and mitochondrial genes found in other species have been
explained by sex-biased dispersal and the particular mode
of mitochondrial gene transmission (Toews and Brelsford
2012), although specific analyses would be necessary to
determine the most likely model for explaining the presence
of Pyrenean haplotypes in the Balkan capercaillies. In
any case, it seems clear from the genomic results that all
Balkan populations belong to the European evolutionarily
significant unit and not to the Iberian unit.

The *Structure* analysis revealed an important admixture
between the two capercaillie evolutionary units, especially
in the Alps (Fig. 3). This admixture is most likely not a
very recent admixture event, as it affects all analyzed
individuals in a similar way. This admixture could be
explained by past gene flow between the Iberian and Alpine
populations, possibly when intermediate populations of the
Massif Central in France were more abundant, as probably

occurred until the great deforestation of the eighteenth and
nineteenth centuries (Williams 2006). Indeed, the IM model
infers gene flow from the Iberian populations into the Alps
(Fig. 6). A lower degree of admixture detected by *Structure*
in the Pyrenean capercaillies, all of which have a small
European component, suggests that some gene flow also
occurred towards the Pyrenees. The presence of European
mitochondrial haplotypes in the Pyrenees is consistent with
this hypothesis. However, IM did not detect gene flow from
the Alps to the Pyrenees, so it could have been small.

In general, the patterns observed in the *Structure* analysis (Fig. 3), the PCA (Fig. 4), and the IM analysis (Fig. 6) were consistent with the geographic location of the specimens, something to be expected for a species with limited dispersal such as the capercaillie (Leclercq and Ménoni 2018; de Juana and Kirwan 2019). However, there are three discordant specimens, which can be appreciated in the PCA (Fig. 4): one from Ukraine, one from Bosnia, and one from Romania, which are genetically far from the populations where they were sampled. Among the various potential explanations, the introduction of capercaillies from the boreal zone into the Carpathians and the Balkans could be behind these mismatches. This interpretation is supported by the fact that these three samples have a large genomic component from the boreal zone, as deduced from the *Structure* analysis with K≥4 (Fig. 3). However, a microsatellite study of a large number of capercaillies from the Carpathians showed that many individuals from the Eastern Carpathians (Romania and Ukraine) had a boreal component, so a common history of Eastern Carpathian and boreal capercaillie populations cannot be excluded (Klinga et al. 2015). As knowing the origin of these samples requires further work, we excluded these three capercaillies from the IM analysis to ensure that the model only reflects the evolutionary history of the species and not possible recent events.

Genetic diversity and inbreeding in the Iberian
capercaillie populations

The genomic heterozygosity results (Fig. 5) point to a low
genetic diversity in the Pyrenean (on both the southern and
northern slopes) and, especially, the Cantabrian population,
compared to those in the rest of Europe, as reported in
previous studies based on microsatellites (Alda et al. 2011;
Vázquez et al. 2012). However, none of the heterozygosity
values obtained are typical of the most genetically
impoverished populations, which usually have values
lower than 100 SNPs/Mb (Bravo et al. 2021; Escoda and
Castresana 2021). In any case, the lower genetic diversity
of the Pyrenean and, particularly, Cantabrian capercaillie
populations reflects the smaller effective population size in
the Iberian Peninsula, as corroborated by the IM analysis,
and indicates the greater vulnerability of these populations.

<figure><img src="image_14.png" /></figure>
In the case of the Pyrenean population, the extensive deforestation in previous centuries (Williams 2006) may have contributed to the decline of this bird in the Pyrenees at the beginning of the twentieth century (Inqram 1915). The even lower genetic diversity in the Cantabrian Mountains can be explained by the reduction of the population over recent decades (Jimenez et al. 2022), without ruling out that the decline could have started even earlier.

The inbreeding coefficients recorded in the Pyrenees seem to be relatively favorable, with most of them being below 0.1, and a maximum of 0.11. There is, however, an upward trend in this coefficient in the specimens sampled since 2005 compared to those sampled in previous decades, although unfortunately the number of samples is not sufficient to prove the statistical significance of this pattern (Fig. 7d). In the case of the Cantabrian population the situation is worse, since several inbreeding coefficient values higher than 0.2 were recorded, indicating that subpopulations in this area are isolated to the degree that mating occurs between close relatives. Similar to the trend noted in the Pyrenees, the highest inbreeding coefficient values in the Cantabrian Mountains were also observed in recent years (Fig. 7b), suggesting a greater genetic erosion of both populations in recent decades, in line with recent demographic studies (Gil et al. 2020; Bal et al. 2021; Canut et al. 2021; Jimenez et al. 2022). In fact, the highest inbreeding coefficient values in the Cantabrian Mountains were found in the eastern subpopulation, where all 3 individuals analyzed presented very high values (Fig. 7b). This subpopulation is now practically extinct (Jimenez et al. 2022). We cannot determine the causes of this local extinction and, in fact, it was not necessarily caused by high levels of inbreeding. However, a thorough monitoring of this parameter could help to detect early deterioration of genetic health and connectivity in any subpopulation, so that mitigating measures could be rapidly implemented. For example, assisted migration could be carried out using individuals from nearby subpopulations that are currently unconnected (Supple and Shapiro 2018; Hohenlohe et al. 2021). In addition, kinship-based connectivity analyses (Escoda et al. 2017, 2019) could be conducted to help gain further insight into the dispersal patterns of the subpopulations that remain more isolated.

## Implications for the conservation of capercaillie populations and future research perspectives

This work represents a first step in the use of genomic information to provide relevant information for the conservation and management of capercaillie populations. The analyses of the genetic structure of the European capercaillie populations based on ddRAD sequences and SNPs suggest that the species is divided into two

evolutionarily significant units: the Iberian unit and the European unit. Nevertheless, in this study we were not able to analyze samples from Asia, so that additional genetic groups in that continent cannot be ruled out. According to the present results, the Cantabrian and Pyrenean capercaillie populations belong to the Iberian evolutionary unit. On the other hand, the populations from the Balkans, Carpathians, Alps and western boreal area belong to the European evolutionary unit. The management of populations should take this structure into account and avoid exchanges between evolutionary units in order to preserve the evolutionary potential of the species. Regarding the possibility of exchanges within an evolutionary unit, populations with the highest genetic affinities should be preferentially used as source populations for genetic rescue to avoid potential problems of outbreeding depression (Funk et al. 2012; Supple and Shapiro 2018; Hohenlohe et al. 2021). Genetic health monitoring showed relatively low levels of inbreeding coefficients for the Pyrenean population, although there seems to be a worrying trend towards an increase in this coefficient in the samples obtained in recent years. This result is consistent with the slow decline and deterioration of reproductive success in this population (Gil et al. 2020; Bal et al. 2021; Canut et al. 2021). In the Cantabrian population, some individuals were found to have very high inbreeding coefficients, especially in the eastern subpopulation, which is now practically extinct (Jimenez et al. 2022). However, there are still subpopulations and individuals with low inbreeding coefficients in the western Cantabrian Mountains. Management plans should favor the use of capercaillie specimens with the lowest inbreeding coefficients and other indicators considered, or at least avoid using the most inbred specimens in different actions.

Our work based on ddRAD has made a significant advance in capercaillie conservation genomics, especially for some of the most threatened populations. However, these studies should be continued in the future using whole genomes and methods that can handle high-density SNPs. These analyses will provide new genomic inferences that were not possible with ddRAD, such as inbreeding measures based on runs of homozygosity, which will add a temporal framework to inbreeding, molecular adaptations that should be conserved, or estimates of genetic load, which have recently been shown to be essential to inform captive breeding, assisted migration, and reinforcement plans (Barbosa et al. 2020; Hohenlohe et al. 2021). All the genetic health information derived from genomes will be crucial to better understand capercaillie conservation problems and to support management plans to try to recover the most threatened populations.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s10592-023-01567-6.
**Acknowledgements** We thank the following people for providing samples: Florian Kunz (Institute of Wildlife Biology and Game Management, Austria), Gabriel de Pedro (Junta de Castilla y León, Spain), Iván Afonso (Conselh Generau d'Aran, Spain), Javier Quesada (Banc de Teixits, Museu de Ciències Naturals de Barcelona, Spain), Jordi Solá (Gobierno de Andorra), Laura Kvist, Jouni Aspi (University of Oulu, Finlandia), María Suárez (Tragsatec, Spain), Pedro de la Cruz (Spain), Rainer Ploner (Südtiroler Landesverwaltung, Land- und Forstwirtschaft, Italy), Rubén Moreno-Opo (Ministerio para la Transición Ecológica y el Reto Demográfico, Spain), the Mountain Galiform Observatory of the French Pyrenees (France), Rudi Suchant and Annette Kohnen (Forstliche Versuchs- und Forschungsanstalt Baden-Württemberg, Germany), Tomasz Strzała (Wroclaw University of Environmental and Life Sciences, Poland), Petar Zhelev (Forestry University, Bulgaria), and Dalibor Ballian (Forestry Faculty, Bosnia). We are also grateful to the following institutions for facilitating samples and support: Servicio de Espacios Naturales, Flora y Fauna de la Dirección General de Patrimonio Natural y Política Forestal de la Junta de Castilla y León, Dirección General del Medio Natural del Principado de Asturias, and Servei de Fauna i Flora de la Generalitat de Catalunya.

**Author contributions** JP, LE and JC conceived the ideas and designed methodology; JP, LP, KF and EM collected the data; LE and JC analyzed the data; LE and JC led the writing of the manuscript. All authors contributed to the final version of the manuscript.

**Funding** Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. This work was commissioned and funded by Ministerio para la Transición Ecológica y el Reto Demográfico, Junta de Castilla y León, and Gobierno del Principado de Asturias. Funding was also provided by the research project PID2020-113586GB-100 of MCIN/AEI/10.13039/501100011033. The funders promoted and financed the work, but had no role in the analysis, decision to publish, or preparation of the manuscript.

**Data availability** ddRAD sequencing data have been deposited in Dryad (https://doi.org/10.5061/dryad.mkkwh714j) and mitochondrial sequences in GenBank (OQ378102-OQ378167).

## Declarations

**Conflict of interest** The authors have no relevant financial or non-financial interests to disclose.

**Ethics approval** No animals were specifically handled or captured for this work and, therefore, this study did not require an ethics approval by a specific committee.

**Consent for publication** All authors have approved the manuscript for publication.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Alda F, Sastre P, De La Cruz-Cardiel PJ, Doadrio I (2011) Population genetics of the endangered Cantabrian capercaillie in northern Spain. Anim Conserv 14:249-260. https://doi.org/10.1111/j.1469-1795.2010.00425.x

Bajc M, Čas M, Ballian D, Kunovac S, Zubić G, Grubešić M, Zhelev P, Paule L, Grebenc T, Kraigher H (2011) Genetic differentiation of the Western capercaillie highlights the importance of South-eastern Europe for understanding the species phylogeography. PLoS One 6:e23602. https://doi.org/10.1371/journal.pone.0023602

Bal G, Bacon L, Ménoni E, Calenge C, Millon A, Besnard A (2021) Modélisation de la dynamique du grand tétras des Pyrénées françaises pour sa gestion adaptative. Technical report MNHN et OFB. https://www.researchgate.net/publication/354150226Modelisation_de_la_dynamique_du_grand_tetras_des_Pyrenees_francaises_pour_sa_gestion_adaptative

Balmori-de la Puente A, Ventura J, Miñarro M, Somoano A, Hey J, Castresana J (2022) Divergence time estimation using ddRAD data and an isolation-with-migration model applied to water vole populations of *Arvicola*. Sci Rep 12:4065. https://doi.org/10.1038/s41598-022-07877-y

Barbosa S, Mestre F, White TA, Pauperio J, Alves PC, Searle JB (2018) Integrative approaches to guide conservation decisions: using genomics to define conservation units and functional corridors. Mol Ecol 27:3452-3465. https://doi.org/10.1111/mec.14806

Barbosa S, Hendricks SA, Funk WC, Rajora OP, Hohenlohe PA (2020) Wildlife population genomics: Applications and approaches. In: Hohenlohe P, Rajora OP (eds). Population genomics: Wildlife. Springer Nature Switzerland AG, Cham. pp. 3-59. https://doi.org/10.1007/13836_2020_83

Bertorelle G, Raffini F, Bosse M, Bortoluzzi C, Iannucci A, Trucchi E, Morales HE, van Oosterhout C (2022) Genetic load: genomic estimates and applications in non-model animals. Nat Rev Genet 23:492-503. https://doi.org/10.1038/s41576-022-00448-x

BirdLife I (2016) *Tetrao urogallus*, Western capercaillie. The IUCN Red List of Threatened Species. https://doi.org/10.2305/IUCN.UK.2016-3.RLTS.T22679487A85942729.en

Bravo GA, Schmitt CJ, Edwards SV (2021) What have we learned from the first 500 avian genomes? Annu Rev Ecol Evol Syst 52:611-639. https://doi.org/10.1146/annurev-ecolsys-012121-085928

Canut J, Garcia-Ferré D, Alfonso-Jordana I (2021) Gall fer *Tetrao urogallus*. In: Franch M, Herrando S, Anton M, Villero D, Brotons L (eds) Atles dels ocells nidificants de Catalunya: Distribució i abundància 2015-2018 i canvi des de 1980. Institut Català d'Ornitologia i Cossetània Edicions, Barcelona, pp 112-113

Catchen JM, Amores A, Hohenlohe P, Cresko W, Postlethwait JH (2011) Stacks: building and genotyping loci de novo from short-read sequences. Genes Genom Genet 1:171-182. https://doi.org/10.1534/g3.111.000240

Coppes J, Kochs M, Ehrlacher J, Suchant R, Braunisch V (2015) The challenge of creating a large-scale capercaillie distribution map. Grouse News 50:21-23

de Juana E, Kirwan GM (2019) Western capercaillie (*Tetrao urogallus*). In: del Hoyo J, Elliott S, Sargatal J, Christie DA, de Juana E (eds). Handbook of the Birds of the World Alive. Lynx Edicions, Barcelona. Doi:https://doi.org/10.2173/bow.wescap1.01

Duriez O, Sachet J-M, Ménoni E, Pidancier N, Miquel C, Taberlet P (2007) Phylogeography of the capercaillie in Eurasia: what is the conservation status in the Pyrenees and Cantabrian Mounts? Conserv Genet 8:513-526. https://doi.org/10.1007/s10592-006-9165-2

Escoda L, Castresana J (2021) The genome of the Pyrenean desman and the effects of bottlenecks and inbreeding on the genomic
landscape of an endangered species. Evol Appl 14:1898-1913.
https://doi.org/10.1111/eva.13249

Escoda L, González-Esteban J, Gómez A, Castresana J (2017) Using relatedness networks to infer contemporary dispersal: Application to the endangered mammal *Galemys pyrenaicus*. Mol Ecol 26:3343-3357. https://doi.org/10.1111/mec.14133

Escoda L, Fernández-González A, Castresana J (2019) Quantitative analysis of connectivity in populations of a semi-aquatic mammal using kinship categories and network assortativity. Mol Ecol Resour 19:310-326. https://doi.org/10.1111/1755-0998.12967

Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol 14:2611-2620. https://doi.org/10.1111/j.1365-294X.2005.02553.x

Frankham R, Ballou JD, Ralls K, Eldridge M, Dudash MR, Fenster CB, Lacy RC, Sunnucks P (2017) Genetic management of fragmented animal and plant populations. Oxford University Press, Oxford. https://doi.org/10.1093/oso/9780198783398.001.0001

Funk WC, McKay JK, Hohenlohe PA, Allendorf FW (2012) Harnessing genomics for delineating conservation units. Trends Ecol Evol 27:489-496. https://doi.org/10.1016/j.tree.2012.05.012

Gil JA, Gómez-Serrano MÁ, López-López P (2020) Population decline of the capercaillie *Tetrao urogallus aquitanicus* in the Central Pyrenees. Ardeola 67:285-306. https://doi.org/10.13157/arla.67.2.2020.ra4

González MA, Fuertes B, Blanco-Fontao B, De Frutos Á (2017) The edge also matters: human threats in a rare Mediterranean habitat for Cantabrian capercaillie. Bird Conserv Internat 27:582-593. https://doi.org/10.1017/S0959270916000587

Hedrick PW, Garcia-Dorado A (2016) Understanding inbreeding depression, purging, and genetic rescue. Trends Ecol Evol 31:940-952. https://doi.org/10.1016/j.tree.2016.09.005

Hey J (2010) The divergence of chimpanzee species and subspecies as revealed in multipopulation isolation-with-migration analyses. Mol Biol Evol 27:921-933. https://doi.org/10.1093/molbev/msp298

Hey J, Chung Y, Sethuraman A, Lachance J, Tishkoff S, Sousa VC, Wang Y (2018) Phylogeny estimation by integration over isolation with migration models. Mol Biol Evol 35:2805-2818. https://doi.org/10.1093/molbev/mys162

Hohenlohe PA, Funk WC, Rajora OP (2021) Population genomics for wildlife conservation and management. Mol Ecol 30:62-82. https://doi.org/10.1111/mec.15720

Inqram C (1915) A few notes on *Tetrao urogallus* and its allies. Ibis 57:128-133. https://doi.org/10.1111/j.1474-919X.1915.tb08183.x

Jahren T, Storaas T, Willebrand T, Fossland Moa P, Hagen B-R (2016) Declining reproductive output in capercaillie and black grouse – 16 countries and 80 years. Anim Biol 66:363-400. https://doi.org/10.1163/15707563-00002514

Jimenez J, Godinho R, Pinto D, Lopes S, Castro D, Cubero D, Osorio MA, Pique J, Moreno-Opo R, Quiros P, Gonzalez-Nuevo D, Hernandez-Palacios O, Kery M (2022) The Cantabrian capercaillie: A population on the edge. Sci Total Environ 821:153523. https://doi.org/10.1016/j.scitotenv.2022.153523

Kämmerle J-L, Braunisch V, Suchant R, Coppes J (2020) Quantifying the effectiveness of habitat management to counter local extinction: A case-study on capercaillie. For Ecol Manag 474:118379. https://doi.org/10.1016/j.foreco.2020.118379

Kardos M, Taylor HR, Ellegren H, Luikart G, Allendorf FW (2016) Genomics advances the study of inbreeding depression in the wild. Evol Appl 9:1205-1218. https://doi.org/10.1111/eva.12414

Kardos M, Akesson M, Fountain T, Flagstad Ø, Liberg O, Olason P, Sand H, Wabakken P, Wikenros C, Ellegren H (2018) Genomic consequences of intensive inbreeding in an isolated wolf

population. Nat Ecol Evol 2:124-131. https://doi.org/10.1038/s41559-017-0375-4

Katoh K, Standley DM (2013) MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol 30:772-780. https://doi.org/10.1093/molbev/mst010

Klinga P, Mikoláš M, Zhelev P, Höglund J, Paule L (2015) Genetic differentiation of Western capercaillie in the Carpathian Mountains: the importance of post glacial expansions and habitat connectivity. Biol J Linn Soc 116:873-889. https://doi.org/10.1111/bij.12643

Kyriazis CC, Wayne RK, Lohmueller KE (2021) Strongly deleterious mutations are a primary determinant of extinction risk due to inbreeding depression. Evol Lett 5:33-47. https://doi.org/10.1002/evl3.209

Lavretsky P, DaCosta JM, Sorenson MD, McCracken KG, Peters JL (2019) ddRAD-seq data reveal significant genome-wide population structure and divergent genomic regions that distinguish the mallard and close relatives in North America. Mol Ecol 28:2594-2609. https://doi.org/10.1111/mec.15091

Leclercq B, Ménoni E (2018) *Le grand tétras*. Biotope Éditions.

Leroy G, Carroll EL, Bruford MW, DeWoody JA, Strand A, Waits LP, Wang J (2018) Next-generation metrics for monitoring genetic erosion within populations of conservation concern. Evol Appl 11:1066-1083. https://doi.org/10.1111/eva.12564

Linck E, Battey CJ (2019) Minor allele frequency thresholds strongly affect population structure inference with genomic data sets. Mol Ecol Resour 19:639-647. https://doi.org/10.1111/1755-0998.12995

Manichaikul A, Mychaleckyj JC, Rich SS, Daly K, Sale M, Chen W-M (2010) Robust relationship inference in genome-wide association studies. Bioinformatics 26:2867-2873. https://doi.org/10.1093/bioinformatics/btq559

Mathur S, DeWoody JA (2021) Genetic load has potential in large populations but is realized in small inbred populations. Evol Appl 14:1540-1557. https://doi.org/10.1111/eva.13216

Moreno-Opo R, Afonso I, Jiménez J, Fernández-Olalla M, Canut J, García-Ferré D, Piqué J, García F, Roig J, Muñoz-Igualada J, González LM, López-Bao JV (2015) Is it necessary managing carnivores to reverse the decline of endangered prey species? Insights from a removal experiment of mesocarnivores to benefit demographic parameters of the Pyrenean capercaillie. PLoS One 10:e0139837. https://doi.org/10.1371/journal.pone.0139837

Moritz CC (1994) Defining “Evolutionarily Significant Units” for conservation. Trends Ecol Evol 9:373-375. https://doi.org/10.1016/0169-5347(94)90057-4

Nam K, Mugal C, Nabholz B, Schielzeth H, Wolf JB, Backstrom N, Kunstner A, Balakrishnan CN, Heger A, Ponting CP, Clayton DF, Ellegren H (2010) Molecular evolution of genes in avian genomes. Genome Biol 11:R68. https://doi.org/10.1186/gb-2010-11-6-r68

Nielsen R, Wakeley J (2001) Distinguishing migration from isolation: a Markov chain Monte Carlo approach. Genetics 158:885-896. https://doi.org/10.1093/genetics/158.2.885

Peterson BK, Weber JN, Kay EH, Fisher HS, Hoekstra HE (2012) Double digest RADseq: an inexpensive method for de novo SNP discovery and genotyping in model and non-model species. PLoS One 7:e37135. https://doi.org/10.1371/journal.pone.0037135

Pew J, Muir PH, Wang J, Frasier TR (2015) related: an R package for analysing pairwise relatedness from codominant molecular markers. Mol Ecol Resour 15:557-561. https://doi.org/10.1111/1755-0998.12323

Pickrell JK, Pritchard JK (2012) Inference of population splits and mixtures from genome-wide allele frequency data. PLOS Genet 8:e1002967. https://doi.org/10.1371/journal.pgen.1002967

Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155:945-959. https://doi.org/10.1093/genetics/155.2.945
Querejeta M, González-Esteban J, Gómez A, Fernández-González A, Aymerich P, Gosálbez J, Escoda L, Igea J, Castresana J (2016) Genomic diversity and geographical structure of the Pyrenean desman. Conserv Genet 17:1333–1344. https://doi.org/10.1007/s10592-016-0865-y

Quevedo M, Bañuelos MJ, Obeso JR (2006) The decline of Cantabrian capercaillie: how much does habitat configuration matter? Biol Conserv 127:190–200. https://doi.org/10.1016/j.biocon.2005.07.019

Robles L, Ballesteros F, Canut J eds. (2006) *El urogallo en España, Andorra y Pirineos franceses. Situación actual* (2005). SEO/Bird-Life, Madrid.

Rodríguez-Muñoz R, Mirol PM, Segelbacher G, Fernández A, Tregenza T (2007) Genetic differentiation of an endangered capercaillie (*Tetrao urogallus*) population at the southern edge of the species range. Conserv Genet 8:659–670. https://doi.org/10.1007/s10592-006-9212-z

Rodríguez-Muñoz R, Valle CR, Bañuelos M-J, Mirol P (2015) Reveal-ing the consequences of male-biased trophy hunting on the main-tenance of genetic variation. Conserv Genet 16:1375–1394. https://doi.org/10.1007/s10592-015-0747-8

Ryder OA (1986) Species conservation and systematics: the dilemma of subspecies. Trends Ecol Evol 1:9–10. https://doi.org/10.1016/0169-5347(86)90059-5

Scott PA, Allison LJ, Field KJ, Averill-Murray RC, Shaffer HB (2020) Individual heterozygosity predicts translocation success in threatened desert tortoises. Science 370:1086–1089. https://doi.org/10.1126/science.abb0421

Segelbacher G, Piertney S (2007) Phylogeography of the European capercaillie (*Tetrao urogallus*) and its implications for con-servation. J Ornithol 148:269–274. https://doi.org/10.1007/s10336-007-0153-1

Stamatakis A (2014) RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30:1312–1313. https://doi.org/10.1093/bioinformatics/btu033

Steiner CC, Putnam AS, Hoeck PEA, Ryder OA (2013) Con-servation genomics of threatened animal species. Annu Rev Anim Biosci 1:261–281. https://doi.org/10.1146/annur-ev-animal-031412-103636

Storch I (ed) (2007) Grouse: Status Survey and Conservation Action Plan 2006–2010. IUCN, Fordingbridge

Storch I, Bañuelos M-J, Fernández-Gil A, Obeso JR, Quevedo M, Rodríguez-Muñoz R (2006) Subspecies Cantabrian capercaillie *Tetrao urogallus cantabricus* endangered according to IUCN criteria. J Ornithol 147:653–655. https://doi.org/10.1007/s10336-006-0101-5

Supple MA, Shapiro B (2018) Conservation of biodiversity in the genomics era. Genome Biol 19:131. https://doi.org/10.1186/s13059-018-1520-3

Toews DPL, Brelsford A (2012) The biogeography of mitochondrial and nuclear discordance in animals. Mol Ecol 21:3907–3930. https://doi.org/10.1111/j.1365-294X.2012.05664.x

van der Valk T, de Manuel M, Marques-Bonet T, Guschanski K (2019) Estimates of genetic load in small populations suggest extensive purging of deleterious alleles. bioRxiv 13. Doi:https://doi.org/10.1101/696831

Vázquez JF, Pérez T, Quirós F, Obeso JR, Albornoz J, Domínguez A (2012) Population genetic structure and diversity of the endan-gered Cantabrian capercaillie. J Wildl Manag 76:957–965. https://doi.org/10.1002/jwmg.346

Wang J (2022) Fast and accurate population admixture inference from genotype data from a few microsatellites to millions of SNPs. Heredity 129:79–92. https://doi.org/10.1038/s41437-022-00535-z

Williams M (2006) Deforesting the Earth: From Prehistory to Global Crisis, an Abridgment. University of Chicago Press, Chicago, United States. https://doi.org/10.7208/chicago/9780226899053.001.0001

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.