The introduction of non-native species is a major threat for biodiversity, negatively impacting native ecosystems if they become successful invaders (Estoup and Guillemaud, 2010; Boubou et al., 2012). Nevertheless, not all species introduced into new localities outside their native ranges become invasive (Boubou et al., 2012; Richmond et al., 2014; Cristescu, 2015). Different characteristics of an invasion determine its failure or success (Estoup and Guillemaud, 2010; Jeschke, 2014). Indeed, critical aspects that must be considered to explain the complex process of establishment and spread of exotic species into new habitats include the original population source(s) where propagules (i.e., first introduced individuals, seeds, eggs) come from, the number of introduced propagules, if several introduction episodes occurred, among others (Boubou et al., 2012; Monzón-Argüello et al., 2015).

Tracing the history of an invasion represents a fundamental question for the understanding of the invasion process itself, which provides key information to build predictive models of secondary spread, to evaluate ecological or economic impacts, and to develop control methods (Audzijonyte et al., 2008; Dlugosch and Hays, 2008; Estoup and Guillemaud, 2010). Nonetheless, we need to rely on indirect information when studying invasive species, because, for instance, the exact source where the invasive population originated, whether one or several introduction episodes occurred, is most commonly unknown; this requires a systematic sampling from both native and invaded regions. Hence, a diverse array of genetic markers and analytical methods has been used to reconstruct the evolutionary history of biological invasions, despite their coarse resolution in some cases (Ficetola et al., 2008; Cristescu, 2015). Indeed, genetic information and molecular approaches have proved to be effective for exploring invasion histories of numerous species, successfully corroborating the source(s) of introduction, invasion routes, number of introductions, among others, at different ecological and evolutionary time scales (Lee, 2002; Dlugosch and Hays, 2008; Boubou et al., 2012). Furthermore, phylogeographic studies of native and nonnative populations have been able to precisely trace back invasion routes, which is essential for deciphering the history of the invasion and of the origin and genetic composition of invading populations (Estoup and Guillemaud, 2010; Boubou et al., 2012; Cristescu, 2015).

Snakes have become significant invasive species due to their naturally efficient predatory habits, their capability to withstand long periods with no food, and because some species are susceptible to being transported by humans for the pet trade (Reed and Rodda, 2009; Reynolds et al., 2013). In the last 50 years, several episodes of successful snake invasions have been described for mainland and island ecosystems across the world: the Brown Treesnake Boiga irregularis on Guam (Richmond et al., 2014), the natricine Natrix maura on Mallorca (Guicking et al., 2008), the Californian Kingsnake Lampropeltis californiae on Gran Canaria (Monzón-Argüello et al., 2015), and the Asian Wolfsnake Lycodon capucinus on Christmas Island (Smith et al., 2012). Controversially, while most of these invasive events involve snake species of less concern in terms of their conservation status, remarkable examples in which a vulnerable species turns into invasive have been reported, for example the Burmese python Python bivitattus in the Everglades (Reed and Rodda, 2009; Dorcas et al., 2012) and Boa constrictor on several Caribbean islands (see below). Both are classified as vulnerable in their native ranges, as a result of local harvesting for food, skin for the leather industry, traditional medicine purposes, and the pet trade (IUCN, 2017).

The Boa constrictor complex is the most widely distributed boid naturally present throughout the Neotropics and inhabiting a variety of vegetation types (Bertona and Chiaraviglio, 2003). Recently, phylogeographic and phylogenetic analyses based on mitochondrial and nuclear DNA showed it comprises three putative species: Boa sigma present along the Pacific coast of Mexico (west of the Isthmus of Tehuantepec), Boa imperator distributed throughout the Gulf of Mexico, Yucatán Peninsula and Central America, and Boa constrictor inhabiting most of South America (Suárez-Atilano et al., 2014; Card et al., 2016). The natural distribution of species of Boa on different islands across its latitudinal range has been historically recognized, for instance B. imperator in Crawl Cays, Bay Islands or Cayos-Cochino (Boback, 2005, 2006; Boback and Siefferman, 2010; Green, 2011) and B. sigma in Islas Marías (Casas-Andreu, 1992). However, during the last four decades it has established populations outside its natural range, particularly along the Caribbean, where it has become a successful invasive species: Aruba (Quick et al., 2005; Bushar et al., 2015), Puerto Rico (Reynolds et al., 2013; Bushar et al., 2015), Cozumel (Martínez-Morales and Cuarón, 1999; Vázquez-Domínguez et al., 2012), and St. Croix (Golden, 2017; Reynolds and Henderson, 2018).

