Evaluating a potential source of founders for ex situ conservation efforts : genetic differentiation between disjunct populations of the endangered red siskin Spinus cucullatus

DNA Amplification and sequencing – For cytochrome B, our PCR reactions included 1x PCR buffer, 1.5mM MgCl2, 0.2mM of each dNTP, 0.8uM of each primer, 1% DMSO, 1 unit Biolase Taq, and ~25ng template DNA in a total volume of 25ul. The PCR thermal profile began with preheating at 95C for 3 min, followed by 35 cycles of denaturing at 95C for 30 seconds, annealing at 50C for 30 sec and extending at 72C for 1 min. The final extension step lasted 10 min.


INTRODUCTION
Captive breeding for reintroduction is a conservation technique that can be important in preventing extinctions (Butchart et al. 2006, Hoffmann et al. 2010).One of the most basic factors that may influence the success of a captive breeding/reintroduction effort is obtaining appropriate founders (IUCN/ SSC 2013).Founders should ideally be representative of populations from the region to be restored, and come from wild stock which has not adapted behaviorally, physiologically, or genetically to captivity (Ballou et al. 2010).
Unfortunately, captive breeding programs are of ten established haphazardly, without foreknowledge of future conservation needs, or under logistical constraints that prevent considering founder origins.When species are long-lived, overcoming early mistakes can be difficult (Russello & Amato 2004, El Alqamy et al. 2012), and, for short-lived species, additional founders may need to be regularly introduced into the captive population to prevent in breeding (Ballou et al. 2010;e.g. Hedrick et al. 2012).Therefore, performing genetic evaluations of reintroduction programs, both before and after conservation actions are taken, can result in valuable management insights (e.g.Tollington et al. 2013).However, when species are rare, obtaining sufficient information for such evaluation can be difficult.
The red siskin Spinus cucullatus, listed by the IUCN as 'Endangered,' is a bird for which captive breeding could be an important safeguard against extinction (BirdLife International 2017), but for which obtaining appropriate founders is a challenge.Intensive and ongoing trapping has decimated this small seedeater throughout its originally known range, which was mainly in northern Venezuela, stretching from border areas of Colombia to Trinidad (Coats & Phelps 1985, Castro & Asuaje 2013) (Fig. 1).It has been protected by regulation in Venezuela since the 1940s (Coats & Phelps 1985), has been listed in the Convention on the International Trade in Endangered Species (CITES) Appendix I since 1975 (CITES 2013), and is protected by more recent national legislation in Venezuela, Guyana, the USA, and other countries worldwide (Venezuela 1996, USFWS 2017).Critically Endangered in Venezuela, it is also threatened by habitat loss, as the tropical dry forests that are an important part of its habitat are also endangered (Rodríguez et al. 2010, Rodríguez-Clark et al. 2015).However, if the trapping threat were mitigated, sufficient habitat remains to support reintroduced populations in the future (Coats & Phelps 1985, J. Miranda & A. Sánchez-Mercado unpublished data).
The recently founded Red Siskin Initiative is an international consortium that aims to promote this species' recovery in the wild (Red Siskin Initiative 2017).Organizations assisting with the captive breeding program include the Smithsonian Institution, which recently established an ex situ colony of red siskins for research and education; ZooMiami, which will establish a colony soon; the National Finch and Softbill Society; the Venezuelan NGO Provita; and Parque Zoológico y Botanico Bararida, a Venezuelan zoo interested in maintaining captive populations to receive confiscated individuals from the illegal trade, to educate stakeholders, and to breed birds for eventual reintroduction.
Until recently, only 2 possible sources of founders for red siskin conservation breeding efforts were known: captive individuals presumed to be of Venezuelan origin, and wild-caught individuals from Venezuela.Obtaining individuals from the wild is a challenge because the species is exceedingly rare, nomadic, and only sporadically present in the few locations where it is currently known (J.Miranda & D. Ascanio unpublished data).Moreover, deliberately bringing wild individuals into captivity further reduces surviving wild populations, and may also alert trappers to their location.Illegally traded red siskins are occasionally confiscated in Venezuela, but whether they are of wild or captive origin generally cannot be determined with certainty.Captive birds, on the other hand, are readily available, as thousands of individuals are held as pets worldwide.However, with captive animals there may be problems relating to domestication, inbreeding, legal status, or hybrid ancestry.Historically, red siskins have  Recently, a third potential source of founders was discovered: a population in Guyana, nearly 1000 km distant from the previously known Venezuelan range (Robbins et al. 2003; our Fig.1).Surveys suggest that this population is healthy and stable, comprising hundreds to thousands of individuals, and protections are in place to safeguard against illegal trade (SRCS 2017, https://www.facebook.com/southrupununiconservationsociety/). The Guyana population thus represents a potentially attractive source of founders for captive breeding and reintroduction in Venezuela.
However, the disjunct nature of this population raises questions about its potential level of differentiation from other populations.Neotropical bird populations frequently have stronger phylogeographic structure than their north temperate counterparts, presumably due to longer residence times (e.g.Smith et al. 2014).In addition, Guyanese red siskins may be physiologically adapted to a different habitat; the Guyanese population inhabits a hot, low-elevation savanna/forest ecotone, while Venezuelan siskins occur at higher elevations in more mesic habitats (Coats & Phelps 1985, Robbins et al. 2003).Thus, the Guyanese population may be sufficiently differentiated to make it undesirable as a source of founders for Venezuela.On the other hand, the entire South American siskin radiation of 10 species appears to have been recent and rapid (Beckman & Witt 2015), and a vicariant origin of the Guyana population based on savanna expansion in the region would imply isolation of just 8000−10 000 years (Van der Hammen 1983).
The Guyana population could also be the result of recent long-distance dispersal or an anthropogenic introduction (Robbins et al. 2003).In either of these scenarios, we would expect the gene pool of the Guyana population to be a subset of the Venezuela gene pool, with the likelihood of reduced diversity due to founder effects.The red siskin is a semi nomadic, flocking species which may therefore exhibit long-distance dispersal.Red siskin populations in Cuba and Puerto Rico may be derived from escaped cage birds (Raffaele 1983, Lever 1987, Collar 1992), and many other feral populations of small finches exist in the Guianas and the Caribbean (Bond 1971).Yet, while Guyana has been a source of bird trafficking for centuries (Hanks 2005), at the time of discovery of the Guyana population, traffic in red siskins was unknown in that country (Robbins et al. 2003).
In the present paper, our aim was therefore to evaluate both nuclear and mitochondrial (mtDNA) genetic differentiation between the Guyanese and Venezuelan populations, in order to explore the possibility of using Guyanese individuals as founders for a captive breeding program to restore populations in Venezuela.Here we report the results of sequence comparisons of 2 mtDNA genes (cytochrome B and control region) and nuclear genetic variation at 312 loci surveyed by the amplified fragment length polymorphism (AFLP) technique.A finding of little to no genetic differentiation would support the use of individuals from Guyana in ex situ conservation efforts aimed at restoring red siskins in Venezuela, while a finding of significant differentiation would caution us in such an endeavor.

Sample collection and DNA extraction
Samples of this species are extremely difficult to obtain both in the wild and from captive flocks.In the wild, birds are nomadic, high-flying, sparsely distributed, and Endangered, and permits for their capture and sampling are a major challenge, often requiring years of effort.Many owners of captive birds do not have proper paperwork and/or are re luctant to allow sampling of valuable, delicate birds.Thus, the samples we were able to obtain were not ideal.However, given the importance of the conservation question at hand, we preferred to use available samples and carefully consider potential sources of bias when drawing what we believe are conservative inferences.
The 5 samples available from Guyana (GU; Table 1) were from wild birds in adult plumage sampled at a single location on 12 April 2000 following their unexpected discovery during an ornithological survey of the Rupununi Savanna, conducted with the permission of the Guyana Environmental Protection Agency and Ministry of Amerinidian Affairs (Robbins et al. 2003;our Fig. 1).Efforts by M.J.B. and M. Robbins, after this discovery, were instrumental in establishing legal protection for this species in Guyana.Samples were frozen in the field in liquid nitrogen and maintained at −130°C or below.
Most or all red siskins in captivity today worldwide are thought to derive from Venezuelan stock, al -though records of their origin are generally unavailable.The 13 individuals used to represent Venezuela ('VE') in the present study came from 5 captive flocks in 3 US states (Table 1).They were donated for this study following death from natural causes, and had no morphological traits suggesting hybrid ancestry.These 'VE' individuals consisted of 2 pairs of known siblings and 9 individuals with no known first-order relationships.However, given the small size of US avicultural flocks and regular transfers among them, some level of relatedness and inbreeding is plausible among all 13 (P.Hansen pers.comm.).Captive specimens were stored at −20°C after death, shipped on dry ice, and subsequently maintained at −80°C.Genetic samples and voucher specimens will be deposited at the US National Museum of Natural History (USNM) of the Smithsonian Institution, in Washington, DC.
Genomic DNA was extracted from all GU samples and 'VE' samples B20014 through B20021 using standard phenol− chloroform extraction.(Sambrook et al. 1989).DNA from the remaining 'VE' samples was extracted on an automated Autogenprep 965 extractor (Autogen) following the manufacturer's instructions using a standard mouse tissue protocol.DNA concentration and purity were assessed using a Nan-oDrop ND-1000 spectrophotometer (ThermoFisher Scientific).

Mitochondrial DNA amplification, sequencing, and editing
We amplified 2 mitochondrial ge nes: cytochrome B, and a portion of the control region.For cytochrome B, we used the primers L14764 (5'-TGR TAC AAA AAA ATA GGM CCM GAA GG-3'; Sorenson et al. 1999) and H16060 (5'-TTT GGY TTA CAA GAC CAA TG-3'; Robbins et al. 2005) to amplify the entire coding sequence and short flanking regions.For the control region, we designed species-specific primers targeting an initial segment of ~670 bp: RSCRL000 (5'-CTC TCT CCG AGA TCT ATG GCC TGA A-3') and RSCRH690 (5'-CAC TTG AAG GGC TTA TTG AAG AGA C-3').All amplicons were sequenced on both strands with additional internal primers, and reads were assembled with Sequencher 5.0 (Gene Codes) to arrive at consensus sequen ces for each individual.Full details of amplification and se quencing protocols are given in the Supplement at www. int-res.com/ articles/ suppl/ n036 p183 _ supp.pdf.No insertions or deletions were detected in either gene, and no stopcodons in the case of cytochrome B, suggesting that our sequences were of mitochondrial origin and not nuclear pseudogenes (Sorenson & Quinn 1998).Genes were concatenated for all analyses.

AFLP scoring
To compare GU and 'VE' individuals across the nuclear genome, we developed a set of AFLP markers (Vos et al. 1995, Bensch & Akesson 2005, Meudt & Clarke 2007).We used the protocol of Kingston & Rosel (2004) with some modifications (see 'Detailed methods' in the Supplement) to screen all 18 individuals for variation with 14 selective primer pair combinations (Table S1 in the Supplement).We used an ABI Prism 3100 genetic analyzer to detect fragment sizes, multiplexing 2 selective PCR products labeled with different dyes in each run.Electropherograms were scored using GeneMapper 4.0, following Kingston & Rosel (2004).Polymorphic peaks were scored as dominant, biallelic markers (Vos et al. 1995).We used strict scoring criteria to minimize errors, only scoring peaks larger than the second smallest size standard (89 bp) and smaller than the second ) scored all fragments in all individuals separately and re moved any loci with discrepancies.In total, we developed 312 loci that could be reliably scored.

