Genetic assessment reveals inbreeding, possible hybridization, and low levels of genetic structure in a declining goose population

Abstract The population numbers of taiga bean goose (Anser fabalis fabalis) have halved during recent decades. Since this subspecies is hunted throughout most of its range, the decline is of management concern. Knowledge of the genetic population structure and diversity is important for guiding management and conservation efforts. Genetically unique subpopulations might be hunted to extinction if not managed separately, and any inbreeding depression or lack of genetic diversity may affect the ability to adapt to changing environments and increase extinction risk. We used microsatellite and mitochondrial DNA markers to study the genetic population structure and diversity among taiga bean geese breeding within the Central flyway management unit using non‐invasively collected feathers. We found some genetic structuring with the maternally inherited mitochondrial DNA between four geographic regions (ɸ ST = 0.11–0.20) but none with the nuclear microsatellite markers (all pairwise F ST‐values = 0.002–0.005). These results could be explained by female natal philopatry and male‐biased dispersal, which completely homogenizes the nuclear genome. Therefore, the population could be managed as a single unit. Genetic diversity was still at a moderate level (average H E = 0.69) and there were no signs of past population size reductions, although significantly positive inbreeding coefficients in all sampling sites (F IS = 0.05–0.10) and high relatedness values (r = 0.60–0.86) between some individuals could indicate inbreeding. In addition, there was evidence of either incomplete lineage sorting or introgression from the pink‐footed goose (Anser brachyrhynchus). The current population is not under threat by genetic impoverishment but monitoring in the future is desirable.


| INTRODUC TI ON
Knowledge of population genetic structure is essential for guiding management and conservation of species. The presence of highly divergent subpopulations with low levels of gene flow could warrant a status of separate management units (MUs; Palsbøll et al., 2007).
Loss of genetic diversity and inbreeding depression, processes that especially affect small and fragmented populations, may lead to a loss of evolutionary potential and contribute to a higher extinction risk (Frankham, 2005). Harvesting can also lead to adverse genetic changes such as alteration of population subdivision (extirpation of local populations), loss of genetic variation, and selective genetic changes (reduction in certain phenotypes targeted by hunting) that may further compromise population viability (Allendorf et al., 2008).
Thus, knowledge of genetic structure and diversity is especially important for managing harvested species and subspecies to ensure sustainable hunting. This has been our motivation to study the taiga bean goose (Anser fabalis; Figure 1), which is hunted over most of its range but has suffered a marked decline during recent decades (Fox et al., 2010). However, very little is currently known about the genetic population structure and diversity of the taiga bean geese in their breeding area.
Most Holarctic goose populations are on the increase (Fox & Leafloor, 2018) to the extent that this is causing conflicts with humans, especially in agriculture (Fox & Madsen, 2017). However, a few goose populations are declining, one of them being the taiga bean goose (Fox et al., 2010;Fox & Leafloor, 2018) breeding in Northern Europe and Asia ( Figure 2). The population size of the taiga bean goose has nearly halved from the 90,000-100,000 individuals in the 1990s (Nilsson et al., 1999) to the current estimate of 52,000 individuals for the total wintering population size (Fox & Leafloor, 2018). Although the exact cause of decline is uncertain, overharvesting is one plausible explanation. Since the taiga bean goose population is still open to hunting, sustainable management of the population is of crucial importance. Within this framework, an International Single Species Action Plan (ISSAP) was developed by The African-Eurasian Migratory Waterbird Agreement (AEWA) to conserve the taiga bean goose (Marjakangas et al., 2015). Hunting of the Northeast/Northwest European population of the taiga bean goose can still be continued within the limits of agreed sustainable use within the ISSAP framework (Marjakangas et al., 2015).
The ISSAP for the taiga bean goose recognizes four subpopulations or flyway management units: Western, Central, Eastern 1, and Eastern 2 subpopulations (Marjakangas et al., 2015; Figure 2). These flyway units are distinguished by isotopic composition of feathers collected from bean geese from Western, Central, and Eastern 1 flyway units  but the units have not been confirmed genetically and it is not known if these units should be further subdivided into smaller management units. For effective management of the taiga bean goose, population genetic structure should be assessed in order to preserve genetic diversity. Previously, the population genetic structure of the taiga bean goose has only been studied within a limited area in Central Scandinavia belonging to the Western flyway management unit (de Jong et al., 2019).
The taiga bean geese are elusive, especially during the breeding period, and observing even neck-banded geese in breeding areas is rarely possible due to the long escape distance of geese (Pirkola & Kalinainen, 1984a)-although camera traps have been shown to be a viable alternative for taiga bean goose monitoring on breeding grounds (Nykänen et al., 2021). Thus, catching geese and sampling blood would be very difficult and stressful for the birds.
Developments in non-invasive genetic assessments have allowed sampling of elusive and endangered species without the need to handle or even observe the animals . For birds, nest material and molted feathers provide a valuable source of DNA (Pearce et al., 1997;Segelbacher, 2002). The taiga bean geese perform a molt from the middle of June to the middle of August, and during this flightless period, they spend time in the wettest part of mires or in the vicinity of ponds, where they leave abundant cues of their presence such as tracks, signs of grazing, feces, and molted feathers (Pirkola & Kalinainen, 1984b). In our study, we utilized this source of non-invasive feather samples collected by volunteers (citizen science) to study the genetic population structure of the breeding taiga bean geese.
Our aims were (1) to evaluate the utility of a citizen-science approach in collecting genetic samples for the elusive taiga bean goose, (2) to assess the broad-scale genetic structure among breeding taiga bean geese in Finland (consisting of a large part of the Central management unit), (3) to estimate genetic diversity within this population, (4) to evaluate the current effective population size and demographic fluctuations, and (5) to investigate if taiga bean geese hybridize with another subspecies of bean goose breeding in Europe, the tundra bean goose (Anser fabalis rossicus), and with the most closely related species, the pink-footed goose (Anser brachyrhynchus). For this, we used highly variable microsatellites first to F I G U R E 1 Taiga bean goose (Anser fabalis fabalis) sampled for this study during ringing. Photo: Tuomo Turunen identify individuals and parentage and then to study variation at the population level. In addition, we defined the genetic structure of maternal lineages by sequencing the most variable part of the mitochondrial control region. This region was also used to verify the subspecies (see Honka et al., 2017).

