A genetic assessment of the human‐facilitated colonization history of black swans in Australia and New Zealand

Movement of species beyond their indigenous distribution can fundamentally alter the conservation status of the populations involved. If introductions are human‐facilitated, introduced species could be considered pests. Characterizing the colonization history of introduced species can therefore be critical to formulating the objectives and nature of wildlife management strategies. The black swan (Cygnus atratus) is native to Australia but is considered a reintroduced species in New Zealand, where the endemic population was reported extinct during the 19th century. After the reintroduction of a small number of individuals from Australia, the New Zealand population expanded unexpectedly rapidly, which was attributed to simultaneous waves of migration from Australia. An alternative, but hitherto unformalized, hypothesis is that local extant populations remained and admixed with introduced individuals. To contribute to our understanding of the reintroduction history of the species, we investigated dispersal patterns and demographic histories of seven populations from Australia and New Zealand, using population genetic inferences from a microsatellite dataset. Our results on genetic structure, dispersal rates, and demographic histories provide mixed evidence on the origin of New Zealand black swans. The hypothesis that reintroduced individuals mixed with remaining local individuals and that the subsequent dramatic population expansion may have been due to genetic rescue of the inbred indigenous population cannot be discarded and needs further investigation.

are carefully monitored, they can provide important insights into evolutionary and ecological processes in real time (Sax et al., 2007).
Documenting the introduction history of species beyond their traditional distribution is critical to understanding their ecological impacts on novel environments, helps predict problematic introductions, and potentially guides proposals to translocate species in the face of climate change. Uncertainty regarding provenance may even lead to controversies about whether a species is native or introduced, a designation that radically alters the management goals (e.g., extermination versus conservation) and attention afforded to species (Grosser et al., 2016).
While archeological evidence or historical records provide important insights into the introduction history of non-native species (Middleton, 2012), genetic studies allow direct hypothesis testing of various evolutionary scenarios (Pool et al., 2010) based, for example, on historical accounts of a species' introduction into a new environment. These scenarios can then be tested via analyses of extant genetic variation to reveal insights uncovered by nongenetic research techniques.
Black swans (Cygnus atratus), a common wetland species throughout much of Australia and New Zealand, represents a challenging case study for molecular and introduced species ecology, with direct consequences for wildlife management. The population status of the species as native in Australia is uncontroversial, although reports from the 17th century suggest that the species was rare in Western Australia, subsequently increasing in abundance after introductions during the 19th century from eastern Australia (Marchant & Higgins, 1990).
Contrarily, no consensus exists on the origins of the current New Zealand population. Fossil records from the early to late Holocene indicate that the species historically occurred in New Zealand, but it was believed extinct by the 19th century (Thomson, 1922;Worthy, 1998Worthy, , 2004. Several attempts of reintroduction from Australia took place 1864-1870, to control aquatic weeds in rivers, and for ornamental and hunting purposes (Thomson, 1922;Williams, 1964). Birds were released always in small groups of 1-8 units, with one big release of 42 individuals in 1867, with a total of 86 individuals reportedly released (Thomson, 1922). After reintroduction, the population grew unexpectedly rapidly, dispersed hundreds of kilometers from liberation sites (Miers & Williams, 1969), and was already common throughout the entirety of New Zealand in the first decades of the 20th century (Thomson, 1922).
Propagule pressure (i.e., the number of individuals introduced or the number of introductions) stemming from a few introduced individuals of a species with a relatively low reproductive output, such as the black swan (Braithwaite, 1981), may not have been sufficient to establish a viable population (von Holle & Simberloff, 2005;Colautti, Grigorovich, & MacIsaac, 2006). It was therefore hypothesized that the quick proliferation of the black swan in New Zealand was due to spontaneous migration from Australia during the same period as the human-driven reintroductions (Kirk, 1895). However, while individuals may move hundreds of kilometers, most remain near their natal site, and adults are highly site faithful (van Dongen et al., 2015;Payne et al., 2012;Williams, 1977). Another possibility is that New Zealand black swans were mistakenly declared extinct and that a native population was still extant at the time of the reintroduction. This population may have then subsequently admixed with the newly arrived individuals leading to the dramatic increase in population growth through genetic rescue of the inbred indigenous population (Tallmon, Luikart, & Waples, 2004). Relevant information on the origin of current New Zealand black swans most likely resides in genetic data and might have important implications for species management. Black swans are protected from widespread hunting in Australia under state legislation (e.g., the New South Wales National Parks and Wildlife Act 1974; Victorian Flora and Fauna Guarantee Act, 1988), although occasional regulated culls occur under a strict licensing system (Donavon, 2008).
In contrast, the species, being considered a pest in New Zealand, as it does in some parts of Australia, causes considerable agricultural damage to crops (in NZ since the early twentieth century; Lever, 2005). It is also declared as game under the New Zealand Wildlife Act (1953) and is commonly hunted (Nugent, 1992). The only exception is for the species on the Chatham Islands (NZ) despite reports that the population there is also nonindigenous, having arrived some time before 1922 (Thomson, 1922). Competing hypotheses, that is, whether the populations are introduced (the prevailing assumption) or whether they may at least in part be indigenous (an untested proposition), regarding the origin of current New Zealand populations have not yet been tested using genetic data.
Here, we conducted genetic typing of 12 short tandem repeat (STR) loci in samples collected from seven different geographic locations, spanning Australia and New Zealand. Using these data, we performed a variety of population genetic inferences to explore whether a recent reintroduction of black swans in New Zealand (i.e., a founder effect followed by a rapid population expansion) is unambiguously the most likely scenario. We show that alternative hypotheses could also explain the observed genetic patterns, thus questioning the current notion of New Zealand black swans being an introduced species. This study is the first large-scale contribution to the understanding of the evolutionary and reintroduction history of the iconic black swan across Australia and New Zealand using genetic data.

| MATERIALS AND METHODS
Black swans are currently distributed throughout New Zealand, and through all of Australia except northern Cape York Peninsula, with an apparently range gap in the central and western deserts (Marchant & Higgins, 1990). It inhabits fresh and low-energy marine wetlands, where it breeds with partial seasonality, mostly in pairs, with a modal clutch size of 5 -6 eggs; longevity, fidelity to partner and site, and capacity to multiclutch are poorly documented (Marchant & Higgins, 1990).
Fragment sizes were scored using an automated binning system in the Beckman Coulter 8000XL fragment analysis software, which were also confirmed visually. This binning system is well established for these loci and has previously been used elsewhere for this species (van Dongen et al., 2015;Kraaijeveld et al., 2004).

| Summary statistics
Genetic summary statistics for samples and loci were calculated with R package adegenet (Jombart, 2008) and R package hierfstat (Goudet, 2005). Observed heterozygosity was tested from departure from expected heterozygosity using a two-sided Student's t test. Pairwise F ST between populations were estimated with adegenet, while pairwise G ST (Hedrick, 2005) and pairwise Jost's D (Jost, 2008) were estimated using R package mmod (Winter, 2012 To test the correlation between individual genotypes and geography, a redundancy analysis was run in the R package vegan (Oksanen et al., 2015) and R package fossil (Vavrek, 2011), where genotypes are the dependent observations and geography is used as the independent explanatory variable. Latitude and longitude were tested together and separately, for the latter case using matrices of linear distances calculated as the pairwise difference of each coordinate. An ANOVA test was run with 10 k permutations to evaluate statistical significance of the fraction of explained variance (Legendre, Oksanen, & ter Braak, 2011). A Mantel test, implemented in R package ape (Paradis et al., 2004), was added to quantify the proportion of genetic distance that may depend on physical distances.

| Population structure
We used STRUCTURE 2.3.3 (Pritchard, Stephens, & Donnelly, 2000) to estimate the number of genetic clusters (K) in our dataset.
STRUCTURE assigns each individual to a cluster with probability q.
Ten independent analyses were performed for 1 ≤ K ≤ 10 using a Markov chain of 600 k iterations, discarding the first 100 k as burn-in.
To avoid potential biases, we assumed the admixture model without including prior knowledge of the origin of samples. Correlated allele frequencies were also assumed.

| Spatial analysis
A spatial analysis was carried out using the R package SpaceMix (Bradburd, Ralph, & Coop, 2016), which estimates admixture due to long-distance movement versus common shared ancestry. The method assumes that covariance among populations should be lower with increasing geographic distance. Higher covariance than expected on the long distance is interpreted as the result of gene flow, and lower genetic covariance than expected at short distances is interpreted as isolation. The SpaceMix algorithm provides information on how gene flow among observed genetic units (in this case our geographic localities) compares to actual geographic distances and gives an estimate of the genetic admixture among populations.
The results are visualized with a "geogenetic" map of "perceived" genetic distances based on the estimated gene flow among populations superimposed on actual geographic distances. Default setting values were used. Total run length was set to 10e7. Five parallel runs gave convergent results.

| Testing geographically driven dispersal
Bayesian estimates of migration rates among localities were obtained using BIMr (Faubet & Gaggiotti, 2008). The software calculates recent gene flow from the gametic disequilibrium that usually derives from immigrant alleles (Faubet & Gaggiotti, 2008), while linkage equilibrium is considered a lack of recent migration. To set initial search parameters for the Bayesian analysis with BIMr, twenty pilot runs of 50 k iterations are performed, starting from the default parameter set (Faubet & Gaggiotti, 2008). Thus, a burn-in of 100 k was performed and the final parameter search was conducted for 5 million iterations.
Ten parallel replicates were performed and output files were analyzed with Tracer_v1.6 (Rambaut & Drummond, 2007) to evaluate the posteriors. Internal convergence of single runs and convergence of parallel runs were tested with R package coda (Plummer et al., 2006).
After conversion checking, five convergent runs estimated the best posterior parameters and quantiles for migration rates and inbreeding coefficients.

| Effective population size and time to the most recent common ancestor
We used a coalescent-based estimator of population size changes over generational time (MsVar; Beaumont, 1999). We used version 1.3, which directly estimates mutation rate (μ), population effective sizes (current and ancestral, or N 0 and N 1 , respectively), the time since Convergence was reached for all independent runs using a running time of 12e10 iterations (SFile 1). Following Beaumont (1999), Bayes factor was used to evaluate the probability of declining versus expanding population size changes. Modes and confidence intervals of posterior parameters were estimated with R package modeest (Poncet, 2012).

| Bootstrap analysis of the dataset
We selected a subset of 70 random samples from our genotypes to be analyzed with MsVar, using the same prior parameters and running time as in our populations. The number of samples per subdataset was the same as in observed populations (Table S1), with each sample repeated 10 times. This analysis allowed us to build general distributions of posterior population parameters for the whole metapopulation of black swans and to use them to test the likelihood of demographic parameters for our populations. We also tested the effect of number of samples on the estimates.
The logical flow behind our statistical approach has been summarized in a paragraph provided in Supplementary Material. The main R scripts used in this study are provided in Dryad repository.

| Genetic variation and population differentiation
A similar level of internal heterozygosity occurred across populations and loci (Tables S2 and S3). No substructure was observed for pooled locations (see Materials and Methods and Table S1) using the find.cluster function in R package adegenet (Jombart, 2008 Fig. S3).
Redundancy analyses of genotypes versus geographic coordinates revealed that geography explains a small but significant fraction of genetic distances (proportion of inertia explained with geography = 0.043, e ≪ 0.001). When geography is summarized into a single factor (i.e., sampling locality), 9% of the variance in the distribution of allele frequencies is explained (p < .001). A Mantel test indicated a lack of correlation between the matrix of geographic distances and the genetic distances (p > .12). However, latitude taken alone seems to explain genetic distances (p < .00), while longitude is nonsignificantly correlated with genetic (p > .189), as shown in plots of genetic distances against physical distances, using distances in kilometers for the two-coordinate system and absolute distance in degrees for latitude and longitude separately ( Fig. S2D-F).

| Genetic structure
Genetic clustering analysis with STRUCTURE detected four genetic clusters among the samples (Fig. S1). Location differentiation detected with F ST (>0.01) is sufficient to obtain a reliable estimate of the genetic structure (Coleman, 2014). The most likely number of clusters in the dataset was estimated by calculating ΔK following the procedures of Evanno, Regnaut, and Goudet (2005), which estimates the most likely number of clusters based on the rate of change in log-likelihood probabilities for each K. ΔK was greatest for four clusters, although ΔK for two clusters was also relatively high (Fig. S1). Following Pritchard et al. (2000), we assumed four clusters as this created a more biologically likely clustering than when two clusters were assumed. In the latter F I G U R E 2 Distruct plot of STRUCTURE posterior cluster assignment probabilities among black swans. Each vertical line represents one individual, and color indicates probability of belonging to each of the four genetic clusters case, samples from WA and QL would be clustered together, which is geographically and ecologically very unlikely. Individuals generally displayed admixed origins; however, QL, WA, and BA appear to be the most internally homogenous and thus displayed greater differentiation compared with other localities (Figure 2). The remaining four localities show similar within-group composition, although WE and WL are more admixed than the two island localities (TS and NZ) where most individuals belong to one cluster (Figure 2). This common ancestry could be a wide southeast component, with exception to be made for Ballarat (see discussion below in this regard). The STRUCTURE analysis is therefore in good agreement with both the F ST and PCA results, which show less differentiation for southeastern Australian samples.

| Admixture and spatial origin of the populations
SpaceMix estimated a very weak degree of genetic admixture among populations (see "geogenetic" map in Figure 3). The map showed that the spatial origin of distant samples (NZ, QL, WA, WE, and WL) is a good approximation of genetic covariance, indicating that greater physical distances and genetic distance can be explained by genetic ancestry. On the other hand, two of the closest localities (BA and TS) are more scattered on the genetic than on the geographic map. This pattern indicates that BA is more differentiated than would be expected under a simple isolation-by-distance model, while TS seems to be more closely related to mainland samples than expected given the physical distances. Therefore, SpaceMix analysis complements the results obtained with our previous analyses, and clarifies why the Mantel test failed to detect a linear correlation between geographic distance and genetic distance matrices.

| Estimates of recent asymmetric migration rate
Bayesian estimates of multilocus migration rates confirmed a lack of recent ongoing gene flow among black swans. As for the STRUCTURE analysis, the F ST values observed among the seven populations fall between 0.01 and 0.05 (Table S5), a range for which the estimator gives acceptable results, although posterior estimates improve with increasing genetic differentiation (Faubet & Gaggiotti, 2008). Five of ten parallel runs were convergent. The results indicate migration from one locality to another to be < 1% (Table 1). These estimates give a quantitative perspective on the degree of current isolation among our samples. Closer localities appeared to be as isolated as more distant ones, which further supports the outcome of the SpaceMix analysis.
Interestingly, the remaining five nonconvergent runs differed from those that converged only for the NZ results. Each of these runs indicated some degree of gene flow into NZ from all other locations, although the total amount of gene flow substantially changed from one run to the other.  (Fig. S4). Given this, a model of population growth has virtually zero likelihood compared to a model of population decline (Beaumont, 1999). As a substructure in the sample may produce a false signal of population decline, Chikhi et al. (2010) proposed that random subsamples of the metapopulation can be analyzed to verify whether the signature is also detected at the metapopulation level to confirm decline in subpopulations. Thus, a bootstrap test was run on 70 sub-datasets, which showed that the signature of strong decline is also detected in the metapopulation (black diamonds in Fig. S4).

| Population size changes over time
The absolute values of current effective population sizes range from ~50 to ~150 individuals ( Table 2). The ratio between current and ancestral effective population sizes (r;  (Thuillet et al., 2002).
F I G U R E 3 SpaceMix analysis plot. Population names are abbreviated as in Table 1. The black names represent the geographic location based on coordinates, while color-coded names represent the "geogenetic" distances (i.e., the distances perceived by the populations depending on their genetic covariances)

| Time to population decline and time of the most recent common ancestor
In MsVar 1.3, the generation time is specified in the input file and the posterior estimates of time since population size change are converted into number of years. Black swans can reproduce twice per year once sexual maturity is reached (18-36 months); their age of first successful reproduction can vary between 2 and 5 years (Coleman, 2014 T is expressed here in number of years, equal to the number of coalescent units outputted by MsVar multiplied by the effective current population size and generational time (Table 3), while r is modal and quantile values of the distribution of the ratio N 0 to N 1.

| DISCUSSION
While the ecology and population biology of black swans have enjoyed substantial research effort, little is known regarding the species' molecular ecology and genetics (see van Dongen et al., 2015;Kraaijeveld et al., 2004). We used a combination of inferential evolutionary tools to document patterns of ongoing gene flow among populations across Australia and New Zealand and examine the long-term population demography of the species. Our results highlight three important findings: (i) a lack of substantial dispersal between populations, (ii) a strong demographic decline across all populations, approximately 1 -2 kya, and (iii) a more ancient origin of New Zealand and Tasmania samples compared to samples from most Australian localities. Our results may question the prevailing notion that the black swan is a highly mobile species that originated in Australia, which owes its high current abundance in New Zealand to a coincidence of human-driven releases and spontaneous, novel, and substantial immigration from Australia.
Our multivariate analyses, coupled with multilocus inference of recent migration rates and a landscape analysis of covariation between genetic ancestry and geographic distances, indicate that dispersal rates of this species are low, despite the capacity of the species to fly long distances (Williams, 1977). In support of this, van Dongen Werribee could act as a sink locality for black swans, such scenario is not mirrored in our estimates of migration rates. On the other hand, Ballarat is characterized for being the highest location of our sampling and in general one of the highest cities in Australia. This may suggest a role of altitude in shaping black swans' genetic diversity, a topic that would deserve further investigation.
Overall, the results suggest that this species does not disperse long distances when local environmental conditions are satisfactory, although long-distance dispersal may be more likely when local conditions deteriorate or are less stable (Williams, 1977). The fact that latitude seems to be more correlated with genetics than longitude was expected from previous studies on bird dispersal, as latitude is more correlated with environmental factors, that is, temperatures and precipitations, than longitude (Masello et al., 2015). Such a strategy is typical of many Australian birds which rely on waterbodies, and is also observed in environments where rainfall can be temporally and spatially unpredictable (e.g., Roshier, Asmus, & Klaassen, 2008;Roshier, Robertson, & Kingsford, 2002). Although we could not formally test an hypothesis relating climatic variables to genetic differentiation in black swans, our results points to an indirect correlation of genetic differentiation with physical distances. Indeed, even if the longitudinal dispersal (1139-3334) of our population is much higher than the latitudinal range (501-1009), latitude turns out to be a better predictor than longitude for genetic differentiation. Thus, climatic differentiation along a south to north gradient could influence genetic differentiation more than physical distances. We thus suggest that this hypothesis could be tested also in other Australian aquatic birds in relation to dispersal and genetic differentiation. Other swan species also tend to remain at the same site year round and not disperse from breeding sites between breeding events (e.g., Włodarczyk et al., 2013). We found pronounced differentiation between Western Australian swans and the other populations despite reports that the population may have at least partly arisen from introductions from eastern states (Marchant & Higgins, 1990).
Our results also indicate that all populations examined have declined over the last 1-2 kya. We found no signature of an increase in population size in New Zealand or Western Australia that would be expected from the reported population increases. This may be due to a lack of statistical power in our analyses to detect such a recent changes in population size. It is also to note that the expansion has occurred over few generations. Alternatively, the decrease in population size over the last 2,000 years may have overwhelmed the signature of the more modest contemporary increase. In the absence of experimental studies, drivers of population change are correlative and speculative, but we propose climate changes and hunting pressure as a possible explanation of detected population decline. Studies on the paleoclimatic changes in Australia during the Holocene (Neolithic) indicate a period of peak moisture 5-8 kya during which wetland water levels were higher than present (Quigley et al., 2010). Aridification commenced approximately 5 kya, where levels of salt lakes in mainland Australia fluctuated between moderate to low and then dropped drastically over the last 1 ky (Bowler, 1981;Quigley et al., 2010). This presumably resulted in severe fragmentation of swan habitat, followed by isolation of habitats and local populations. Indeed, the incomplete population structure and the lack of recent gene flow we report indicate that the current populations might have resulted from recent geographic isolation followed by population divergence. The effects of habitat change on population structure may be continuing in the present day, with widespread wetland loss or conversion (Anonymous, 1988;Goodrick, 1970).
In addition to the onset of aridification in Australia, the last 5 kya in Australia has also been characterized by a rapid acceleration in human population growth (Johnson & Brook, 2011). Indigenous Australians hunted black swans for food (e.g., Clarke, 2003), and the rapid expansion of the European population may have accelerated the species' population decline via increased hunting pressure. The apparent coincidence of these two processes, coupled with our coalescent inferences of population decline in black swans, would suggest a reasonable ecological explanation for the decline of swans. Similarly, the decline in black swans in New Zealand may be at least partially attributed to the arrival of human populations approximately 800 years ago (Wilmshurst, Anderson, Higham, & Worthy, 2008). The arrival of the Maori leads to rapid landscape transformations (McWethy et al., 2010) and was coupled with a high extinction rate of indigenous avifauna due to hunting pressure (Duncan & Blackburn, 2004), including the hunting of black swans (Horn, 1983).
From a population genetics perspective, estimates of population size changes over time may often be biased by the confounding effect of population structure (Chikhi et al., 2010;Mazet et al., 2016), as most methods designed to estimate this parameter, including MsVar (Beaumont, 1999), assume panmictic populations. However, no population substructure could be detected among our population units. In addition, a bias toward detecting population bottlenecks with MsVar has been reported by Chikhi et al. (2010), but a simulation study by Girod et al. (2011) showed that the method can reliably detect reductions in population sizes when the bottleneck was strong, meaning that the population size reduced by a factor of at least 10e3, which seems to be the present case. Excess immigration into the population could also misleadingly lead to a signal of recent bottleneck (Excoffier & Heckel, 2006). However, all our analyses indicate little to no gene flow among populations and the population reduction detected did not begin recently, as supported by simulation studies (log τa < 1; Girod et al., 2011). This is further corroborated by our simulation approach based on 70 random datasets, following Chikhi et al. (2010), rendering our inferences of metapopulation decline highly reliable.  Bouzat et al., 2009;Hogg, Forbes, Steele, & Luikart, 2006;Madsen, Shine, Olsson, & Wittzell, 1999;Tallmon et al., 2004). The population's sudden increase may also have been related to other factors such as the simultaneous, substantial, and novel immigration of swans from Australia (Kirk, 1895) or changes in ecological conditions (e.g., changes in land use) that encouraged population growth (Crookes & Soulé, 2001 Overall, the results presented here may contradict the general consensus that the New Zealand black swan is an introduced (pest and game species) from Australia. Instead, the New Zealand population appears to be at least partly derived from a natural, and not a feral, population; thus, it may constitute an evolutionary significant unit (ESU; Moritz, 1994). A test of mtDNA reciprocal monophyly of NZ, TS, WE, and WL samples could provide more insightful evidence on the definition of ESU (Moritz, 1994). The designation of species as invasive and non-native is complex and the nuances of population histories need to be considered (Crees & Turvey 2015). In particular, approaches to actively managing "native" invaders require addressing social and political considerations (Carey, Sanderson, Barnas, & Olden, 2012). Our results suggest that, in the nuanced framework of Crees and Turvey (2015), the black swan in New Zealand is "xenonative" (rather than "returned" or "restored" native), that is, native but restocked from external populations. The management strategies and legislation regarding the conservation and hunting status of the swan in New Zealand may accordingly need to be reviewed. While swans are currently numerous in New Zealand, extirpation would lead to the loss of genetic variation attributable to more ancient endemic swans, and as such should perhaps be avoided, or at least a more nuanced decision made regarding desirable outcomes for the species' population.

| FUTURE DIRECTIONS
Our study is the first in-depth assessment of the population demographic history of the black swan and provides the first convincing evidence that the current black swan population in New Zealand may be partly derived from indigenous populations. However, the complex evolutionary scenarios that we propose in the present work will require additional analyses to provide further clarification. The availability of a precise historical record on the reintroduction of black swans in New Zealand allows the collection of samples from exact localities where the individuals were captured (Australia) and released (New Zealand). Such a sample could then be typed for higher-resolution genomic loci than microsatellites, such as radSeq data, or the whole mitochondrial genome, which was published in 2010 for C. atratus (Jiang et al., 2010) and might reveal unambiguous phylogeography of Australian and New Zealand lineages. Such dataset will likely have enough power to run in-depth hypothesis tests to verify our indications of a New Zealand or Tasmanian origin of C. atratus, including the possibility to detect past sporadic migratory event between Australia, New Zealand, and Tasmania and their intensity.

AUTHOR CONTRIBUTIONS
P-JG, RAM, RWR, and MAW conceptualized the project and study design. P-JG, RAM, and WFDvD conducted field work. P-JG conducted the laboratory work. VM, WFDvD, and MC analyzed the data. P-JG, RAM, RWR, MAW, VM, and WFDvD interpreted the findings, WFDvD, VM, and MAW wrote the manuscript. All authors read and contributed to the manuscript.

DATA ACCESSIBILITY
The microsatellite data used for this study and an R script for the data analysis are available from Dryad Digital Repository: https://doi. org/10.5061/dryad.ct8j4. Any specific question about data analysis will be answered upon request to the corresponding author.