Data analyses
Data from all individuals (5 GU, 13 'VE') were used for the analyses described below.mtDNA To examine variation in mitochondrial DNA sequences, we constructed a median-joining network using Network (Fluxus Technology Ltd. 2009).We also calculated total haplotype diversity within each population (Ha t ), and nucleotide diversity (π) as well as divergence with F ST and Φ ST as implemented in Arlequin (Excoffier et al. 2005).

AFLP
Population allele frequencies for all AFLP loci were estimated from observed fragment frequencies in the GU and 'VE' samples using a Bayesian approach implemented in AFLP-SURV (Zhivotovsky 1999).We applied a non-uniform prior of allele frequencies computed by combining sample size and the number of individuals without fragment presence to take into account small sample sizes.We estimated allele frequencies assuming Hardy-Weinberg equilibrium (HWE; F IS = 0), as well as assuming average F IS values ranging from 0.0625 (individuals related on average at the level of second cousins) to 0.75 (the equi valent of 4 generations of full-sib mating).In order to assess the effect of possible size homoplasy on our estimates of genetic divergence with these markers, we also calculated the average fragment size and the Pearson correlation coefficient (r) between fragment size and frequency, along with its significance (Vekemans et al. 2002).Finally, to better separate and understand possible sources of bias, we calculated relatedness (r ab ) of each individual with respect to individuals from the country of origin of that individual, following Lynch & Milligan (1994), and compared them with levels of relatedness between known sib pairs, to test assumptions about probable levels of F IS .
The percentage of polymorphic AFLP loci (at 5% or above; PL), Nei's gene diversity (H j , equivalent to expected heterozygosity), and F ST between the 2 populations were computed with AFLP-SURV (Vekemans 2002), using the range of allele frequency estimates described above.F ST values were tested for significance by comparing observed values with the distribution of values in 10 000 random permutations of individuals among groups, calculated on the basis of expected heterozygosity of dominant marker loci (Lynch & Milligan 1994, Vekemans 2002).We also searched for significant linkage disequilibrium (LD) among all locus pairs using an algorithm for dominant markers (Li et al. 2007).LD analysis can also reveal aspects of population structure not evident in individual-locus analyses, since extensive non-random associations of allele frequencies across many loci can indicate recent founder events and/or bottlenecks that would be expected if, for example, the Guyana population had a recent origin from few captive individuals.
Population structure in the AFLP data was further explored using a variety of methods.To visualize differentiation, we used NTSYSpc version 2.2 (Rolf 2008) to first create a matrix of band-sharing between individuals as measured by the Jaccard similarity value, a metric appropriate for dominant loci be cause it makes no assumption of homology among band-absent genotypes (Ajmone-Marsan et al. 2002).We then represented the relationships revealed by these values using an ordination technique, i.e. nonmetric multi-dimensional scaling (NMDS; Rolf 2008).
Private alleles and diagnostic loci are expected to accumulate in isolated populations over time, with in creasing numbers indicating increasing divergence (Schönswetter et al. 2004).Private alleles are those for which fragment presence is observed in just 1 population, and diagnostic loci are those which distinguish all individuals of a population or group from all individuals of other populations or groups.We estimated the number of private alleles (N p ) and diagnostic loci (N d ) from our AFLP data, and calculated the probability of the observed numbers occurring by chance using randomization (Manly 1997).We reconfigured our dataset 1000 times (considering only 1 randomly-selected sib per known sibset), reassigning individuals each time randomly to a group (GU or 'VE').For each reconfiguration, we tallied the number of private alleles and diagnostic loci between those 2 groups.We then calculated a p-value as the proportion of null model iterations where N p and N d exceeded or were equal to the observed values.Calculations were performed with Excel (Microsoft).
In addition to classical analyses of predefined groups, we used 2 clustering methods to investigate population structure in our samples.We first applied model-based Bayesian clustering analyses using algorithms appropriate for dominant loci as implemented in STRUCTURE 2.3.4 (Pritchard et al. 2000, Bonin et al. 2007, Falush et al. 2007).These analyses assume HWE, an assumption likely to be violated by the known and potential familial relationships present in our samples.However, although family relationships and inbreeding can lead to an overestimation of the number of distinct population clusters K, there appears to be little effect on the correct assignment of individuals to populations for a fixed K (Falush et al. 2003, Pritchard et al. 2010).We chose not to include prior information about sample origin in STRUCTURE models, because the geographic origin of 'VE' samples is presumed rather than known.We used a standard admixture model to allow for the possibility of mixed ancestry, assumed that allele frequencies in each population were correlated, and conducted unsupervised runs with K from 1 to 6 groups.Although the correlated frequencies model can overestimate K in the presence of family relationships, it is more appropriate for populations that may share ancestry (Pritchard et al. 2010).We chose to run more chains for shorter periods, in order to examine variation among runs (Evanno et al. 2005), and used equal burn-in and data collection periods of 10 000 iterations each, for 20 independent runs per model.For each run, we recorded the estimated posterior probability of the data given the assumed model, and used STRUCTURE Harvester (Earl & vonHoldt 2012) to calculate ΔK, an ad hoc statistic based on the secondorder rate of change of the likelihood function with respect to K (Evanno et al. 2005).We took the model corresponding to the modal value of the distribution of ΔK as indicating the likely number of populations, and then plotted each individual's estimated membership coefficients in those populations (Q) from a representative run (as variation among runs was minimal).To consider the case of K = 1 (i.e.lack of genetic structuring), we examined variation in α, a model parameter indicating the extent of admixture: variation among iterations beyond a range of 0.2 units or more in a single run indicates a lack of true structure.We also examined individual assignment and Q values for models corresponding to peak ΔK, because roughly equal numbers of individuals assigned to each putative popu lation and a majority of admixed individuals also indicate a lack of true structure (Pritchard et al. 2010).Finally, we also considered which value of K had the highest ln Pr(X |K ), which is recommended as an additional indicator of the true value of K (Janes et al. 2017).
We also used a model-free iterative reallocation method, FLOCK 3.1 (Duchesne & Turgeon 2012) to estimate the number of populations, K.This method is robust to population inbreeding and non-zero relatedness among sampled individuals because it creates clusters based on maximizing multilocus genetic similarity rather than minimizing deviations from HWE and LD.In this method, samples are initially partitioned randomly into K clusters (K ≥ 2), allele frequencies are estimated for each of the K clusters, and each individual is then reallocated to the cluster that maximizes its likelihood score.Twenty repeated reallocations are performed within each run, and 50 runs are carried out for each K. Strong consistency among runs, resulting in 'plateaus' of identical mean log likelihood difference (LLOD) scores, is used to indicate the most likely number of clusters (Du chesne & Turgeon 2012).Although it is not run ex plicitly with K = 1, FLOCK does test for K = 1.In short, K = 1 is the default hypothesis, and is retained if no plateau of length ≥6 is found for any K ≥ 2.
Once reference populations have been correctly identified, allocation programs take advantage of this information and so are generally less prone to misallocations than are cluster programs.Thus, we also performed reallocation procedures using the method and software designed for AFLP data by Duchesne & Bernatchez (2002;AFLPOP) to reallocate individuals to populations ('VE', GU) based on the allele frequencies.We used the default settings (fixed correction value for 0 frequencies = 0.001, minimal LLOD to allocate specimens = 0, number of artificial genotypes to compute p-values = 500).AFLPOP calculates the LLOD score for each genotype (the difference between the log likelihood of the most likely reference for the genotype and that of its second most likely reference) and the mean LLOD (MLLOD) over all genotypes.Higher differentiation between references will tend to produce higher MLLOD scores.Marker loci with minor allele frequency of < 5% were considered monomorphic and uninformative for the purpose of our FLOCK and AFLPOP analyses.We therefore discarded loci with either 17 or 18 'band present' phenotypes over all 18 genotypes.Because FLOCK and AFLPOP do not accept loci with missing scores, those were also removed.Thus 78 loci were retained among the 312 loci originally developed and scored.
Although the 5 GU samples were from wild adults, given the fact that they all came from the same location, it is possible they were related.Similarly, our 13 'VE' samples could have additional family relationships unknown to their breeders.With small sample sizes, close relatedness among individuals could produce distorted estimates of allele frequencies, mimicking population structure.Therefore, we considered the following hypotheses to explain the presence of GU and 'VE' genetic clusters: (1) GU and 'VE' belonged to 2 distinct populations, (2) GU and 'VE' belonged to the same population, but the GU genotypes were strongly inbred (related at the level of full sibs) so that enough differentiation was generated to be picked up by various algorithms, (3) GU and 'VE' belonged to the same population, but the 'VE' genotypes were the inbred ones (related at the level of half sibs).
We designed a specific procedure in order to test hypotheses that moderate or very high levels of relatedness among the 'VE' or GU genotypes might have been sufficient to explain the high MLLOD score obtained from running the reallocation procedure of the AFLPOP program.Essentially, we kept the empirical (real) set of either 'VE' or GU genotypes and generated 100 simulated sets of the other genotypes, using allelic frequencies based on the assumption that all 18 actual specimens ('VE' + GU) originated from the same population.To test hypothesis 2, we simulated sets of 5 genetic full sibs.To test hypo thesis 3, we simulated sets of 13 half sibs.Each simulated set stood in place of the empirical GU or 'VE' genotypes, respectively, depending on the simu lation.We thus ran the reallocation procedure of AFLPOP with each set of the 5 simulated full sibs and 13 'VE' genotypes, or the 13 simulated half sibs and the 5 GU genotypes.For each of the 100 simulations, the MLLOD score of the reallocation result was calculated (see Fig. S1a,b in the Supplement).To obtain p-values, we located the MLLOD score from the reallocation of empirical GU and 'VE' genotypes within each distribution of the 100 MLLOD scores from both simulation procedures.