| Study species
The taxonomy of the bean goose and the related pink-footed goose (A. brachyrhynchus) has been controversial (Ruokonen & Aarvak, 2011). Currently, the bean goose is either split into two species (A. fabalis and Anser serrirostris; Sangster & Oreel, 1996) or treated as one species (A. fabalis) with four subspecies (Anser fabalis fabalis, A. f. rossicus, Anser fabalis middendorffii, and Anser fabalis serrirostris; Mooij & Zöckler, 1999). Based on a recent study using genome-wide data, A. f. fabalis and A. f. rossicus should be classified as subspecies (Ottenburghs et al., 2020) and thus we follow here the classification of one species with four subspecies.
From here on, we focus on the western subspecies and refer to A. f. fabalis as the taiga bean goose and A. f. rossicus as the tundra bean goose.
The taiga bean goose breeding distribution covers the forested taiga area from Scandinavia to Western Siberia (Scott & Rose, 1996; Figure 2) and in Finland the "aapa" mire zone (see Laitinen et al., 2007) with core breeding area extending from Lapland to Northern Ostrobothnia (Pirkola & Kalinainen, 1984b; Figure 3). In Finland, the taiga bean goose is listed as vulnerable (the Red List of Finnish species; Lehikoinen et al., 2019) with the breeding population size estimated to be 1700-2500 pairs (i.e., 3400-5000 breeding individuals; Valkama et al., 2011). Subadults and failed breeders perform a molt migration to Novaya Zemlya Piironen et al., 2021) and thus are not counted. In addition to the taiga bean goose, a very small number of tundra bean geese may breed in

| Sampling and DNA extraction
Bean goose feathers were collected from nests or brood-rearing/  (Piironen et al., 2021), and their molted feathers were not collected in this study. The samples collected in the present study are hereafter referred to as the "Finnish population." We used mostly a citizen-science approach-in which the public is involved in scientific research-to collect feather samples. As single feathers cannot be aged reliably, our sample might include adults, juveniles, and goslings with unknown proportions. We also created outgroups for our Finnish population by obtaining samples of: (1) Swedish breeding taiga bean geese (2014; n = 7), (2) migrating Russian taiga bean geese hunted from Finland (Honka et al., 2017) or Estonia (2010; n = 7), (3) Norwegian breeding tundra bean geese (2002, 2006; n = 7), and (4) Iceland breeding pink-footed geese (Ruokonen et al., 2005, the Natural History Museum of Reykjavik; n = 7) ( Figure 3). All feathers were stored in paper envelopes at room temperature prior to DNA extraction. Sampling of feathers was performed in a laboratory in which no PCR products are handled.
DNA from a calamus and from a blood clot (when visible ;Horváth et al., 2005)  F I G U R E 2 Breeding and wintering distributions of the taiga (Anser fabalis fabalis) and tundra (Anser fabalis rossicus) bean goose shown with different colors, and the autumn migration routes of (a) the taiga bean goose shown with green arrows and (b) the tundra bean goose shown with blue arrows. The flyway management units used in the International Single Species Action Plan (Marjakangas et al., 2015) for the taiga bean goose are also shown with dashed ellipses. Maps redrawn from BirdLife International (2018) and Marjakangas et al. (2015) (a) (b) Also, feathers belonging to other bird species based on morphology (e.g., feather shape typical of common crane Grus grus) were omitted from DNA extraction. The extracted DNA was stored in −20°C.

