Levels of genetic differentiation and gene flow between four populations of the Scaly-naped Pigeon, Patagioenas squamosa: implications for conservation

ABSTRACT Island-endemic columbid species are particularly vulnerable to environmental degradation, extreme climatic events, and interactions with exotic species. The situation might be even more critical in the case of exploited species, where legal hunting and poaching can severely affect population dynamics. Here we document for the first time the genetic structure of the Scaly-naped Pigeon, Patagioenas squamosa, a Caribbean-endemic columbid species of cynegetic interest, over a large part of its range. Using both mitochondrial DNA and nuclear markers (microsatellites), we investigated gene flow, genetic diversity, and genetic structure among four islands populations originating from Puerto-Rico, Guadeloupe, Martinique and Barbados. We found evidence for a significant genetic differentiation only between the Barbados and the three other populations, consistent with the fact that the Barbados population originated from a few captive individuals escaped from a rooftop aviary in Bridgetown about 100 years ago. Given the absence of genetic differentiation between Puerto Rico and the French Antilles, our results suggest that, apart from Barbados, the species may mainly consist of a single large, homogeneous population. We discuss the relevance of our findings in relation to management and conservation issues.


Introduction
It is now widely accepted that genetic diversity is crucial for the persistence of wildlife populations as low diversity results in both low individual fitness and reduced adaptability of populations in the face of environmental changes (Frankham 2005;Evans & Sheldon 2008). Accordingly, over the last 30 years, molecular and conceptual tools from population genetics and phylogeography have been increasingly used in conservation biology (Soulé & Mills 1998;Hedrick 2001;Russello et al. 2020). Population genetics and phylogeography are particularly helpful to determine the evolutionary and geographic boundaries of species, subspecies and populations (Haig et al. 2006(Haig et al. , 2011Garnett & Christidis 2007;Haig & D'Elia 2010). Defining such units is of prime importance to assess extinction risk, to prioritize conservation actions and identify the relevant spatial and temporal scales at which they should be conducted (Haig et al. 2011). This is especially true for species occurring in archipelagos, for which different populations may face different selection pressures because of contrasted anthropogenic, biotic or abiotic influences between islands (Clegg et al. 2002;Illera et al. 2007;Hoeck et al. 2010).
In addition, island endemic species are often characterized by restricted habitat range and smaller population sizes resulting in lower genetic diversity than mainland species, such that they are more exposed to inbreeding depression and, ultimately, face a higher risk of extinction (Frankham 1998(Frankham , 2005Brooks et al. 2002). Consequently, ecological disturbances and environment changes that are globally occurring, such as deforestation, urbanization, introduced pathogens, and competition and predation pressures from invasive species, particularly affect island endemic species (Walsh et al. 2012;Treglia et al. 2013;Ferrer-Sánchez & Rodríguez-Estrella 2014;Ortega et al. 2015;Dornburg et al. 2016;Palmas et al. 2017;Turvey et al. 2017). This is particularly true of bird species, with 80% of all avian species having become extinct since 1600 being island endemics (Manne et al. 1999). The situation might be even more critical in the case of exploited species, where legal and illegal hunting can severely affect population dynamics (Juillet et al. 2012;Carvalho et al. 2015). In that context, genetic information, including genetic diversity and genetic structure, is of prime importance for delimiting meaningful management units, and developing relevant conservation policies of island game bird species accordingly.
Endemic columbid species can play an important ecological role in island ecosystems through their capacity to disperse seeds over long distances (Shanahan et al. 2001;Bucher & Bocco 2009). However, they are particularly vulnerable to environmental degradation, extreme climatic events, predation, and hunting pressure (Walker 2007;Carvalho et al. 2015;Collar 2015;Stirnemann et al. 2018). Preservation and restoration of forest habitats, sustainability of hunting pressure and control of introduced predators are therefore the main priorities for the management of island columbid species (Walker 2007). However, at a regional scale, information about the degree of connectivity between different island populations is also necessary to understand population structure (Young & Allard 1997;Monceau et al. 2013), movements (Strong & Bancroft 1994), and the role of columbids in the intra-and inter-island transfer of seeds and ecosystem functioning (McConkey et al. 2004;Buelow et al. 2018).
Here we document the extent of genetic structure and gene flow in the Caribbean-endemic Scaly-naped Pigeon (Patagioenas squamosa) using both mitochondrial DNA and nuclear markers (microsatellites). The species' natural range includes both the Greater and the Lesser Antilles and extends to Florida in the north and to the islands off the coast of Venezuela in the south (Del Hoyo et al. 1997). Although the species is a year-round resident throughout much of its range (Nellis et al. 1984;Levesque et al. 2011;García-Quintas & Isada 2014;Rodríguez Batista et al. 2014;Madden et al. 2015), it occurs only seasonally on some islands (Paice & Speirs 2010) and is a vagrant on Jamaica. In addition, individuals can move between islands, particularly following extreme climatic events (Rivera-Milán 1995). The Scaly-naped Pigeon was introduced in Barbados in the early twentieth century, after about 25 individuals escaped from captivity from a rooftop aviary in Bridgetown (Buckley et al. 2009). This arboreal and frugivorous pigeon has a diverse diet and is supposed to play an important role in forest regeneration (Pérez-Rivera 1978). Although the species is currently considered of least concern in the IUCN red list, data on population size, demographic trends and movements of individuals between islands are scarce. Still, the species is exposed to intensive hunting pressure over a large part of its distribution area and is under the threat of ongoing habitat loss and fragmentation (Wiley 1985;Brooks et al. 2002;Acevedo & Restrepo 2008). According to Raffaele et al. (1998), the species has largely declined in the West Indies, except in Puerto-Rico thanks to the regeneration of secondary forests from abandoned agricultural lands (but see Case & Hughes 2011). In addition, climate change in the Caribbean, with a likely increase in extreme climatic events, such as hurricanes (O'Brien et al. 1992;Wiley & Wunderle 1993;Boose et al. 2004) and an overall rainfall reduction (Bhardwaj et al. 2018), may also put the species at risk in the near future. Therefore, using both mitochondrial and nuclear makers, we assess the extent of genetic differentiation and that of gene flow between four different islands spread over the geographical range of the species. In addition, we compare our results with previous estimates of the genetic structure of other columbid species of contrasted conservation status.

Sample collection and DNA extraction
Biological samples were obtained between 2009 and 2012 from four different islands: Puerto-Rico (PR), Guadeloupe (GUA), Martinique (MAR) and Barbados (BAR) (Figure 1 and Table S1). Birds in GUA and MAR were collected only during the hunting season, such that birds were not necessarily local breeders. Depending on the island, various types of biological samples were taken from live individuals (blood, feathers) or from individuals shot by local hunters (liver, legs). Liver samples and blood were placed in a storage buffer (70% ethanol and 30% Tris-EDTA buffer pH 8), whereas legs were deep-frozen and feathers were stored in envelopes with silica-gel.
From liver and toe samples, DNA extraction was carried out from small pieces (about 1 cm2) of tissues. In the case of feathers samples, the calamus tip (2-3 mm in length) was used. For each blood sample, 100 µL were centrifuged at 4000 rpm, at 4°C, during one minute in order to separate the plasma (the supernatant) from the red blood cells. Then, all samples were overnight incubated with 200 µL of Queen Lysis Buffer (Tris 10 mM, EDTA 10 mM, NaCl 10 mM, n-lauroylsarcosine 1%; Seutin et al. 1991) and 5 µL of Proteinase K at 55°C. DNA was thereafter extracted according to a standard phenol-chloroform method, as described by Hillis et al. (1996).

Mitochondrial DNA amplification and sequencing
We amplified a mitochondrial DNA (mtDNA) gene coding for the Cytochrome c Oxidase subunit I (COI) from 111 scaly-naped pigeons, by using reverse and forward primers described by Kerr et al. (2007): BirdF1 (5ʹ-TTCTCCAACCACAAAGACATT GGCAC-3ʹ) and VertebrateR1 (5ʹ-TAGACTTCTGG GTGGCCAAAGAATCA-3ʹ). The 50 µL PCR mixture contained 200 µmol/L dNTP, 1X Buffer (HotMaster TM Taq Buffer, 5 PRIME), 0.25 U Taq (HotMaster TM Taq DNA Polymerase, 5 PRIME) and 200 nmol/L R and F primers. The PCR reaction, performed within a Biorad DNA Engine Peltier Thermal Cycler (Bio-Rad Laboratories, Inc.), began by an initial denaturating at 94°C for 1 min. 30 s, following by 34 thermal cycles consisted of: 94°C for 30 s 53°C for 45 s, and 65°C for 45 s, ended by a final extension at 65°C for 10 min. The quality of amplifications was assessed on 2% agarose electrophoresis gel. Thereafter, PCR products were purified by adding two enzymes: Exonuclease I (EXOI) at a final concentration equals to 2 U per PCR product, and Shrimp Alkaline Phosphatase (SAP) at a final concentration equals to 1 U per PCR product. These enzymes were incubated for one hour at 37°C to be activated, and were deactivated through an incubation at 80°C during 10 min. All samples were externally sequenced by Macrogen Inc. (Republic of Korea) on an Applied Biosystems 3730xl DNA Analyzer sequencer (Applied Biosystems), carried out according to a protocol of Big Dye Sequencing.

Mitochondrial data analysis
Sequences were manually aligned using MEGA 5 software (Tamura et al. 2011). Molecular diversity indices such as haplotype diversity (h), average number of nucleotide differences between haplotypes (π), number of polymorphic sites (S), Nucleotide diversity (k), and pairwise island differentiation (Φ ST ) were estimated with ARLEQUIN 3.5 software (Excoffier & Lischer 2010). Φ st is a F ST analogue for mtDNA sequences and corresponds to the molecular distance between populations based on the number of pairwise differences between haplotypes. These values were estimated with 50 000 permutations and 5% nominal level was adjusted for multiple comparisons with Benjamini-Yekutieli step up procedure (BY; Benjamini & Yekutieli 2001;Narum 2006), resulting in α BY = 0.0204. A minimum spanning network was constructed through POPART 1.7 software (Bandelt et al. 1999;Leigh & Bryant 2015) based on pairwise absolute distances between haplotypes. Neutrality and demographic history tests such as Tajima's D test (Tajima 1989) and Fu's F S test (Fu & Li 1993) were performed, as well as the more conservative square differences statistic (SSD, mismatch distribution), to compare the observed distribution of the number of nucleotide differences between haplotype pairs to that expected under the null hypothesis of sudden expansion, using 10 000 bootstraps (Rogers & Harpending 1992;Harpending et al. 1993;Fahey et al. 2012). These analyses were also performed with ARLEQUIN 3.5.

Microsatellite data analysis
MICROCHECKER 2.2.3 was used to identify potential genotyping errors due to null alleles, scoring errors due to stuttering, and large allelic dropout, using 10,000 iterations and Bonferroni correction (Van Oosterhout et al. 2004). The number of alleles per locus (N a ), observed heterozygosity (H 0 ), expected heterozygosity in Hardy-Weinberg equilibrium (H E ; Nei 1978) and the deviation from Hardy-Weinberg equilibrium (HWE; Guo & Thompson 1992) were performed with Arlequin 3.5, considering 1,000,000 Markov chain steps and 100,000 dememorisation steps. Inbreeding coefficient (F IS ) was calculated in FSTAT 2.9.3 (Goudet 1995) by independently taking each island and at the regional scale. Since some loci included null alleles (Table S2), islands paiwise F ST value were first computed using FREENA (Chapuis & Estoup 2007). However, corrected values of pairwise F ST , which excluded null alleles, were not significantly different than those including them (Table S3). Uncorrected values of pairwise F ST were therefore used to visualize pairwise island differentiation, using FSTAT 2.9.3. Nominal significance level (5%) was adjusted with BY's correction.
Population structure was also visualized using an individual-based Bayesian clustering method, as implemented in STRUCTURE 2.2.3 (Pritchard et al. 2000), and using principal coordinate analysis (PCoA), as implemented in GENALEX 6.41 (Peakall & Smouse 2006). STRUCTURE assigns individuals into K optimal region groups according to their genetic similarity. To that end, and to determine the optimal number of genotype clusters (K), we ran simulations using the admixture and correlated allele frequencies models by programming ten independent runs for K values ranging from 1 to 6 and a burn-in length of 30,000 followed by 100,000 Markov chain Monte Carlo. Thereafter, the optimal number of genotype groups (K) was estimated using the R (version 3.6.2; R Core Team 2019) package CORRSIEVE (version 1.6-8; Campana et al. 2011) by following the ad hoc statistic ΔK method described by Evanno et al. (2005). Clusters were finally visualized using DISTRUCT 1.1 (Rosenberg 2004).
Scaly-naped pigeons were not historically present in Barbados but were imported by local pigeon fanciers (Buckley et al. 2009). Therefore, we used BOTTLENECK 1.2.02 (Piry et al. 1999) to detect distortions in the distribution of alleles frequencies due to a recent founder effect (Luikart et al. 1998), considering the single-step mutation model (SSM), known to be more appropriate for microsatellites loci (Di Rienzo et al. 1994;Piry et al. 1999).

Mitochondrial polymorphism
We obtained 624 bp long COI sequences for the 111 samples of Scaly-naped Pigeon, revealing five haplotypes (H 1 -H 5 ) associated with polymorphism at four single nucleotides (Tables 1 & Tables 2). Haplotype and nucleotide diversities at the local and regional scale are presented in Table 2. The haplotypic network ( Figure  2) presents the relationship between the five different mitochondrial haplotypes detected across the 111 studied pigeons. GUA had the highest haplotype richness (H R = 4), followed by MAR (H R = 3). PR and BAR presented only one haplotype (H 1 ), which was the most abundant at the local and regional scale (Tables 1 &  Tables 2 and Figure 2). Exclusive haplotypes were found in both GUA (H 3 and H 4 ) and MAR (H 5 Figure 2). Pairwise island mtDNA Φ ST values (Table 3), revealed differentiation between BAR and both GUA and MAR, and between MAR and PR.
We first carried out independent demographic scenario tests for islands showing polymorphism (GUA and MAR), and thereafter carried out the same tests at the regional scale. Results from these tests were in accordance with a demographic stability for MAR (Tajima's D = −0.0697, P = 0.239; Fu's F S = −0.6240, P = 0.224 and SSD = 0.00870, P = 0.046). For GUA, results from Fu's F S test supported a within-island demographic expansion (Fu'F S = −2.194, P = 0.021), whereas both Tajima's D and mismatch distribution tests were in agreement with demographic stability (Tajima's D = −1.350, P = 0.069 and SSD = 0.00228, P = 0.032). All of the demographic scenario tests clearly supported a demographic expansion at the regional scale (Tajima's D = −1.544, P = 0.015; Fu's F S = −4.919, P = 0.001 and SSD = 0.000322, P = 0.101).
All island pairwise F ST values for BAR were significantly different from zero (F ST ranging from 0.0496 to 0.0594), whereas all other values were not, suggesting a moderate level of differentiation between BAR and the three other islands (Table 3). This result was in accordance with results based on STRUCTURE analysis ( Figure 3) and PCoA (Figure 4). In both analyses, GUA, MAR and PR seemed to form a homogenous group, whereas BAR appeared to be genetically isolated. PCoA result was based on the first and second axes that, respectively, explained 22.66% and 21.23% of the overall variance in microsatellite data.
Based on BOTTLENECK analysis, we found that BAR was significantly in deficiency of heterozygosity under the single-step mutations model (Wilcoxon rank test: P = 0.0391; Table 4). At the regional scale, no excess or deficiency of heterozygosity was highlighted.