mtDNA variation
All 18 individuals examined had cytochrome B (MT-CYB) sequences consistent with a previously published sequence for this species.No previously published sequence for the control region (CR) was found.Out of a total of 1813 bp sequenced across the 2 genes (1143 bp in MT-CYB and 670 bp in CR), there were 8 variable sites, resulting in 3 haplotypes for MT-CYB and 4 haplotypes for CR, with a total of 4 haplotypes when both genes were considered together (Table 2).One haplotype was present in all 5 GU birds and absent in 'VE' birds (Table 2, Fig. 2).That haplotype had 2 sites in MT-CYB that distinguished all GU from all 'VE' individuals.The 3 'VE' haplotypes observed differed from the GU haplotype by 3 to 6 substitutions.Total haplotype diversity was high, at 0.75, and was due entirely to diversity within 'VE' birds and differences between GU and 'VE' birds, while nucleotide diversity was low, at 0.008, again due entirely to variation within 'VE' birds (Table 3a).

AFLP profiles and relatedness
For the 312 AFLP loci scored, individuals had on average 222.8 bands present.In total, 170 loci were variable (54.5% of all loci).Average fragment size 189 Origin ID Variable site number MT-CYB MT-CR 1 2 3 4 5 6 7 8 Complete ID numbers available upon request to safeguard the location of the GU population was 279.5 bp (SD 115.5).The correlation between fragment size and frequency was negative, and was significant when calculated with respect to all individuals (r = −0.16,p = 0.0047); for individuals within each country of origin, it was smaller and not significant ('VE': r = −0.09,p = 0.18; GU: r = −0.08,p = 0.17).
Estimates of LD were low and not widespread; each locus had significant disequilibrium with, on average, just 8.3 other loci, or 3% of pairs, with the median number just 1 locus.Estimates of relatedness from AFLP genotypes were consistent with, but slightly upwardly biased with respect to, known family relationships, suggesting the presence of additional family relationships (Table 4).The 2 known sib pairs B20016/17 and B20023/24 had estimated levels of relatedness (r ab values) of 0.54 and 0.68, respectively, when calculated only with respect to individuals with the same country of origin.Within GU, average relatedness was low but significantly different from 0, at 0.09 ± 0.04, while within 'VE', average relatedness was significantly higher, at 0.45 ± 0.04.Values were highest within the 'Oregon1' captive population, among whose individuals were many estimates of r ab approaching or above 0.5 (Table 4).The individual from the 'Florida2' captive population, B20021, also had high r ab values with respect to 'Oregon1' birds.Values were higher but similar when calculated with respect to all individuals (data not shown).