Oceanic or ecological islands are particularly vulnerable to biological invasions, mainly due to their high levels of endemism and frequent small species pools and population sizes in comparison with their mainland counterparts (Frankham, 1997; Bellard et al., 2016). Cozumel Island, on the Mexican Caribbean, is an oceanic island (i.e., with no land bridge ever connecting with the mainland) of approximately 478 km², located 17.5 km off the Yucatán Peninsula, in southern Mexico (Cuarón, 2009). The island was formed between the Oligocene and the Pleistocene, separated from the continental mainland by the Cozumel channel, more than 400 m deep, with yearly predominant strong currents. The dominant vegetation type is semi-evergreen tropical forest (Romero-Nájera et al., 2007; Cuarón, 2009). Boa individuals were introduced into Cozumel and were set free on the island after the making of the movie El Jardín de Tía Isabel in 1971 (Martínez-Morales and Cuarón, 1999). Anecdotal information indicates that the introduction included a small number (2–30) of adult individuals. This snake rapidly became the most abundant terrestrial vertebrate on the island (Romero-Nájera et al., 2007). Being a key predator, the boa's presence has been directly linked with the collapse of the populations and the precarious conservation situation of several endemic taxa like the Pygmy Raccoon Procyon pygmaeus, Cozumel Coati Nasua nelsoni, Cozumel Harvest Mouse Reithrodontomys spectabilis, Cozumel White-Footed Mouse Permyscus leucopus cozumelae, Cozumel Thrasher Toxostoma guttatum, Cozumel Curassow Crax rubra grioscomi, amongst others (Martínez-Morales, 1999; Martínez-Morales and Cuarón, 1999; Cuarón et al., 2004, 2009; Romero-Nájera et al., 2007; Espindola et al., 2014). Notably, the phylogenetic relationship of the Cozumel boa within the Boa constrictor complex is still unknown. Vázquez-Domínguez and collaborators (2012) analyzed the genetic diversity of the Cozumel boa populations based on microsatellite loci. They found the island population is structured as two genetic clusters, with moderate levels of genetic diversity and signal of recent genetic bottlenecks. Their results are in agreement with a recently introduced population, founded by a few individuals likely from localities from the Gulf of Mexico and Yucatán Peninsula, rendering Cozumel an ideal system to decipher the history of invasion and to corroborate the anecdotal and genetic information about the invading boa population. Importantly, the absence of boas on Cozumel before and during the Spanish Conquest has been confirmed based on archaeological and historical evidence (Hamblin, 1984).

Hence, our aim was to reconstruct the invasion history of the non-native Cozumel boa population, by deciphering the original source(s) of the Cozumel population, routes of invasion, and likely number of propagules, while also estimating historical genetic and demographic parameters. To this end, we applied a comprehensive set of analytical tools including phylogenetic inference and haplotype networks reconstruction methods, and Approximate Bayesian Computational (ABC) algorithms, combining nuclear and mitochondrial genetic markers. In addition, Álvarez-Romero et al. (2008) indicated that the source of the introduced individuals was from Morelos, Mexico; this specific locality was phylogeographically assigned to the distribution of Boa sigma (Suárez-Atilano et al., 2014; Card et al., 2016); thus, in order to adequately identify the identity of the source populations as a first step, we needed to clarify the genealogical relationships of the Cozumel boa in relation to the main clades present in Mexico (Boa imperator and Boa sigma, Suárez-Atilano et al., 2014; Card et al., 2016). Accordingly, we evaluated the population genetics of B. constrictor (sensu lato) in Cozumel, based on mitochondrial and nuclear sequences and nuclear microsatellite loci, and we included the widest geographic sampling along the species entire distribution. Our premise was that if the Cozumel population was derived from a reduced number of individuals and/or from different source populations, we would find demographic parameters accurately in agreement with this scenario of invasion, namely historical reduced population sizes and bottleneck signals, in combination with contemporary population growth and recent genetic and temporal divergence from the putative source populations.

MATERIALS AND METHODS

Sample collection and molecular procedures.—We selected a subsample of 16 individuals of Boa constrictor from Cozumel (Fig. 1), collected by Romero-Nájera et al. (2007) and genotyped for microsatellite loci by Vázquez-Domínguez et al. (2012). We amplified the mitochondrial cytochrome b gene (cyt b) with the primers L14919 (5′–GACCTGTGATMT-GAAAACCAYCGTTGT–3′) and H16064 (5′–CTTTGGTTTA-CAAGAACAATGCTTTA–3′; Burbrink, 2002) and a nuclear intron of the ornithine decarboxylase gene (ODC) using the primers OD-F (5′–GACTCCAAAGCAGTTTGTCGTCT-CAGTGT–3′) and OD-R (5′–TCTTCAGAGCCAGGGAAGC-CACCACCAAT–3′; Noonan and Chippindale, 2006). Those novel sequences were deposited in GenBank (cyt b: MN325976–MN325990, ODC: MN325991–MN326004). Amplification protocols and PCR conditions were followed exactly as in Suárez-Atilano et al. (2014). For the ODC dataset, alleles were inferred with a maximum-likelihood method implemented in Phase v.2.1.1 (Stephens and Donnelly, 2003).

Fig. 1

Map depicting sampling localities for Boa imperator from (A) continental mainland and (B) Cozumel Island. Squares indicate individuals from the Gulf of Mexico clade (n = 14; see Results and Fig. 2), circles from the Yucatán Peninsula clade (n = 16), and diamonds from Cozumel (n = 16). Colors follow the ODC haplotype network population nomenclature (Fig. 3).

Phylogenetic inference.—In order to establish the genealogical relationship of the boa from Cozumel within the Boa constrictor complex, we performed phylogenetic analyses based on mitochondrial cyt b. For this, we included 254 sequences with a verified geographic origin (Hynková et al., 2009; Reynolds et al., 2013; Suárez-Atilano et al., 2014; Bushar et al., 2015; Card et al., 2016; Tables S2, S3; see Data Accessibility). The dataset encompassed the entire distribution of the complex, from the northernmost portion of Mexico throughout the Caribbean to northern South America. We performed phylogenetic analyses with maximum-likelihood and Bayesian inference methods.