Genetic diversity
Based on mtDNA analysis, we found an overall relatively low haplotype (h = 0.154) diversity across the four Caribbean islands populations, compared to other island-endemic columbid species that was estimated using the control region gene.    represent results from simulation to determine the optimal number of genotype clusters (K) according to the method describes by Evanno et al. (2005). Error bars represents standard deviations over ten runs. (c) represents results from STRUCTURE analysis carried out with maximum-likelihood centered at K = 2 genotype clusters. Each individual is represented by a single vertical bar, partitioned into K colored segments that represent its estimated population membership fractions. Individuals are grouped by population on the x-axis. Population codes are given in Figure 1.

Genetic structure
Genetic diversity metrics based on both mtDNA and nuclear markers were indicative of genetic depletion in Barbados compared to the three other studied populations, suggesting that this population experienced a recent historical demographic event such as a bottleneck or founder event, a more or less longterm geographical isolation or a combination of these phenomena. This hypothesis was supported by the BOTTLENECK analysis (conducted using the stepwise mutation model), which provided evidence for heterozygosity deficiency, suggestive in the present case, of a founder effect. This recent historical demography event broadly explained the observed genetic structure. Indeed, Barbados was strongly differentiated from Guadeloupe and Martinique based on haplotype frequencies (0.112 and 0.242, respectively). In contrast, Barbados was not differentiated from Puerto-Rico, since both populations exclusively consisted of the same haplotype (H 1 ). For the same reason, Puerto-Rico was differentiated from Martinique, whereas, using microsatellite markers, differentiation was only observed between Barbados and the other islands. This pattern was in accordance with results obtained with STRUCTURE and the PCoA analyses, from which two genetic groups were identified, separating individuals from Barbados from those from Martinique, Guadeloupe and Puerto-Rico. This genetic structure could be explained by: (i) the emergence of the Barbados population from a few individuals imported by local pigeon fanciers (i.e., founder effect), (ii) the geographic context of the Barbados island, which is relatively isolated compared to other islands from the Antillean Arc (i.e. restricted gene flow) and (iii) genetic drift, which may be intensified by the two previous phenomena.
In contrast, Martinique, Guadeloupe and Puerto-Rico were not differentiated or structured, what suggests that these islands could be inter-connected by effective migrants and regular movements. This pattern contrasts to some extent with what has been observed in two other Caribbean-endemic columbid species. In  the Zenaida Dove, limited gene flow was observed between eight Caribbean islands, even at very short distance like, for instance, between Martinique and Saint-Lucia (Monceau et al. 2013). Interestingly, the latter study also observed a strong divergence between Barbados and the other islands. In the congeneric Plain Pigeon, P. inornata, restricted gene flow was observed between Puerto-Rico and the Dominican Republic, although the extent of genetic differentiation between the two populations suggests that they may correspond to two sub-species (Young & Allard 1997). However, the observed pattern of genetic structure is in accordance with several other studies on columbids, which highlighted the peculiar aspect of the evolution of this family in a biogeographic context of insularity (Santiago-Alarcon et al. 2006;Kirchman & Franklin 2007;Seki et al. 2007;Goldberg et al. 2011;Ando et al. 2014). Indeed, most columbid species have conserved a strong flying ability contrary to most other island-endemic terrestrial birds (Bollmer et al. 2005;Petren et al. 2005;Kawakami et al. 2008). For instance, the White-crowned Pigeon, P. leucocephala, the phylogenetically closest species to the Scaly-naped Pigeon (Johnson et al. 2010), is known to perform daily movements between breeding areas and foraging sites by flying over 50 km, and can fly up to 150 km over water to reach another Caribbean island in order to track food resources and/or suitable environmental conditions (Wiley 1979 (Wiley & Wunderle 1993;Ries et al. 2018). Hurricanes are known to affect avian community and population dynamics by mainly disturbing feeding and breeding habitats (Rittenhouse et al. 2010;Jenouvrier 2013). Due to the inherent strong ability to fly of most of Columbidae, the high frequency of hurricane events in the Caribbean region may induce movements of individuals, potentially homogenizing the allelic pool over the different island populations, as observed in the present study. Therefore, if movements of individuals between islands are frequent and philopatry is low, it is possible that the whole species actually consists of a single large population.
At the geological timescale, results from historical demographic analyses suggest a recent demographic expansion of the Scaly-naped Pigeon population at the regional scale. The Scaly-naped Pigeon might have had a restricted range during the glacial maxima, while population expanded after the Last Glacial Maximum. The same pattern of expansion was found for the New Zealand Pigeon (Wallis & Trewick 2009;Goldberg et al. 2011).

Implications for the management of the species
The Scaly-naped Pigeon is facing several threats such as habitat loss, legal and illegal hunting, exotic predators and climate change. There is a relatively high hunting pressure within its range of distribution, but monitoring data are lacking (Raffaele et al. 1998;Latta et al. 2010;Latta 2012;Rivera-Milán et al. 2014). In Puerto-Rico, where studies are particularly concentrated, it has been argued that the species requires the development of a sustainable harvest strategy at the local level (Case & Hughes 2011). However, given our results and the strong flying ability of Columbiformes, the insular Caribbean (to the possible exception of Barbados) could be considered as an evolutionarily significant unit (ESU) or management unit (MU; Moritz 1994;Paetkau 1999;Zink 2004), emphasizing the need for conservation effort to be undertaken at the regional rather than local scale. To that end, further sampling could be carried out in other Caribbean islands to confirm the genetic structure of the species at the regional scale evidenced in the present study. In addition, other methods could be used to distinguish between historic and contemporary gene flow, and to identify recently diverged population for an effective conservation management of the species (Paetkau 1999;Palsboll et al. 2006;Whitehead 2010;Esteban et al. 2016). For instance, Paetkau (1999) criticized the unique use of genetic data for MUs recognition and advised to combine it with movement data such as Mark-Recapture, telemetry, or GPS tracking data. Using such methods would indeed be valuable to detect contemporary migration/dispersal due to, for instance, current anthropogenic disruptions (Esteban et al. 2016).

Lay summary
• No genetic differentiation between Puerto-Rico, Guadeloupe and Martinique. • Possibly one large population of Scaly-naped Pigeon at the Caribbean regional scale. • Significant founder effect in the Barbados population, consistent with its introduced origin.
to capture Scaly-naped Pigeons in Barbados, and Nicole Atherley and Sébastien Motreuil for help in catching and processing birds. We thank Francisco J. Vilella for providing samples from Puerto-Rico, and Patrick Asselin de Beauville, Jean-François Maillard, Anthony Levesque, Georges Tayalay (Fédération des chasseurs de Martinique) and local hunters for providing samples from Martinique and Guadeloupe. We also thank Pierre Leclercq for help with genetic analyses at an earlier stage of the study.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
CC was funded by a doctoral grant from the Conseil Régional de la Guadeloupe and Caribaea Initiative.