Population differentiation
Estimated allele frequencies at AFLP loci varied somewhat depending on assumptions about de viations from HWE, but resulting estimates of gen etic diversity varied only slightly (Table 3b).Furthermore, across the entire range of possible values of F IS , estimates of F ST between GU and 'VE' birds were large and highly significant, varying from 0.15 to 0.24 (Table 3b).
Results were similar but more extreme for mitochondrial loci, with an F ST of 0.709 and a Φ ST of 0.605 (Table 3a).An NMDS plot of AFLP data re flected these 2 distinct and clearly separated clusters, each formed only of individuals of the same country of origin (Fig. 3).
Of 170 variable AFLP loci observed, 65 had alleles that were private to either one putative population or the other, and 4 were diagnostic of a country of origin.The observed number of both private and diagnostic loci was significantly higher than expected by  Bayesian model-based cluster analyses of AFLP data with STRUCTURE also signaled that population structuring was present.According to Evanno's ΔK, genotypes were separated into 2 clusters (Fig. 4), with additional signs that structuring was not an artifact: α never varied more than 0.2 units, only 3 individuals had nonhomogeneous Q scores, assignment was invariant between replicates, and ln Pr(X |K ) was highest for K = 2 (Fig. 4).Individuals with homogeneous Q scores included GU birds in 1 cluster, and in the other they included all individuals from the Oregon1 captive flock as well as the Florida2 individual.The 3 individuals with non-homogenous Q scores, i.e.B20014, B20020, and B20022, all came from different captive flocks (Florida1, Oregon2, California1).
FLOCK decisively identified 2 genetic clusters by producing 49 identical solutions out of 50 runs (pla teau length = 49) with K = 2.One run aborted.This solution allocated 12 out of the 13 'VE' genotypes to one cluster and all 5 GU genotypes to the other.The one apparently misallocated genotype (B20020) was also identified by AFLPOP as very likely not belonging to either the 'VE' or the GU populations (p < 0.002 for both).However, AFLPOP did allocate B20020 to the 'VE' group, albeit with a much lower LLOD score than the other 'VE' genotypes.B20020 was also intermediate in STRUCTURE analyses (Fig. 4) and in the NMDS plot (Fig. 3).Thus, both STRUCTURE and FLOCK identified 2 genetic groups corresponding, with the possible exception of B20020, to 'VE' and GU origins.All simulations in both sets of 100 runs produced MLLOD scores below the MLLOD obtained when reallocating the empirical GU and 'VE' genotypes (see the Supplement).Therefore, we rejected (p < 0.01) the hypotheses that the GU and 'VE' samples were differentiated solely because of a high degree of relatedness among the GU specimens or the 'VE' specimens, and concluded that the 2 samples very likely originated from 2 genetically distinct populations.