Sequences alignment was done with ClustalW v.2.1 (Larkin et al., 2007) in BioEdit v.7.2.5 (Hall, 1999). We used the Akaike information criterion scores from jModelTest v.2.1.9 (Darriba et al., 2012) to select the best fitting model of sequence evolution for our final dataset (see results below), which was used for the following phylogenetic analyses. Maximum-likelihood was performed with PhyML v.3.1 (Guindon et al., 2010), based on NNI+SPR for branch lengths and topology optimization. Clade support was assessed with 1,000 non-parametric bootstrap replicates. Bayesian inference was done with MrBayes v.3.2 (Ronquist et al., 2012), using default settings and four chains sampled every 1,000 generations for 75 million generations; 25% of generations were discarded as burnin. We estimated the 50% majority-rule consensus topology and posterior probabilities for each node with the remaining trees. We also analyzed the mtDNA dataset with the Bayesian method implemented in Beast v.1.8.4 (Drummond et al., 2012), considering a Yule Process tree prior, the uncorrelated lognormal relaxed molecular clock method, and 50 million generations sampled every 1,000^th generation. Convergence and stationarity for both Bayesian analyses were visualized with Tracer v.1.6 (Rambaut et al., 2014). In addition to phylogenetic tree-based inferences, we constructed a rooted haplotype network using the Neighbor-Net algorithm with SplitsTree 4.14.2 (Huson and Bryant, 2006), based on the Patristic Distance corrected by the GTR+I+G model of evolution. Sequences of Corallus annulatus (KC750011) and Corallus ruschenbergerii (HM348838) were used as outgroups in all phylogenetic analyses.

Considering the phylogenetic position of Cozumel samples within the Boa imperator clade (see Results below, Fig. 2) and with the aim to analyze the genealogical relationships of Cozumel individuals at a population level with respect to the rest of the B. imperator clade, we constructed unrooted networks of unique haplotypes/alleles separately for the mitochondrial cyt b and the nuclear intron ODC, using the Median-Joining algorithm from PopArt v.1.7 (Leigh and Bryant, 2015).

Fig. 2

Graphical representation of the demographical reconstruction based on ABC analyses implemented in DIYABC (Cournet et al., 2014) considering different molecular datasets and priors schemes (see text for a detailed explanation).

Population origin and demography.—For the demographic analyses that follow (see Fig. 3 for a graphical description), we used two datasets which vary in terms of sample size and kind of molecular markers: a) Original Reduced Dataset (ORD): this included seven polymorphic microsatellites, ∼1,060 bp of cyt b, and ∼600 bp of ODC intron, amplified in 16 individuals from Cozumel, 14 from the Gulf of Mexico, and 16 from the Yucatán Peninsula; all these individuals were analyzed by Vázquez-Domínguez et al. (2012) and Suárez-Atilano et al. (2014). b) Polymorphic Microsatellite-based Dataset (PMD): we built a dataset that included a higher sample size (57 individuals from Cozumel, 30 from Gulf of Mexico, and 54 from Yucatán Peninsula), based only on the seven polymorphic microsatellite loci genotyped across all samples.

Fig. 3

(A) Bayesian maximum clade credibility tree based on 220 unique cytochrome b haplotypes for the genus Boa inferred with Beast (Drummond et al., 2012). Black dots represent the main clades obtained. Posterior probabilities and bootstrap support values are shown above and below the diagonal, respectively. (B) Rooted network based on the same 220 unique cyt b haplotypes, inferred with SplitsTree (Huson and Bryant, 2006). Arrows indicate the position of the Cozumel unique haplotypes within Boa imperator. See Data Accessibility for tree file.

The Approximate Bayesian Computation method (ABC; Beaumont et al., 2002), implemented in the software DIYABC v2.04 (Cornuet et al., 2014), was developed to perform comprehensive analyses of population history. It includes an efficient Bayesian model choice framework founded on linear discriminant analyses of summary statistics from both empirical and simulated scenarios. Importantly, it allows evaluating complex population histories, including any combination of population divergence events, admixture events, and past demographic changes (population size changes). This software can be used to compare specific evolutionary scenarios and quantify their relative support, and also to estimate different temporal and demographic parameters for one or more scenarios, at the same time providing a means to evaluate the level of confidence that can be put into the various estimations; it can be based on a variety of molecular markers (including microsatellite and DNA sequence data; Cornuet et al., 2014).

In our study, the ABC approach was used to infer a) the most likely ancestral source(s) of the non-native boa population on Cozumel, b) the time, expressed as number of generations, when the introduction likely occurred, and c) the most probable number of propagules. ABC analyses require one to follow a certain set of steps; first, we established three hypothetical scenarios based on the phylogenetic position of the Cozumel haplotypes within B. imperator (see Figs. 2, 3), and the genetic similarity with mainland populations found with microsatellites in Vázquez-Domínguez et al. (2012): Scenario 1, the Cozumel population originated from an admixture of individuals from the Gulf of Mexico and Yucatán Peninsula; Scenario 2, it originated solely from Gulf of Mexico populations; or Scenario 3, it originated solely from the Yucatán Peninsula populations.