| Microsatellite genotyping
We genotyped the samples for 28 microsatellite loci (Kleven et al., 2016;Noreikiene et al., 2012) in four multiplex reactions (A-D) ( We performed fragment analysis with an ABI 3730 and scored the alleles using GeneMapper 5 (Applied Biosystems).
We used a stepwise approach to amplify the microsatellite panels. First, we amplified the samples with panel B (five loci) to screen for quality and to reduce the number of duplicate individuals, as follows. We performed individual identification by visual inspection of the genotypes and by using the "Regroup genotypes" option in the program Gimlet v.1.3.3. (Valière, 2002). If one sampling site included several identical genotypes, these were assumed to be replicates of the same individual and all except one was excluded. Next, we amplified panel A (seven loci) from the apparently good-quality samples selected in the first step. We performed the individual identification similarly as after the panel B, as genotyping errors could have been interpreted as different individuals in the first step and again excluded all but one identical genotype (taking into account the genotyping errors). Panels C and D were amplified with the samples selected in the second step and the individual identification was performed once more in order to test for any remaining errors in the data and again kept only one of the identical genotypes.
Due to the inherent presence of allele dropouts and false alleles in microsatellite data , we regarded multilocus genotypes with up to three mismatches to belong to the same individual. Mismatches in 4-5 loci between multilocus genotypes were carefully checked and if the pattern was consistent with allele dropout due to poor-quality template DNA, the samples were excluded.
Consensus genotypes for each individual were created manually by choosing the majority genotype for each locus or, if only two replicates were genotyped from an individual, the heterozygote, as allele dropouts were more likely than false alleles (see Table S1). In addition, individuals with more than 25% missing data were excluded, resulting with 491 individuals. Three individuals were excluded from the analysis of the taiga bean geese because these samples had an mtDNA sequence of a different subspecies or were possibly still in spring migration (see below). Thus, 488 individuals were included in analyses unless otherwise stated. We calculated the unbiased probability of identity (P ID ) and probability of identity of siblings (P ID SIB )  from the identified individuals using the Gimlet program.
We performed the PCR reactions and fragment analyses similarly as

F I G U R E 3
The breeding distribution of the taiga bean goose (Anser fabalis fabalis) in Finland and the approximate location of the sampling sites of the current study. The exact sampling sites are not shown in order to protect the breeding locations and the sites closer than 16 km were merged. Samples are divided here into four geographical regions shown with different symbols and colors. Samples of tundra bean goose (Anser fabalis rossicus) mtDNA are indicated with a star and one dead goose that was not ascertained as breeding is indicated with a black dot. The sampling sites of outgroups are also shown. The background breeding distribution map was created by interpolating the breeding category indices (1 = unlikely breeding, 2 = possible breeding, 3 = probable breeding, and 4 = confirmed breeding) of the

| Microsatellite analyses
The data were divided into four geographical regions ( Figure 3) and all analyses were performed on these groups. We used the pro-  ) were calculated using the program GenAlEx 6.503 (Peakall & Smouse, 2006, and allelic richness (A R ) was calculated using the program FSTAT 2.9.4 (Goudet, 1995). Departure from the Hardy-Weinberg equilibrium (HWE) and the degree of linkage disequilibrium (LD) were determined using the program Genepop 4.7.0 (Rousset, 2008) with a sequential Bonferroni correction (Rice, 1989) applied to these tests. We evaluated the power of the microsatellite markers to detect a signal of population differentiation by using the simulation program POWERSIM 4.1 (Ryman & Palm, 2006).
We inferred parentage and sibships using the program Colony 2.0.6.5 (Jones & Wang, 2010) by setting "monogamy" for both females and males (as geese are known to form stable pair bonds) and performed three iterations of the long run with inbreeding model with a full-likelihood method. All taiga bean goose individuals (n = 488) were placed in the candidate offspring category, females (n = 237) and geese of unknown sex (n = 62) (based on molecular sexing) were placed in the candidate mother category (total n = 299), and males (n = 189) and geese with unknown sex (n = 62) (based on molecular sexing) were placed in the candidate father category (total n = 251). We excluded each individual from being the mother or father of itself. We used the calculated ADO and FA rates (Table   S1) as the marker error rates. We calculated relatedness (r) between goose dyads using the program ML-relate (Kalinowski et al., 2006) taking account of the null alleles. As it has been shown that sampling of close relatives biases genetic structure analyses, especially in the program Structure (Anderson & Dunham, 2008;Rodríguez-Ramilo & Wang, 2012), we created a dataset from which the inferred parents, all but one sibling and individuals with r > 0.55, were removed resulting in a subset of data with 376 non-kin individuals. All analyses were performed on this non-kin dataset except for F ST (fixation index), ɸ ST , and N e (effective population size) that were calculated with the full dataset because the precision of these genetic estimates may suffer with the purging of all siblings (Waples & Anderson, 2017). In addition, analyses in the program Structure (Pritchard et al., 2000) were performed for both the non-kin and the full dataset. Due to complex nature of our data collected by a citizen-science approach without knowledge of the identity of the samples-thereby potentially including multiple cohorts and a family structure-our data should be considered non-random and family correlated.
We inferred population structure using the program Structure 2.3.4. (Falush et al., 2003;Hubisz et al., 2009;Pritchard et al., 2000) with the LOCPRIOR option, with individuals within 16 km treated as coming from one location. We used an admixture ancestry model and correlated allele frequencies with a burn-in of 100,000 and a run length of 1,000,000. Five iterations were performed with the possible number of clusters (K) set from 1 to 10. The ad hoc approach of Evanno et al. (2005) was used to infer the most likely number of K clusters in the data as implemented in the program Structure Harvester (Earl & vonHoldt, 2012). We used the program Clumpak 1.1. (Kopelman et al., 2015) to visualize the Structure results and to create consensus among the different iterations. In addition, we used Discriminant Analysis of Principal Components (DAPC; Jombart et al., 2010) implemented in the R package (R Core Team, 2018) "adegenet" 2.1.1 (Jombart, 2008;Jombart & Ahmed, 2011) to assess population structure. The lowest number of non-kin samples was in the Western Finland region (n = 50) and, accordingly, we randomly chose an equal number of samples in each geographical region (n = 50 in each population, total n = 200) as the DAPC could be biased due to unequal samples sizes. Cross-validation (xval command) was used to evaluate the number of PCs to retain. DAPC was also used to compare the outgroup samples and the Finnish population-by performing DAPC first on the outgroups and then importing the "unknown" Finnish population into this same framework using the "predict" function. Pairwise F ST -values were calculated using the ENA correction implemented in the program FreeNA due to the possible presence of null alleles in the data. To test for isolation by distance, we performed a Mantel test implemented in the R package (R core team, 2018) "ade4" 1.7.13 (Bougeard & Dray, 2018;Chessel et al., 2004; using Euclidean distances. We also performed a spatial autocorrelation analysis (Smouse & Peakall, 1999) using the program GenAlEx to test for fine-scale geographic patterns. The significance of the analysis was tested using a heterogeneity test (Smouse et al., 2008).
In addition, the autocorrelation analysis was performed on females and males separately to estimate sex biases in relatedness (Banks & Peakall, 2012).
The effective population size (N e ) was estimated with a linkage disequilibrium model (Hill, 1981;Waples, 2006;Waples & Do, 2010), assuming monogamy and using 0.05 as the critical value for allele frequency as implemented in the program NeEstimator v2.1 (Do et al., 2014). The N e estimate was compared with the estimate based on sibship assignment (Wang, 2009) (Dib et al., 1996;Ellegren, 2000;Sun et al., 2012) and a generation time of 5-7.5 years (Dillingham, 2010 We studied hybridization between the taiga bean goose and either the pink-footed goose (taiga bean goose x pink-footed goose) or the tundra bean goose (taiga bean goose x tundra bean goose) by using a simulation study. First, we simulated 100 pure parental individuals (selected by q > 0.99 by a preliminary NewHybrids run) and 100 individuals in different hybrid classes (F1, F2, and backcrosses) using the program HybridLab1.0 (Nielsen et al., 2006). We used these simulated individuals as an input for the program NewHybrids 1.1 (Anderson & Thompson, 2002) with Jeffreys-like priors for both mixing proportions and allele frequencies, setting burn-in to 20,000 sweeps and the chain length to 100,000 MCMC sweeps with the z option. We compared these results to the Structure program run with the admixture ancestry model and correlated allele frequencies with burn-in set to 10,000 and run length to 100,000.
Five iterations were performed with the number of K set to 2. As the Structure program performed better than the NewHybrids program (Figures S1-S4), we secondly run Structure analysis with the simulated parentals and the Finnish population including also the two tundra bean geese and the plausible pink-footed goose individual (n = 491). In this second analysis, the aim was to search for admixture with other subspecies or species in the Finnish population.

| Mitochondrial DNA analyses
The hypervariable portion of the mitochondrial control region domain I (210 bp) was amplified using primers AdCR1F and AdCR2R (Honka et al., 2018) from the individuals identified with microsatellite genotyping (n = 491). The primers were designed to contain mismatches to Numts (nuclear sequences of mitochondrial origin; Lopez et al., 1994), which are problematic in genetic studies if not accounted for (Bensasson et al., 2001;Sorenson & Quinn, 1998 et al., 2000), and greylag geese (Anser anser) as an outgroup (AF159961; Ruokonen et al., 2000) using the program BioEdit 7.2.5 (Hall, 1999). A median joining network (Bandelt et al., 1999) was constructed using the program PopART (Leigh & Bryant, 2015), and number of haplotypes (H), haplotype (h) and nucleotide (π) diversities, Tajima's D (D), and Fu's Fs (Fs) were calculated using program DnaSP v. 6.12 (Rozas et al., 2017). The presence of genetic structure on individual and region levels was tested by an analysis of variance framework using analysis of molecular variance (AMOVA), which is based on hierarchical variance of gene frequencies. We calculated pairwise ɸ ST -values between regions and performed an AMOVA analysis using the program Arlequin 3.5.2.2. (Excoffier & Lischer, 2010) with the substitution model of Jukes-Cantor (Jukes & Cantor, 1969). This was the best substitution model based on Akaike information criteria (AIC) and Bayesian information criterion (BIC) values in the program MEGA X (Kumar et al., 2018) and is supported in the Arlequin program. A sequential Bonferroni correction (Rice, 1989) was applied to the ɸ ST statistics.

| Microsatellite genotyping success
Microsatellite genotyping was successful for 1805 feather samples (at least one multiplex panel successful) with a success rate of 84%.
After omitting possible redundancy, this amounted to a total of 491 individuals in our dataset. One individual that may have been on spring migration and had pink-footed goose mtDNA (individual found dead, see Figure 3) was omitted from the Finnish breeding population analyses but included in the mtDNA haplotype network, maps, and hybridization analyses. Similarly, we also omitted two tundra bean goose individuals, leaving us with a final total of 488 taiga bean goose individuals. The error rate in the whole dataset (28 loci) was 0.019 per allele over all loci and 0.039 per locus over all loci.
The mean ADO rate was 0.041, the mean FA rate was 0.009, and the mean null allele rate was 0.034 (Table S1). The cumulative unbiased probability of identity and the cumulative probability of identity of siblings were low (P ID = 5e−25, P ID SIB = 4e−10; Table S1), indicating that our markers can separate different individuals with high confidence.
No evidence of large allele dropout was present in the loci, but the program MicroChecker suggested that some loci show stuttering. However, the stuttering was not consistent over the geographical regions (see Figure 3 for the regions). We removed the loci Abra14 and Abra29 due to low polymorphism, Abra9 due to a high rate of missing data (43%) and Abra15 and Abra43 due to a high frequency of null alleles (>10%) ( Table S1). Locus Afa18 was not in Hardy-Weinberg equilibrium in any of the four geographical regions and was removed. Therefore, a final total of 22 loci were kept for further analyses. Several other loci showed deviations from Hardy-Weinberg equilibrium and some were in linkage disequilibrium after Bonferroni correction, but these were not consistent between regions, so these loci were kept in the analyses. Using these 22 loci, their frequency distributions, and the numbers sampled in each region, the simulation program POWSIM indicated that the power of the markers is fairly high. For example, a true F ST of 0.001 would be detected with a probability of 100%, indicating that these loci are suitable for detecting significant divergence, had it occurred.

| Molecular sexing and relatedness
Molecular sexing using the HINTZ/W gametologs showed congruent results with the cloacal examination except for one goose out of seven. Based on the cloaca, this goose was assigned as a female while based on molecular sexing it was a male. This sample was replicated and the same result was obtained by molecular sexing.
We identified 237 females and 189 males in the Finnish population which was significantly different from equality (χ 2 = 5.40, P = 0.02).
The sex could not be determined for 62 individuals due to lack of amplification of this marker.
Most of the geese dyads studied were unrelated (

| Genetic diversity
The number of alleles varied between 171 and 194 with a mean of 184 over the geographical regions (Table 1). Allelic richness, which takes into account the sample size, was almost the same in all regions (A R = 7.4-7.6; Table 1). The Eastern Finland/Kainuu region had the largest number of private alleles (PA = 9), followed by Lapland (PA = 7;

| Population structure
The Structure analyses indicated little difference between the regions (mapped in Figure 3) or indeed between any individuals within the Finnish breeding taiga bean goose population ( Figure S5a), with some structure introduced when using the full dataset including kin ( Figure S5b). The optimal K was 4 based on the ad hoc statistics ΔK, but the log likelihood values did not reach a clear plateau ( Figure S6), TA B L E 1 Summary statistics for the 376 non-kin taiga bean geese (Anser fabalis fabalis) genotyped by 22 microsatellite loci and grouped by four geographical regions (see Figure 3) instead all samples were assigned to all clusters with high admixture proportions indicating that the K would actually be 1.
The DAPC analyses showed that slight genetic structuring might be present between Western Finland and Lapland regions, as the samples mostly did not overlap (Figure 4a). Samples from Northern Ostrobothnia/Southern Lapland and Eastern Finland/Kainuu region, on the other hand, overlapped with the other two regions, although several samples from Northern Ostrobothnia/Southern Lapland were distinct from the other populations (Figure 4a). The DAPC analysis performed to the outgroups showed that all outgroups (Swedish taiga bean geese, Russian taiga bean geese, tundra bean geese, and pink-footed geese) clustered to their own groups except for minor overlap between Swedish and Russian taiga bean geese (Figure 4b).
When the Finnish samples were fitted to this framework, most of the samples clustered with the Russian taiga bean geese. However, some samples showed genetic affinity to Swedish taiga bean geese, tundra bean geese, and to a smaller extent, to the pink-footed geese ( Figure 4b).
Pairwise F ST -values were always very low, with mean values between regions of less than 0.005 ( Table 2). The Western Finland region was the most differentiated from all of the other regions, although none of the values were statistically significant.
We found no isolation-by-distance pattern using the Mantel test (r = 0.019, P = 0.18; Figure S7) or the spatial autocorrelation analysis ( Figure 5a). No sex-specific differences were evident in the spatial autocorrelation analysis either (Figure 5b).

| Effective population size and bottlenecks
The effective population size (N e ) of the full dataset including kin was estimated to be 1127.6 individuals (95% confidence intervals 3790-5876) ( Figure S8).

| Hybridization
The NewHybrids and the Structure analysis based on simulated data showed that identification of different hybrid categories (F1, F2, and backcrosses) between the taiga bean and the pink-footed goose and between the taiga bean and the tundra bean goose is difficult using microsatellite markers due to the extent of variation within each hybrid category ( Figures S1-S4). However, the pure simulated parentals were readily identifiable. We also run a analysis showed that the simulated pure taiga bean geese and pink-footed geese could be separated and that some Finnish taiga bean geese showed admixture with pink-footed geese (Figure 6a).
No pure pink-footed geese were present in our dataset according to this analysis. The analysis for the simulated taiga and tundra bean geese showed similar results. Some of the Finnish taiga bean geese showed admixture with tundra bean geese (Figure 6b).
The two tundra bean geese sampled from Finland in this study were identified as pure tundra bean geese with >0.99 probability (Figure 6b, arrows).  Figure  3) with non-kin individuals in parentheses (n = 375). Pairwise ɸ ST -values for mitochondrial data (n = 447) above diagonal with non-kin individuals in parentheses (n = 343). Statistically significant values after Bonferroni correction indicated with an asterisk.

F I G U R E 5
Correlogram from spatial autocorrelation analysis for (a) nonrelated taiga bean geese (Anser fabalis fabalis; n = 376) and (b) females and males separately (n = 328). The autocorrelation coefficient (r) was plotted against the function of distance class (50 km). The 95% confidence interval (dashed lines, U = upper limit, L = lower limit) was determined by 999 permutations with the null hypothesis of no population structure, and 95% error bars were determined by bootstrap resampling of 1000 replicates The most common haplotype among the Finnish bean geese was Fa3 (n = 261), followed by the haplotype FAB1a/FAB1b/Fa2 (n = 159).
The slashes between haplotype names denote identical haplotypes based on the sequenced region. These haplotypes, however, differ based on the whole control region (Honka et al., 2017;Ruokonen et al., 2008). The two most common haplotypes were distributed throughout Finland (Figure 8a (Table 4). In addition, 93.8% of the total variation was within region variation and 6.2% was among region variation (ɸ ST = 0.062; P < .001) according to AMOVA.

| DISCUSS ION
We did not detect clear population structure within the Finnish breeding taiga bean geese using microsatellite markers. All analyses suggested close to a panmictic population, except the DAPC that indicated slight structuring among Western Finland, Northern Ostrobothnia/Southern Lapland, and Lapland (see Figure 3).
Presence of geographically localized mtDNA haplotypes and higher ɸ ST -values for females than in males, however, suggested at least some maternal genetic structure. It was unforeseen to find such little genetic structure within such a large geographic area, but the F I G U R E 7 Median joining haplotype network for the Finnish breeding bean geese, different bean goose subspecies (Anser fabalis fabalis, Anser fabalis rossicus, Anser fabalis serrirostris, and Anser fabalis middendorffii), and pink-footed goose (Anser brachyrhynchus). An mtDNA sequence of greylag goose (Anser anser) was used to root the haplotype network. The sizes of the circles are proportional to the frequency of each haplotype and tick marks across branches indicate the number of mutational differences. Forward slashes between haplotype names denote identical haplotypes based on the sequenced fragment but differing based on the whole control region sequence pairing system of geese, in which pair formation occurs already in common wintering or spring staging areas, can explain these results.
We found moderate genetic diversity and signs of inbreeding within the Finnish taiga goose population.
Surprisingly, we also found that a pink-footed goose mtDNA haplotype is widespread (although at low frequency) in the taiga bean goose. This could indicate hybridization between the taiga bean goose and the pink-footed goose and admixture was also evident in the microsatellite data. In addition, we confirmed breeding of tundra bean geese in the northernmost Lapland and the presence of a vagrant tundra bean goose in Southern Finland. The microsatellite data suggested introgression between the taiga and the tundra bean goose as well.

| Citizen science
This study proved that a citizen-science approach to feather collection was highly efficient for the elusive bean goose. About 89% of the samples were collected by a citizen-science approach, while the rest were collected by persons with an affiliation to a research institute. Even though bean goose feathers, especially the down feathers, can be mixed with other large birds in similar habitats, such as the common crane (Grus grus) and the whooper swan (Cygnus cygnus), only a few feathers had to be excluded as belonging to other species. This was determined either based on feather morphology or because no PCR product amplified from pristine feathers not identifiable based on morphology. Benefits of citizen-science approach for this study included the large number of samples and the broad geographic range of sampling.

| Relatedness, inbreeding, and genetic diversity
Interestingly, we identified one family in which one of the goslings was not the offspring of the social parents, thus representing a possible case of gosling adoption or intraspecific nest parasitism, behaviors which have been observed among other goose species (Zicus, 1981;Weigmann & Lamprecht, 1991;Choudhury et al., 1993;Larsson et al., 1995;Nilsson & Kampe-Persson, 2003;Anderholm et al., 2009; for review, see Kalmbach, 2006).
Overall, a large majority of the individual goose dyads were unrelated or related to a very low degree. However, positive inbreeding coefficients (F IS ) were detected and a few of the goose dyads showed high relatedness values (r = 0.60-0.86), indicating potential inbreeding. This was surprising as the detrimental effects of inbreeding to individual fitness, known as the inbreeding depression, are widely documented (Keller & Waller, 2002), and usually animals avoid inbreeding through several mechanisms (Pusey & Wolf, 1996). However, a single case of sibling pairing has been observed in the Canada goose (Branta canadensis; Lebeuf & Giroux, 2013). Heterozygote deficiency, that is, positive F IS could be due to Wahlund effect, which is caused by the merging of populations with different allele frequencies (Wahlund, 1928). Even though we detected only a very low level of population structure, it is possible that the population structure is very fine scaled or there is population structure in the wintering sites instead of breeding sites. We were unable to discern exact family relationships beyond parentoffspring and siblings, which might affect the Hardy-Weinberg equilibrium and cause positive F IS -values.
Nuclear genetic diversity appeared not to be reduced within

| Population structure and demography
We did not observe any population structure when using micro-  (Table 2). We did not observe isolation by distance either, further indicating lack of spatial genetic structuring ( Figure 5a). These findings indicate that the taiga bean goose population is close to panmictic.
On the other hand, some of the rarer mtDNA haplotypes were  The reason for this is unknown as Northern Ostrobothnia/Southern Lapland is a central region and thus there should not be restrictions for movement as might be the case for peripheral regions. This region is a core breeding area for the taiga bean goose due to a large number of aapa mires, so potentially this region is a more favorable breeding habitat and thus represents a source population. In addition, an unknown portion of the data could be from juvenile individuals, which might affect the data in this region. Also, we found a difference in haplotype composition between females and males as the FAB3 haplotype was only found in females and haplotype Fa7 only in males. An AMOVA analysis indicated that only about 6.2% of the total variation was between regions, thus mtDNA genetic structuring is still rather limited.
However, nuclear spatial genetic patterns did not indicate sexspecific differences ( Figure 5b). As opposed to male philopatry observed in most other birds, geese show female philopatry to natal areas (Greenwood, 1980;van der Jeugd et al., 2002). Therefore, genetic structure may be promoted especially in the maternally inherited mitochondrial DNA or the female-specific W chromosome, as was seen here as higher ɸ ST -values in females than in males. It seems, however, that the dispersal of male bean geese is so high that it completely homogenizes the nuclear genome. Geese pair in their wintering grounds or in the spring staging areas and males follow females to the female's natal area (Rohwer & Anderson, 1988), which allows geese from even distant breeding areas to pair and thus mediate gene flow in the biparentally inherited nuclear DNA.
Long-distance dispersal of a few males could lead to panmixia as only one migrant per generation in an ideal population is enough to prevent population differentiation due to drift (Mills & Allendorf, 1996;Wang, 2004, but see Vucetich & Waite, 2000. For example, in greylag and brent geese (Branta bernicla hrota), most males breed close to their natal site, but a minority of the geese undergo longdistance dispersal (Harrison et al., 2010;Nilsson & Persson, 2001).
On the other hand, a lack of sex-specific dispersal differences has Both the linkage disequilibrium and the sibship method produced similar estimates of the effective population size with 1128 and 1134 individuals, respectively. The Finnish population size is estimated to be 1700-2500 breeding pairs (i.e., 3400-5000 individuals) based on survey data (Valkama et al., 2011), thus excluding nonbreeders and juveniles. The ratio between the N e estimated based on microsatellite markers (this study) and the estimate of the number of breeding individuals is 0.23-0.33. Accordingly, N e is just about a quarter or a third of the estimated breeding population. Low N e /N ratios (0.11-0.14) are commonly reported among animals (Frankham, 1995;Palstra & Ruzzante, 2008) and based on N e = 1128, the total population size (N) would be thus 8057-10,255 individuals, a surprisingly large estimate. However, our estimate of N e from only Finnish samples is perhaps somewhat misleading due to a continuous population unrestricted by national borders.
We found no indication of past population bottlenecks, instead, we inferred population expansion starting around 22,000-32,000 years ago. This was unexpected since the taiga bean goose population is in decline (Fox et al., 2010;Fox & Leafloor, 2018

| Hybridization
We also compared the Finnish population with other taiga bean geese (Russian and Swedish), the tundra bean goose and the pinkfooted goose, using DAPC. As expected, the Finnish geese mostly grouped with the Russian geese (Figure 4b), which should belong to the same flyway management unit (Central; see Figure 2).
However, some Swedish geese (Western unit) showed genetic affinity to the Russian and Finnish geese and thus gene flow between the different flyways could be present. Some of the Finnish bean geese clustered closely with either the tundra bean goose or the pink-footed goose ( Figure 6). This could be due to a lack of resolution in our microsatellite markers to discriminate between different populations or an indication of possible hybridization between these populations. Whole-genome resequencing of the taiga and the tundra bean geese has shown that the genomes of these subspecies are homogenous expect for a few "islands of differentiation" due to extensive gene flow 60,000 years ago (Ottenburghs et al., 2020). Thus, the close affinity of the subspecies is probably due to the past hybridization. Two individuals were shown to be tundra bean geese based on tundra bean goose mtDNA haplotypes ( Figure 7) and by microsatellite markers (Figure 6). One of these individuals was sampled from northernmost Finland, confirming that this subspecies breeds in Finland (Figure 8a). The other tundra bean goose was sampled from a metropolitan area in Helsinki (Figure 8a,c), outside of the breeding range of either subspecies and was thus a vagrant bird.
In this study, we showed that the pink-footed goose mtDNA is widespread (although at low frequency, 4% of the studied population) in the Finnish taiga bean goose population (Figure 8a-c).
This was unexpected because the pink-footed goose breeds in Greenland, Iceland, and Svalbard and no breeding attempts have been recorded in Finland. Although the pink-footed goose can be found as a vagrant bird in Finland, and has started to regularly migrate through the Western Finland in recent decades (Heldbjerg et al., 2019), we do not believe that we have found 16 pink-footed geese, but instead admixed individuals. Even though the species do not share breeding grounds, they share common wintering areas in the United Kingdom, Denmark, and the Netherlands, which allow the species to come to contact. The microsatellite data supported admixture ( Figure 6), indicating possible hybridization and introgression between the taiga bean goose and the pink-footed goose. The amount of nuclear admixture seems to be rather high and even more than the 4% suggested by mtDNA. Geese show a high propensity for hybridization (for review, see , and in the genomics era, increasing number of studies have identified that ancient hybridization, adaptive introgression, and hybrid speciation are much more common than previously thought (Ottenburghs et al., 2017;Taylor & Larson, 2019). Interestingly, the mtDNA of the tundra bean goose had not introgressed into the taiga bean goose according to our results. The mtDNA of the pink-footed goose may convey adaptive benefits as hybridization often leads to the introgression of adaptive genetic variation (Arnold & Kunte, 2017), but also incomplete lineage sorting can explain the presence of the pink-footed goose mtDNA (Degnan & Rosenberg, 2009). The pinkfooted goose has been treated as a separate species as it has formed a monophyletic group based on mitochondrial DNA (Ruokonen et al., 2008), although more recently a sister species relationship was suggested between the tundra bean goose and the pink-footed goose (Ottenburghs, Megens, et al., 2016). Further studies are needed to elucidate the phylogenetic position of the pink-footed goose and the possible admixture scenario.

| Management implications
We did not find evidence to divide the Finnish bean goose population into smaller management units or subpopulations as there was no strong genetic structuring within Finland. Therefore, the flyway management units outlined in the International Single Species Action Plan (ISSAP; Marjakangas et al., 2015) seem to be justified based on our study (see Figure 2). The genetic diversity was found to be moderate and effective population size fairly large, thus the Finnish breeding taiga bean geese are not under immediate threat by genetic impoverishment. However, as the genetic diversity was lower compared to other bean goose subspecies and widespread goose species, further reductions in genetic diversity should be avoided to maintain the evolutionary potential of this subspecies, especially since inbreeding was detected in the population. Even though we found evidence of possible hybridization and introgression, naturally occurring hybridization does not pose a threat to populations as it is natural part of species evolution (Allendorf et al., 2001;Taylor & Larson, 2019). Thus, natural hybridization should not disqualify species from conservation programs and protection (vonHoldt et al.,

2018).
Non-invasive genetic sampling could be used for genetic monitoring of the taiga bean goose in the future to provide estimates of local population census sizes, survival, recruitment, and temporal variation in genetic diversity in order to ensure the genetic viability of the population (Schwartz et al., 2007). Besides microsatellites and mtDNA, SNP markers could be used in genetic monitoring, as SNPs do not require laborious calibration between different laboratories enabling flyway-wide monitoring, and are well suited for noninvasive samples as amplicon lengths are short and SNPs are less error prone than microsatellites (Carroll et al., 2018). A SNP panel has been already developed for subspecies identification for the greater white-fronted goose (Wilson et al., 2019) and similarly a SNP panel could provide a feasible alternative for the genetic monitoring for the taiga bean geese.

ACK N OWLED G M ENTS
We are grateful to the numerous volunteers who collected feathers and enabled this study. We also thank Risto Karvonen, Ingar Jostein Øien, Allan Hamari, and the Natural History Museum of Reykjavik for samples. Special thanks to Antti Piironen and the Pohjois-Karjalan metsähanhimiehet for feather samples from ringed geese. We are grateful for the help with field sampling from Toni Laaksonen, Mervi Kunnasranta, Petri Timonen, and Tuomas Seimola (The Natural Resource Institute of Finland) and the crew of HF Helicopters. We also thank everyone who helped advertise the feather collection with special thanks to Mikko Alhainen from the Finnish Wildlife Agency. We thank Kari Merikanto for coding the error rate calculation program. In addition, we thank the Finnish Cultural Foundation's North Ostrobothnia Regional Fund (grant number: 60182022) and the Finnish Cultural Foundation (grant number: 00200356) for providing funding for JH. We thank the two anonymous reviewers for their helpful comments for improving the article.

CO N FLI C T O F I NTE R E S T
None declared.