Genetic differentiation between red siskin populations
All sources of genetic evidence revealed a level of differentiation between red siskins sampled from Venezuelan-descended captive stock versus wild individuals sampled in Guyana that was consistent with an origin from different populations.Divergence between samples was consistently significant and large as measured by F ST , even when very high le vels of population inbreeding were assumed (Table 3).Clustering analyses consistently separated GU from 'VE' individuals, both using model-based (Fig. 4) and model-free algorithms (FLOCK, AFL -POP); that separation could not be explained by high levels of individual relatedness.Furthermore, measures of differentiation between GU and 'VE' birds based on mitochondrial loci were also large and significant (Table 3).While mtDNA sequence divergence between GU and 'VE' was low on an absolute scale, it was substantial with respect to divergence among the other Neotropical species in the genus, which frequently share mtDNA haplotypes (Beckman & Witt 2015).Although 'VE' and GU birds had contrasting levels of high and low within-population mtDNA diversity, respectively, this is a common phylogeographic pattern in central ('VE') and peripheral (GU) populations.It is furthermore unsurprising given the much larger and fragmented distribution of 'VE' populations, which were not individually sampled for this study, versus the single known GU pop- ulation, which is confined to a relatively small area geographically.What is noteworthy is that the divergence from 'VE' is a shared pattern among all GU samples.Furthermore, diversity within GU and 'VE' for nuclear loci was consistent with measures for similar species with similar markers (Bensch & Akesson 2005), and there was no sign of reduced diversity in nuclear markers in GU samples that might be expected with a recent founder event (Table 3).Thus, although we found no mitochondrial variation in our small sample of GU birds, normal levels of AFLP variation combined with significant differentiation from 'VE' make a recent anthropogenic introduction event in GU seem unlikely.