Next, in order to assess the best possible optimization of priors to test those scenarios, we did a preliminary DIYABC run testing only the genealogical origin of the Cozumel population, without any assumption about bottlenecks or effective number of founders (‘dummy trial'). This run included a uniform prior distribution of probabilities across the three scenarios and the default parameters regarding range of population size (10–10,000 individuals) and divergence times (10–10,000 generations ago). We generated 3×10⁶ simulated datasets; the posterior distribution of parameters from this preliminary computation was used to set priors for the next set of simulations, which are compared with our empirical dataset, considering the following:

Relaxed simulation: We tested a specific range of priors (based on the dummy trial posterior ranges), where we set the invasive population as sampled t1 generations ago, a bottleneck duration of db generations before reaching a larger effective population size (t1 > tdb) of 2 to 100 generations, and an effective number of founders (Nb) between 2 and 50 individuals. We assumed a uniform prior distribution of probabilities for the three scenarios, such that each had roughly equal representation in a reference table, which consisted of 6×10⁶ simulated datasets (Tables 1, 2).

Table 1

Posterior parameter estimations from the best-fit Approximate Bayesian Computation (ABC) demographic model (Scenario 1) for Boa imperator, based on the Original Reduced Dataset (ORD), which includes seven polymorphic microsatellite loci, ∼1,060 bp of cyt b, and ∼ 600 bp of the nuclear ODC intron amplified over 16 individuals from Cozumel, 14 from the Gulf of Mexico, and 16 from the Yucatán Peninsula.

Table 2

Posterior parameter estimations from the best-fit Approximate Bayesian Computation (ABC) demographic model (Scenario 1) for Boa imperator, based on the Polymorphic Microsatellite-based Dataset (PMD), which includes seven polymorphic microsatellites loci genotyped over 57 individuals from Cozumel, 30 from the Gulf of Mexico, and 54 from the Yucatán Peninsula.

Constrained simulation: For this, we included information considering the typical generation time of three years per generation reported for continental boa populations (Reed and Rodda, 2009). To directly test this, the time since population establishment (t1) was set to 2–18 generations and the bottleneck duration (db) to 1–10 generations. These values are set considering ∼40 years since introduction and a bottleneck of ∼28 years (from time of introduction in 1971 to time of frequent sightings in 1999), that divided by three years per generation yield the upper values within the range used. The rest of the parameters were the same as those used for the relaxed simulations (Tables 1, 2).

Next, for both the relaxed and constrained final simulation schemes, we compared the fit of the three scenarios by estimating their posterior probabilities: with the obtained reference tables from each scenario, we ranked the simulated datasets in order of increasing distance to the observed data considering direct and logistic approaches (Beaumont et al., 2002; Cornuet et al., 2014). Distance between datasets was based on normalized Euclidean distances between the summary statistics estimated from the empirical and the simulated sets. The summary statistics (within-population and between-populations pair statistics) for the microsatellites loci included: the mean number of alleles (N_A), mean expected heterozygosity (H_E), mean allelic size variance (Var_AS), mean Garza-Williamson's M, F_ST, shared allele distance (D_S), and Goldstein's genetic distance (δµ²). Summary statistics for sequence-based markers were: Tajima's D, number of haplotypes (h), number of segregating sites (Sn), mean of the pairwise differences (k), variance of the number of pairwise differences (Var_k), and private segregating sites. For the microsatellites dataset, we used a stepwise-mutation model (Estoup et al., 2001), while the HYK model of evolution (Hasegawa et al., 1985) was used for cyt b and ODC data (Table 1).

Finally, we performed a pre-evaluation step using a principal components analysis (PCA) to ensure that at least one (or more) scenarios and priors combination would produce simulated datasets close enough to the empirical data, allowing for model comparisons. The PCA was based on the space of parameter values and summary statistics for 5,000 simulated datasets, generated from the prior distributions of parameters specified in the different scenarios. We assessed the performance of the parameters of the three scenarios estimating type I (false positives) and type II (false negatives) rate errors, by computing the relative bias across those 5,000 simulated data sets.

Bayesian assignment analyses.—With the aim of assessing the admixture proportions between the Cozumel and the inferred source mainland populations (see Results), we performed a Bayesian assignment analysis using Structure v.2.3.4 (Pritchard et al., 2000), including 30 individuals from the Gulf of Mexico, 54 from the Yucatán Peninsula and 14 from Cozumel, in which cyt b and microsatellites were amplified. The ODC intron was not used due to the relatively low resolution and the limited number of samples in which this was sequenced (Fig. S2; see Data Accessibility). We calculated the probability of individual assignment into population clusters (K) without prior information on the origin of individuals (POP-FLAG = 0); correlated allelic frequencies and admixture and non-admixture model was set for microsatellites and cyt b, respectively. This method handles both nuclear and mitochondrial data, where the former relies on allele frequencies and the latter on polymorphic sites. We conducted several tests using different numbers of population clusters (Pr[X|K], K = 1 to K = 5), with length value set at 1,000,000 iterations after a burn-in period of 100,000 MCMC repetitions and ten simulations per value of K. Each test yielded a log-likelihood value of the data (lnprobability for K) that was used to estimate the maximum number of clusters with the Evanno's ΔK test (Evanno et al., 2005) available in Structure Harvester v.0.6.94 ( http://taylor0.biology.ucla.edu/struct_harvest/). To evaluate the geographic distribution of the assignment patterns, we reclassified the probabilities of individual assignment for those individuals with P > 0.8 (K for each marker separately) into categories of assignment to the different clusters and interpolated those values across uniformly spaced ∼10 km grids. We used the Inverse Distance Weighted interpolation procedure (IDW; Watson and Philip, 1985) in the Spatial Analyst extension of ARCGIS v.10.2.1 (ESRI), based on spatially unique points (65 for microsatellites and 80 for mtDNA); the interpolation was masked within the predicted geographical distribution of B. imperator (Suárez-Atilano et al., 2017).

RESULTS

Phylogenetic inferences.—We successfully amplified a ∼1063 bp fragment of cyt b from the 16 Cozumel individuals. We found two haplotypes (h = 0.143), with a nucleotide diversity of p = 0.00054 and an average number of nucleotide differences between haplotypes of k = 0.57. The combined dataset including Boa imperator and Cozumel had 220 unique haplotypes, for which the TIM2+I+G (pinv = 0.48; α = 0.7) was the selected evolutionary model; the closest equivalent, GTR+I+G, was used in all subsequent analyses. The estimated Tr/Tv ratio was 3.6, with mean nucleotide frequencies of A: 34.64%, C: 29.15%, G: 11.33%, T: 24.89%. Regarding the ODC intron, we amplified a ∼611 bp fragment from our 16 Cozumel selected individuals, which yielded five alleles (h = 0.581), nucleotide diversity of p = 0.0021, and average number of nucleotide differences between alleles of k = 1.29.

All phylogenetic analyses based on mtDNA recovered a similar and highly supported topology, in which the Boa constrictor complex was split into three main clades as suggested by Card et al. (2016): Boa constrictor, Boa imperator, and Boa sigma. Also, all haplotypes from Cozumel were nested within the B. imperator clade, encompassing a well-supported clade with samples from the Yucatán Peninsula (Fig. 3A). A similar position for the Cozumel haplotypes, closer to the Yucatán Peninsula clade, was found in the rooted network (Fig. 3B).

The haplotype network based only on cyt b haplotypes from the B. imperator clade showed that the most abundant haplotype from Cozumel is also widely distributed across Yucatán, Belize, and Guatemala (Fig. S1; see Data Accessibility). On the other hand, the network based on 41 unique ODC alleles showed that the four inferred alleles are present in Cozumel (Fig. 4). The most abundant haplotype in Cozumel was also widely distributed across Campeche and Guatemala, whereas the other three also occur in Veracruz, Tabasco, and Puebla (Figs. 1, 4).

Fig. 4

Neighbor-net haplotype network based on 41 unique nuclear ODC haplotypes within Boa imperator.

Approximate Bayesian computations.—Model comparisons for the Relaxed simulation based on the Original Reduced Dataset (ORD) indicated that Scenario 1, in which the Cozumel population originated from a mixed contribution from the Gulf of Mexico and Yucatán Peninsula individuals (Fig. 4A), had the highest probability values both for the preliminary (direct P = 0.29; logistic P = 0.56) and final (P = 0.54; P = 0.92) computations (Fig. S3H, J; see Data Accessibility). Pre-evaluation of model scenarios via PCA showed that the empirical (observed) data fall within the cluster of simulated data points based on the three demographic scenarios (Fig. S3D; see Data Accessibility). The confidence test for scenario choice showed that Type I error (i.e., probability in which the true scenario is rejected) was 0 for the logistic approach of the final constrained run; Type II error values (i.e., probability of accepting a scenario when it is not true), based on scenarios 2 and 3 simulated datasets, were 0.23 and 0.15, respectively, for the logistic approach of the final run. Finally, the goodness-of-fit test using the posterior predictions indicated that parameter values and summary statistics from the simulated datasets based on Scenario 1 closely matched the empirical data (Fig. S3F; see Data Accessibility).

Model comparisons achieved with the Relaxed simulation and Polymorphic Microsatellites-based Dataset (PMD) also showed that Scenario 1 had the highest probability values both for the preliminary (direct P = 0.32; logistic P = 0.56) and final (direct P = 0.52; logistic P = 0.58) computations (Fig. S3I, K; see Data Accessibility). Also, the pre-evaluation of model scenarios via PCA indicated that the empirical (observed) data fall within the cluster of simulated data points based on the three demographic scenarios (Fig. 5E). The confidence test for scenario choice showed that Type I error was 0.06 and values of Type II error were 0.56 and 0.11 (simulated scenarios 2 and 3, respectively), while the goodness-of-fit test also indicated that parameter values and summary statistics from the simulated datasets based on Scenario 1 closely matched the empirical data (Fig S3G; see Data Accessibility).

Fig. 5

(A–C) Main demographic scenarios modeled with DIYABC (Cournet et al., 2014) based on two molecular datasets and the Constrained Optimization scheme of Boa imperator from Cozumel, the Yucatán Peninsula, and the Gulf of Mexico. Note that time is not to scale. (D–E) Prior distribution of parameters corresponding to the three invasion scenarios tested, where the first two axes are shown. (F–G) Posterior predictive test of model fit corresponding to Scenario 1. (H–I) Direct and (J–K) logistic approaches of the three scenario model comparisons. In all cases, left panels (D, F, H, J) and right panels (E, G, I, K) show results from the Original Reduced Datasets (ORD) and the Polymorphic Microsatellite-based Dataset (PMD), respectively.

The demographic estimations for the Relaxed simulation were similar between datasets, where results indicated for ORD and PMD, respectively, an average invasive population size (Nb) of 31.5 (CI 95%: 8.2–57.5) and 37.5 (CI 95%: 8.22–73.6) individuals, a mean time of introduction into Cozumel (t1) of 53.9 (CI 95%: 14.3–95.3) and 46.9 (CI 95%: 17.6–81.5) generations, with a bottleneck duration (db) of 16.8 (CI 95%: 3.56–29.5) and 23.2 (CI 95%: 11.3–46.34) generations (Table 2).

Model comparisons based on the Constrained simulation (Fig. 5) also indicated the best model selection for Scenario 1, for both ORD (direct P = 0.36; logistic P = 0.58) and PMD (direct P = 0.39; logistic P = 0.64) datasets. The confidence test for scenario choice showed that Type I and Type II error values were 0.075 and 0.66 for ORD, and 0.055 and 0.62 for PMD. The goodness-of-fit test exhibited similar parameter values and summary statistics from the simulated datasets based on Scenario 1 as those observed with the Relaxed simulation. The ORD computation results (Table 1) estimated an average invasive population size (Nb) of 31 individuals (CI 95%: 26.3–35.5), the mean time of introduction (t1) as 15 generations (CI 95%: 9.78–18), with bottleneck duration (db) of seven generations (CI 95%: 6.98–7.98). For the final estimation based on PMD (Table 2), the estimated average invasive population size (Nb) was 28.4 individuals (CI 95%: 22.1–29.4), the mean time of introduction into Cozumel (t1) 17.1 generations ago (CI 95%: 12.3–18.0), with bottleneck duration (db) of 6.94 generations (CI 95%: 6.87–8.0).

Bayesian clustering.—Structure results were similar for both markers, where the maximum number of clusters detected by the Evanno's test was K = 3 for microsatellites [Ln Pr(K = 3) = –3011.94] as well for cyt b [Ln Pr(K = 3) = –1990.33]. The K = 3 results showed a similar clustering pattern, in which the Cozumel population represents a single cluster that is well differentiated from the mainland populations. Notably, the admixture results based on the mitochondrial gene indicate that Cozumel has higher ancestry from regions of the Yucatán Peninsula (Fig. 6B, D); in comparison, the microsatellites exhibit higher ancestry from Tabasco (Gulf of Mexico), Yucatán, and Quintana Roo (Fig. 6A, C).

Fig. 6

Bayesian clustering for two main continental clades of Boa imperator (Gulf of Mexico and Yucatán Peninsula) and the Cozumel population using Structure (Pritchard et al., 2000), based on the resulting K = 3 for (A) microsatellite loci and (B) cytochrome b haplotypes. Maps show spatial interpolations of individual assignment into the three genetic clusters obtained with (C) microsatellites and (D) cyt b. Circles indicate individuals from the Gulf of Mexico clade (n = 30), triangles from the Yucatán Peninsula clade (n = 54), and diamonds from Cozumel (n = 16).

DISCUSSION

The presence of populations of this snake on Cozumel has been considered controversial due to the island's geographic proximity to the Yucatán Peninsula and the boa's natural distribution across the Yucatán mainland and several Caribbean islands, suggesting that its occurrence on the island could result from a natural expansion (Álvarez-Romero et al., 2008). Nonetheless, the historical (Hamblin, 1984), ecological (Martínez-Morales and Cuarón, 1999; Romero-Nájera et al., 2007), and genetic (Vázquez-Domínguez et al., 2012) evidence has shown signatures of a recently founded population. Indeed, a natural origin of the Cozumel population, which would have occurred either via rafting of individuals from the nearby mainland (Scenario 3) or rafting from mainland sources farther away (Scenario 2), was refuted by the human-mediated introduction (Scenario 1). Hence, our study supports the introduced origin of this species on Cozumel Island, unraveling specific aspects of the invasion history that were unknown or only anecdotal.

Notably, we were able to confirm two fundamental hypotheses proposed to explain the historical spread and success of the boa on Cozumel, namely its origin from a relatively reduced number of propagules and its admixed origin from multiple localities. Despite the fact that the shallow molecular variation found prevented the identification of specific source-localities, we detected that the Cozumel population of Boa originated from an admixture of individuals from lineages phylogeographically recognized as Boa imperator (Gulf of Mexico and Yucatán Peninsula; Suárez-Atilano et al., 2014). A multiple-sources propagules scenario is not a common account when considering other non-native boa populations (see Quick et al., 2005; Reynolds et al., 2013; Bushar et al., 2015) or even other invasive snake species on islands (Smith et al., 2012; Richmond et al., 2014; Monzón-Argüello et al., 2015; Silva-Rocha et al., 2015). In those examples, authors found that propagules originated from a single source or lineage, showing essentially no genetic variation for different nuclear and mitochondrial markers.

Phylogenetic history and genetic diversity patterns of the Cozumel boa.—Our results show key phylogenetic patterns from the mitochondrial and nuclear genes, where both confirm that the Cozumel individuals belong to the B. imperator clade, and the nuclear information supports multiple source populations. The Cozumel individuals are included within B. imperator, a species naturally distributed across the Yucatán Peninsula, Gulf of Mexico, and Central America (Suárez-Atilano et al., 2014; Card et al., 2016). We found two cyt b haplotypes on the island: one, most abundant, is also widely distributed in the Mexican states of Yucatán and Quintana Roo, as well as in Belize and Guatemala (i.e., the Greater Yucatán Peninsula), while the other haplotype is exclusively present in our sampling within the island. Both the low cyt b haplotype diversity (h = 0.143) and the low average number of nucleotide differences between haplotypes (k = 0.57) suggest a founder effect, as a result of a small pool of founders, with haplotypes little differentiated by four transitions among sequences.

Invasive populations of Boa occurring on other islands have been shown to have close and well-resolved phylogenetic relationships to proximal coastal populations, low haplotype diversity, and low genetic divergence, similar to the patterns detected in our study. The population in Aruba has only three cyt b haplotypes closely related with Brazilian populations, exhibiting a genetic divergence lower than 1.35% (Bushar et al., 2015); also, the only two cyt b haplotypes found in Puerto Rico belong to B. constrictor (South America) and B. imperator (Panamá), respectively (Reynolds et al., 2013). On the other hand, boid species naturally distributed within the Caribbean region exhibit higher levels of haplotype diversity and genetic divergence within their populations; for example, the Jamaican Boa Chilabothrus subflavus (h = 0.79, F_ST = 0.76; Tzika et al., 2008, 2009), the Puerto Rican Cave Boa C. inornatus (h = 0.92, ϕ_ST = 0.20; Puente-Rolón et al., 2013), the Turks Island Boa C. chrysogaster (h = 0.79, ϕ_ST = 0.61; Reynolds, 2011), and B. imperator across Bay Islands and Cayos-Cochino (11 haplotypes; Green, 2011). Native vertebrate species on Cozumel show deep mitochondrial divergences between the island populations and their continental counterparts, like Procyon pygmaeus (Flores-Manzanero, 2014) and Crocodylus acutus (Pacheco-Sierra et al., 2018).

In contrast with the mitochondrial results, four alleles were found in Cozumel with the nuclear ODC gene. One allele is distributed across the Yucatán Peninsula, while the other three have a wide mainland distribution on the Mexican states of Veracruz, Tabasco, and Puebla, localities included in the Gulf of Mexico clade in our phylogenetic results (Fig. 3). Considering the higher substitution rate of the mitochondrial gene in relation to the nuclear loci (Avise, 1994), a higher nuclear than mitochondrial diversity is unusual within boids. For instance, Reynolds (2011) found nine different cyt b haplotypes and only two alleles based on the NT3 and c-mos nuclear loci in the native Turks Island Boa C. chrysogaster. Nonetheless, it would not be unexpected that for introduced populations like those on Cozumel, the nuclear gene with its lower substitution rate has likely retained the diversity (ancestral polymorphism) of the original (source) individuals/populations (from the states of Campeche, Tabasco, Veracruz, and Puebla). At the same time, cyt b diversity could have been lost due to the genetic bottleneck resulting from introduction, although this would need to be tested with a wider sampling from the mainland.

Bayesian clustering/assignment methods are frequently used to evaluate admixture proportions and source populations, based on the principle that the invasive population will cluster with one of the potential source populations (Roman and Darling, 2007; Rius and Darling, 2014; Cristescu, 2015). If no clustering is observed, it is likely because individuals from the invasive population share ancestry with various population sources (evidence for an admixture origin of the invasive population), unsampled source populations, a severe drift process during and after introduction, or insufficient number of markers (Estoup et al., 2001; Darling et al., 2008). Our assignment results with both mitochondrial and nuclear loci delimited the Cozumel boa population as a differentiated genetic cluster from the Yucatán Peninsula and Gulf of Mexico, suggesting Cozumel has a particularly defined genetic composition, distinct from that of the sampled mainland populations. The admixture origin evidenced by the nuclear and mitochondrial genes is further supported by the assignment results, showing in more detail that regions from the Yucatán Peninsula (cyt b), and from Tabasco, Campeche, Yucatán, and Quintana Roo (microsatellite loci) are likely original sources of the Cozumel population (Fig. 6). This pattern of nuclear ancestry from the southernmost portion of the Gulf of Mexico region, also solved by the ODC intron based network (Fig. 3), was confirmed in the clustering analysis based on microsatellites (Fig. 6A, C).

Invasion history of the non-native boa population in Cozumel.—It is recognized that coalescent-based simulation analyses (including ABC simulations), which are based on neutral theory, cannot always accurately reconstruct human-mediated invasions (Estoup and Guillemaud, 2010). Hence, they have to be used with caution due to the recent and restricted timescales of most introduction processes, typically much smaller than the timescale of mutation and genetic drift (Fitzpatrick et al., 2012). Also, invasive species do not always experience the expected reduction in genetic diversity, as some studies have shown (Dlugosch and Parker, 2008; Hundertmark and van Daele, 2010; Richmond et al., 2014; Cristescu, 2015). On the other hand, ABC specifically is considered adequate for complex demographic and historical models where bottlenecks, serial introductions, and/or genetic admixture are present, as in biological invasions (Estoup and Guillemaud, 2010).

Our ABC simulations (Fig. 5, S3; see Data Accessibility) significantly exhibited that the scenario in which the Cozumel population derived from an admixture of individuals from the Yucatán Peninsula and the Gulf of Mexico was the most statistically supported. Importantly, the geographically mixed origin scenario was consistent regardless of prior optimization strategies (Relaxed and Constrained) or dataset analyzed (ORD and PMD), even though the latter varied both in terms of sample size and of number and kind of genetic markers. It has been noted that increasing the number of markers, high levels of polymorphism within markers, and selection of an adequate set of summary statistics when comparing empirical and simulated scenarios has a greater effect in the probability of selecting the true demographic model than increasing the number of individuals in a sample (Aeschbacher et al., 2012; Gautier et al., 2013; Robinson et al., 2014; Buzbas and Rosenberg, 2015). Hence, our results show that the genetic information included in the original dataset (microsatellite loci, nuclear, and mitochondrial sequences) likely counteracted the relatively small size; on the other hand, the higher sample size and high level of microsatellite polymorphisms in the PMD was enough to support the same demographic scenario. Furthermore, the occurrence of a population bottleneck was strongly supported when changes in population size were considered in the simulations, thus confirming previous findings based solely on microsatellites (Vázquez-Domínguez et al., 2012), which showed a strong signal of a contemporary bottleneck, suggesting a few individuals were introduced, and also that these originated from multiple geographical sources.

Our relaxed and constrained ABC simulation schemes were statistically appropriate to simulate and evaluate our empirical data; namely, both schemes resulted in equally good-fit estimations for the ORD and PMD molecular datasets. Although the temporal estimations obtained with the relaxed simulation were three times higher than the constrained, when the results for number of generations were transformed into year-based time, considering 1–1.5 years per generation for the former and three years for the latter simulations, results were remarkably similar. That is, for the constrained simulation, in which we are explicitly considering three years per generation (Reed and Rodda, 2009) and, in accordance, a maximum number of generations of 18, the Cozumel population merge time (t1) is 45–51 years ago and a bottleneck duration (db) of ∼21 years (Tables 1, 2). These estimates are concordant with the information about the boa's invasion: April of 1971 is when boa individuals were released on two different sites on Cozumel Island (Martínez-Morales and Cuarón, 1999; Vázquez-Domínguez et al., 2012). The simulations estimated differences between the genetic markers that may be related with their coalescence times, where microsatellites have higher effective population sizes and coalescence times than mtDNA (Zink and Barrowclough, 2008). In addition, the time needed for the population's establishment, the time from introduction to when individuals started to reproduce regularly (i.e., within their generation time) and spread across the island (lag phase), must be considered; thus, the lower estimates are likely closer to the introduction date. Our estimations of the bottleneck are also in agreement with the time at which boas actually became conspicuous and common to see across the island (1990s).

A potential flexibility in generation time after invasion needs be considered, which could result from different factors: a) the arrival of sexually mature boas and likely gravid females to the island, favoring a rapid population expansion (Vázquez-Domínguez et al., 2012) and also affecting the bottleneck duration; b) the capability of the introduced adult propagules of reproducing repeatedly during their lifetime (e.g., iteroparity; Hughes, 2017); furthermore, those propagules could be the only reproductive individuals during the early years after their arrival on Cozumel; c) overlapping of generations; d) the time lag for the mating of individuals born in Cozumel and delayed contribution to the gene pool of the island population; e) the absence of predators (naturally present on mainland); and f) a continuous source of food. We also need to take into account that generation time in snakes has been typically difficult to measure in natural populations (Rossman et al., 1998), jointly with the ample parameter ranges (they are not exact numbers) that are obtained by ABC analyses (Estoup and Guillemaud, 2010). Accordingly, for the relaxed simulation in which the maximum number of generations were estimated by the Bayesian optimization of DIYABC, an average boa generation time of 1–1.5 years matches the known historical and ecological evidence of the boa's arrival on Cozumel (Martínez-Morales and Cuarón, 1999; Cuarón et al., 2009), resulting in a Cozumel population merge time (t1) of 47–81 years ago and 17–35 years for the bottleneck duration (db).

Regarding demographic estimates, the average number of propagules we obtained with both the relaxed (31.5–37.5 individuals, ODC and PDM, respectively) and the constrained simulations (31.2–26.4) were highly concordant with the number of individuals being released on Cozumel (2–30; Martínez-Morales and Cuarón, 1999). Hence, although it is difficult to establish with absolute certainty the number of propagules and population sizes, based on anecdotal (historical) estimates (Martínez-Morales and Cuarón, 1999), preliminary genetic information (Vázquez-Domínguez et al., 2012), and the fact that several adult animals were used during the filming (as shown in scenes in the movie), we consider our genetic evidence to clarify and reconcile the ecological and historical record.

Interestingly, our results suggest complementary scenarios based on the mitochondrial haplotypes and the nuclear OCD alleles for the potential sex contribution to the original propagules. In one scenario, a few closely related females from the Yucatán Peninsula were introduced along with a few more unrelated males from a wider geographic origin, including diverse localities across the Yucatán Peninsula and Gulf of Mexico. Another scenario is that the females were gravid during introduction, had stored viable sperm from previous mating events, and produced neonates soon after release (Vázquez-Domínguez et al., 2012). A final possibility is differential reproductive success. However, a higher sample size from within the island population and further analyses are needed to corroborate these scenarios.