Possible biases
A comparison of estimated levels of relatedness r ab (Table 4) with the known family relationships among our samples (Table 1) allowed us to weigh the relative importance of several sources of possible bias in our other estimates.Estimates of relatedness between 2 known sib pairs were slightly above their theoretically ex pected level of 0.5, by 0.11 units on average, even when calculated only with respect to birds between which true immediate pedigree relationships were plausible (presuming that short-lived birds in Guy ana could not share immediate family relationships with those from random captive colonies in the USA).Possible causes of this upward bias among known relatives included size homoplasy, additional hidden population structuring, and inbreeding (Wang 2011).Size homoplasy was present, as evinced by the negative and significant correlation between fragment size and frequency (Caballero et al. 2008), and additional hidden population structuring among 'VE' birds may have been present if individuals were drawn from across the historic Ve nezuelan range, from populations with natural substructuring.However, even if both of these caused the upward bias in relatedness among known sibs, additional within-popu la tion inbreeding clearly was present particularly among the Oregon1 captive populations; many pairwise levels of relatedness within populations ex cee ded the ~0.1 level of bias present among known sibs.
Nevertheless, regardless of the level of inbreeding (F IS ) assumed among birds from a given country of origin, estimates of F ST were always large and significant (Table 3).Furthermore, simulations revealed that inbreeding/drift effects alone -even if ex -treme -were insufficient to cause a spurious signal of differentiation (Table S2a,b).Even if size homoplasy was an important source of bias for relatedness calculations, its effect is to cause underestimates of genetic differentiation between populations (Innan et al. 1999, Vekemans et al. 2002, Caballero et al. 2008); this implies that our estimates of F ST would have been even higher had size homoplasy not been present.Additionally, F ST estimates may be biased by non-equilibrium conditions, at least for our 'VE' samples (Crow & Kimura 1970).The Bayesian clustering analyses implemented in STRUCTURE were also vulnerable to this violation of assumptions.
The explanation for nonhomogeneous Q scores in 3 captive individuals from separate locations is currently still unclear.One possibility is true admixture, i.e. that wild red siskins have been traded from Guyana, and not only Venezuela.This would highlight the need for effective continued monitoring of the illegal bird trade in both countries.However, this lack of homo geneity may more likely result from additional, poorly sampled structuring within 'VE' birds, or from a lack of full resolution between GU and 'VE' due to small sample sizes.Indeed, 1 individual (B20022) had the highest number of loci with missing data, although the other individuals had no missing data.
Three final potentially confounding factors in our evaluation of genetic differentiation include the effects of cryptic hybridization with other species, temporal biases, and possible domestication.Cryptic hybridization is potentially the more problematic with our 'VE' samples from US captive flocks, for which very limited information is available about breeding history.Red siskins were historically brought into captivity to hybridize and backcross with canaries to produce red canary varieties (Birkhead 2003); the instability of many of these varieties drives a continued demand for wild individuals (McCarthy 2006).Other species have also been re ported to have been hybridized with red siskins (McCarthy 2006); however, avicultural groups pre sently tend to prize conservation and efforts at purity over those of hybridization (e.g.Porter 2017), and continued trafficking means that wild birds continue to be incorporated into captive stocks (R. Weil 2013, https://www.academia.edu/35413094/El_bachaquero_ilustrado).Understan ding the extent of genetic introgression into captive populations is thus an important area of future research.Temporal biases may exist because birds were sampled at different points in time from the 2 populations.Domestication in 'VE' birds may be a final source of spurious differentiation (Frankham 2008).However, conti n ued incorporation of wild indi viduals would slow the genetic effects of adaptation to captivity.

Possible origins of red siskins in Guyana
In spite of the multiple possible sources of bias outlined above, a signal of genetic differentiation likely due to true population structuring is evident in our data.Our initial hypotheses about the presence of red siskins in Guyana included a recent anthropogenic introduction or a natural origin from a historical vicariant or dispersal event.While our data were insufficient to distinguish definitively between these alternatives, ancillary evidence leads us to favor the latter hypothesis.
An anthropogenic introduction could produce the observed level of differentiation only in 2 scenarios.The first alternative is if it was associated with an extremely severe bottleneck event (i.e. the introduction was a single pair of individuals, with an extremely low population size for many generations).However, an accidentally introduced population would be very unlikely to establish, and more importantly, such a scenario would be expected to create strong and widespread LD, and reduced variation, which we did not observe among AFLP loci.The second alternative is that the GU individuals were descendants of a recent introduction from a previously strongly differentiated Venezuelan population not represented in the 'VE' samples included in the present study.However, the range of red siskins was formerly continuous across northern Venezuela, poaching has been widespread across the entire historic range of the species in Venezuela (Coats & Phelps 1985), and our 'VE' samples included individuals from 5 flocks in 3 US states.These facts give us no particular reason to suspect that major Venezuelan variants would not be represented in captive flocks.A recent anthropogenic origin of the GU population from a natural, differentiated 'VE' population could be examined in the future by analyzing sequence variation in geo referenced museum specimens of wild-caught individuals from across the former Venezuelan range.
The disjunct coastal−inland distribution implied by a natural origin for Guyana birds is not unknown in other bird (and even non-avian) species in the region (Mees 2000, Lew et al. 2006, Lim et al. 2010) although the biogeographic reasons behind these patterns are unclear.Bird species often share some aspect of habitat preference; thus, past expansion and contraction of savanna/forest ecotone habitat areas may tend to produce similar distributional patterns.Fleshing out this hypothesis with possible causal factors will require additional sampling and in-depth study of distribution and natural history of populations in Guyana.

Founders for a captive breeding program aimed at recovering populations in Venezuela
Our evidence is still insufficient to determine definitively the origin and history of the Guyanese red siskin population.However, genetic differentiation is significant between birds of Guyanese and Venezuelan origin by multiple measures, and their disjunct distribution makes ongoing demographic exchange between these countries unlikely.These observations, coupled with preliminary observations of differences in habitat and behavior in the 2 countries (Robbins et al. 2003, J. Miranda, unpublished data), along with the precautionary principle (Groom et al. 2006), support treating these groups as separate elemental conservation units (Wood & Gross 2008) until additional evidence can be produced that de monstrates otherwise.
Future research producing this additional evidence will be crucial, as well as research focused on using captive birds as founders for such recovery efforts, including understanding risks such as possible cryptic hybrid ancestry, inbreeding/drift due to generations at small population size, and domestication.For the time being, however, we conclude that although the wild population of red siskins in Guyana is crucial for mitigating the threat of global extinction of this species in the wild, present evidence points to significant genetic differentiation from birds of Venezuelan origin, making GU birds less attractive for use in a captive breeding program aimed at the recovery of this species in Venezuela.

Fig. 1 .
Fig. 1.Inferred historical distribution of the red siskin Spinus cucullatus (after Robbins et al. 2003) and presently persisting populations, including a recently discovered one in Guyana.Note that the precise number and location of extant Venezuelan populations are indicated arbitrarily to safeguard against potential poaching

Fig. 2 .
Fig. 2. Median-joining network among Spinus cucullatus mtDNA haplotypes (cytochrome B and control region).Substitutions are shown as tick marks on each branch.Numbers within circles are the observed number of individuals with each haplotype; specimen IDs are indicated alongside.Complete ID numbers for the Guyanese samples are available upon request to safeguard the location of the Guyanese population

Fig. 3 .
Fig. 3. Non-metric multi-dimensional scaling (NMDS) analysis of amplified fragment length polymorphism (AFLP) loci genotyped in Spinus cucullatus individuals sampled in Guyana (empty circles) or US captive flocks of presumed Venezuelan origin (filled circles).Complete ID numbers for the Guyanese samples are available upon request to safeguard the location of the Guyanese population

Fig. 4 .
Fig. 4. Bayesian clustering analyses of amplified fragment length polymorphism (AFLP) data from Spinus cucullatus individuals sampled in Guyana (GU) or US captive flocks of presumed Venezuelan origin ('VE'), using model-based algorithms as implemen ted in STRUCTURE.(a) Likely number of genetic clusters as indicated by ΔK (gray line) and mean ln Pr(K) (above-axis numbers, with the highest value indicated in bold), and (b) membership coefficients (Q) in those clusters.Complete ID numbers for the Guyanese samples are available upon request to safeguard the location of the Guyanese population

Table 1 .
Spinus cucullatus specimens used in the present study.'VE': of presumed Venezuelan origin; GU: of Guyanese origin; USNM: US National Museum of Natural History, Washington, DC (USA) largest (508 bp) or the last monomorphic peak, whichever was smallest.Peaks had to conform to the following additional criteria, following Bonin et al. (2007): fluorescence intensity above 100; low baseline fluorescence; a clean negative control; a clear, single base width profile without 'shoulders;' no peak at the same location in the co-loaded PCR product or size standard; no closer than 3 bp from another fragment; and strong sample amplification across the entire size range of fragments.Two coauthors (K.M.R.C. and B.D. a Complete ID numbers available upon request to safeguard the location of the GU population ({Lindenmayer & Scheele 2017)

Table 4 .
Relatedness (r ab ) of each sampled Spinus cucullatus individual, with respect to individuals from the same country.Those from the Oregon1 captive population are highlighted in gray boxes, and relatedness values among all individuals from that population are demarcated by dashed lines.Values for known sibling pairs are highlighted in black boxes; values consistent with first-order relationships (> 0.5) are highlighted in bold a Complete ID numbers available upon request to safeguard the location of the Guyanese population