This historical genetic distinctiveness of the Cozumel boa population is likely due to the differentiation effects of genetic drift that result from the founder effect it experienced. Serial introductions or introductions from multiple sources may mitigate the drastic reduction of genetic diversity expected in recent invasive populations (Kolbe et al., 2007; Dlugosch and Hays, 2008; Jeschke, 2014; Richmond et al., 2014), a factor also often correlated with successful establishment due to the combined benefit of maintaining a high population size and high genetic diversity (Estoup et al., 2001; Cristescu, 2015). In Cozumel, just over two decades after its introduction, the boa drastically increased its population numbers, rapidly spreading across the entire island (Martínez-Morales and Cuarón, 1999), while currently still exhibiting moderate levels of genetic diversity. Moreover, the fine spatial scale study within Cozumel we previously performed (Vázquez-Domínguez et al., 2012) allowed us to identify an emergent pattern of genetic structuring on the island, where the boa population shows two distinct genetic clusters associated with the presence of anthropogenic landscape features, specifically main roads and urban settlements. Hence, we have been able to discern both the historical process of invasion and the contemporary genetic structuring of the boa in Cozumel. The propagules' multiple-sources origin and the unique genetic composition of this invasive population rendered a genetic signature, detectable and traceable by different phylogenetic and clustering methods. Similar patterns of genetic structuring at different ecological and evolutionary scales have been found in other invasive vertebrates, like the black rat Rattus rattus in Senegal (Konečný et al., 2013) and R. rattus and the house mouse Mus musculus in Cozumel (Borja-Martínez, 2017).

Boa's duality in conservation terms.—Our results confirm the non-native origin of Boa imperator, where a reduced number of propagules, belonging to different localities from mainland Mexico, were intentionally released on Cozumel from the early 1960–70s. Information based on a geographically extended sampling of putative original ranges from both native and invaded distributions, in combination with a robust methodological approach and use of molecular markers with different inheritance and polymorphism models (Boubou et al., 2012; Cristescu, 2015), allowed us to discern the invasion history of the boa population. We provide relevant information to be considered for the management of the Cozumel boa, taking into account the recent partition of Boa into three differentiated species (Boa constrictor, B. imperator, B. sigma; Suárez-Atilano et al., 2014; Card et al., 2016). The latter furthermore concerns the status of this boa because of its duality in conservation terms: as a threatened species protected in Mexico (Semarnat, 2010) and other countries, which nationally hinders the possibility of control or eradication programs, but at the same time, an introduced predator threatening a global conservation priority site (e.g., an Alliance for Zero Extinction, Endemic Bird Area, and Key Biodiversity Area site; Ricketts et al., 2005; IUCN, 2016; BirdLife International, 2017) and a critically endangered insular ecosystem. Notwithstanding, we need to work with the Mexican authorities in order that our evidence of the boa's introduced nature and its negative impact on the endemic and native biota is taken into account, hence opening the possibility that this critically endangered species—forbidden by law to be managed in any way—can be controlled or eradicated from Cozumel. Notably, the boa's admixed origin limits the utility of Cozumel individuals for conservation breeding programs for the species and hinders the possibility of reintroduction or translocation of individuals into the population sources, strategies that have been suggested as potential actions.

Finally, it is important to highlight that the analytic approach we followed can be applied to track propagule sources, infer routes of invasion, and design management and conservation actions for different invasive species systems, including scenarios where the information about the invasion history is incomplete, anecdotal, or blurred despite a recent occurrence.

DATA ACCESSIBILITY

Supplemental material is available at  https://www.copeiajournal.org/cg-18-102.