Local dispersal pathways during the invasion of the cactus moth, Cactoblastis cactorum, within North America and the Caribbean

Cactoblastis cactorum, a species of moth native to Argentina, feeds on several prickly pear cactus species (Opuntia) and has been successfully used as a biological control of invading Opuntia species in Australia, South Africa and native ruderal Opuntia species in some Caribbean islands. Since its introduction to the Caribbean its spread was uncontrolled, invading successfully Florida, Texas and Louisiana. Despite this long history of invasion, we are still far from understanding the factors determining the patterns of invasion of Cactoblastis in North America. Here, we explored three non-mutually exclusive explanations: a) a stepping stone model of colonization, b) long distance colonization due to hurricanes, and/or c) hitchhiking through previously reported commercial routes. Genetic diversity, genetic structure and the patterns of migration among populations were obtained by analyzing 10 nuclear microsatellite loci. Results revealed the presence of genetic structure among populations of C. cactorum in the invaded region and suggest that both marine commercial trade between the Caribbean islands and continental USA, as well as recurrent transport by hurricanes, explain the observed patterns of colonization. Provided that sanitary regulations avoiding human-mediated dispersal are enforced, hurricanes probably represent the most important agent of dispersal and future invasion to continental areas.

During the last decades, biological invasion studies have strongly benefited from the use of neutral molecular markers to disentangle routes of invasion [1][2][3] . This has helped to identify source populations, frequency of invasion events, and the geographical patterns and demographic consequences of invasive species spread 4,5 . Reconstruction of past events of invasion can help in the recognition of the mechanisms of dispersal to prevent the prevalence of invasion or further spread to additional areas 1,6-8 . However, understanding local dispersal remains a central challenge to prevent and control the economic and biodiversity costs of biological invasions.
Initially, molecular markers were used to identify sources and frequency of dispersal events of invasive species to non-native regions [1][2][3]9 . Whereas in some cases, these patterns explain the lower genetic variation of invasive species relative to native areas 1,10-12 , combination of migration events from different source populations have sometimes increased genetic variation within invaded regions, potentially increasing the risk of fast adaptation to novel conditions [13][14][15][16] . Once invasive species are established within non-native areas, understanding further local dispersal is a major challenge to identify potential environmental barriers and recommend control management programs 7,17 . It is possible to use highly variable molecular markers for understanding local scale dispersal and the entangled nature of the species invasion phenomena.
The cactus moth, Cactoblastis cactorum, is an oligophagous herbivore during its larval stage that consumes the inner tissue of cladodes of plants within the genus Opuntia, negatively affecting the plant survival 18 . This moth, native to South America, was used in 1924 for biological control purpose against exotic Opuntia species with successful results in Australia. Later, was intentionally introduced to South Africa (1933), New Caledonia (1933) DnA extraction and microsatellite analyses. We performed total DNA extraction with DNEasy blood and tissue kits (Qiagen, MD, USA) following the manufacturer protocol. Population genetic variation was determined using nuclear microsatellites developed by Genetic Marker Services (http://www.geneticmarkerservices. com). We removed loci that did not amplify among the 29 primers pairs tested. We chose 14 potential polymorphic loci (Table 2) and labeled with fluorophores (Applied Biosystems) for fragment analysis. Each multiplex PCR mixture (5 µL) contained 2.5 µL RadyTaq (Qiagen cat. 206143) and 1 µL DNA template (20 ng) 5pmol for each fluorescent labeled primers. PCR were performed through touchdown reaction starting with initial heat activation at 95 °C for 10 min followed by 6 denaturation cycles of 94 °C for 1 min, annealing for 1 min and 1 min of extension at 72 °C; annealing cycle temperature began at 60 °C or 57 °C and decreased 1 °C every cycle. The PCR reaction ended with two stages of 12 cycles each (57° and 56°, respectively) and final elongation of 72 °C for 5 min. The PCR product was diluted and run on an ABI 3730xl automated capillary sequencer. The allele size was manually scored using a Liz 600 size standard (Applied Biosystems) in GeneMarker (V2.2.0) (Soft Genetics LLC, State College, Pennsylvania, USA).
Basic population estimators of genetic variation. Genetic variation within sampled populations was characterized using the mean number of alleles per locus (NA), the allelic richness (AR), the mean expected-(He) 34 and observed heterozygositiy (Ho), the Fixation Index (F is ) and estimation of F ST 35 between each pair of sampling site using FSTAT 2.9.3. 36 . Mean allelic richness (AR) was calculated using the rarefaction method of Leberg 37 . Exact test for deviation from Hardy-Weinberg equilibrium were performed with GENEPOP 38 . We used FreeNA to determine the frequency of null alleles using EM algorithm 39 . Pairwise population differentiation was tested using only those loci that were in H-W equilibrium in more than 50% of the populations.
Genetic clustering of populations was examined using STRUCTURE 40 (v 2.3.3). We chose the admixture model with correlated allele frequencies. Two analyses were performed, one including all native and invaded population samples, and one including only populations from the invaded area. The first analysis examined whether genetic differentiation has occurred after the human-mediated journey of C. cactorum from its native Argentina to the Caribbean (1957). The second analysis was performed to explore genetic structuring within the invaded Caribbean only. Because groups of larvae were collected from discrete sites, sampling locations were used as prior information 41 . Each run consisted of a burn-in period of 10 5 Markov Chain Monte Carlo (MCMC) iterations, followed by 10 6 iterations. 20 replicated runs were carried out for each value of the potential number of clusters (K) set, between 1 and 9. STRUCTURE HARVESTER 42 (available at http://taylor0.biology.ucla.edu/structure-Harvester/) was used to collating output results from STRUCTURE and to determine the uppermost level of www.nature.com/scientificreports www.nature.com/scientificreports/ structure using the method of Evanno et al. 43 . Because STRUCTURE assumes the absence of null alleles and H-W equilibrium for all loci, analyses were performed using 10 out of 14 microsatellites that fulfilled the required assumptions.
isolation genetic and connectivity in the invaded region. The hypothesis that genetic differentiation between a pair of populations is a linear function of the geographic distance between them was examined through a Mantel test that correlates genetic and geographic distance matrices. Paired genetic distances were estimated as [F ST / (1-F ST )] and geographic distances were estimated as the logarithm of the Euclidean distances (b log ) using GenAlEx 44 (v.6.4). Isolation was also tested considering that hurricanes may also shaped the genetic structure of populations in this area. To explore the relationship between genetic distances and the frequency of hurricanes incidence across pairs of populations, a resistance matrix was constructed using the CIRCUITSCAPE program 45 (Fig. 1). Based on a previous subdivision of the geographic area using a grid with the probability of incidence of hurricanes (categories 1 to 5) following NOAA [46][47][48] . Pairs of populations connected by cells with high incidence of hurricanes were assigned a high conductance value (i.e., low resistance to migration). This rationale was applied to all pairs of populations to construct a resistance matrix (indicated by light areas in Fig. 1). To correlate the resistance matrix against the genetic distance matrix controlling for the geographic distance among populations, a partial Mantel test with 10,000 permutations was implemented.
Approximate Bayesian computational analysis. To provide a quantitative evaluation of the dispersal routes of the cactus moth into Florida (USA), an Approximate Bayesian Computational Analysis was performed with DIYABC 49 . Analyses were performed with ten microsatellites that were at Hardy-Weinberg equilibrium. This approach uses genetic information provided by microsatellites data and historical data given by dates of first observation within invaded regions. Three hypothetical invasion scenarios were examined (Fig. 2), and parameter values were drawn from prior distributions (Table 3). In all three cases we used the native Argentinian population of Ayuí as the initial source of invasion. This population was genetically close in average with the invasive populations (Paired F ST = 0.225). Since moths were taken from Argentina to Australia, and then to South Africa, the analysis considered a "ghost" population representing a non-sampled population that may represent this step. Given the location and date of introduction to the Caribbean, the sampled population of Antigua was considered the source population that initiated the invasion in the region. Antigua was preferred over Nevis because of the larger number of individual larvae in our collection, the high genetic similarity between both islands (Paired F ST = 0.012, N.S.) and the short difference in time of introduction between them (1 year). The following step in the invasion process was characterized by the population La Romana (Dominican Republic), given that the moth   www.nature.com/scientificreports www.nature.com/scientificreports/ was detected in the island of Desecheo, located between Puerto Rico and Dominican Republic in 1963. After this step, hypotheses differ in how the cactus moth reached Florida. In the first hypothesis, Cuba, represented by the sampled population of Santiago, is the source of the Floridian invasive population (Fig. 2, Scenario 1). Due to US embargo to Cuba, any dispersal of the moth from Cuba to Florida is more likely related to environmental agents (such as hurricanes) rather than commercial transportation. The second hypothesis considers that moths belonging to La Romana (Dominican Republic) directly invaded Florida (Highlands) (Fig. 2, Scenario 2). This route constitutes the recorded trajectory of commercial traffic to Florida of ornamental cacti infested with C. cactorum from Dominican Republic and Puerto Rico detected in 1989 50 . However, we cannot rule out the possibility that hurricanes also contributed to dispersal among these areas. The third hypothesis considers that the population of Highlands (Florida) was the result of the admixture of genotypes belonging to Cuba and Dominican Republic (Fig. 2, Scenario3).
Given that C. cactorum usually has two generations per year 27,28,51 , this generation time was used to scale coalescent time in all scenarios. Dates of introduction to Australia (Tg) and Antigua (T2) were given as fixed values based on historical data of the intentional introduction events 18 . In the other areas (T 1 (Argentina), T 3-5 (Dominican Republic, Cuba and USA); Table 3), a minimum and maximum value was assigned to each time www.nature.com/scientificreports www.nature.com/scientificreports/ parameter, with the dates of first observation of the insect as the lower boundary (Table 3). We simulated three million microsatellite data sets (one million for each scenario). Subsequently, the probability for each scenario was inferred by polychotomous logistic regression on the 1% of the simulated data sets closest to the observed data set 49 . The selected scenario was chosen as that with the highest posterior probability value. To evaluate the robustness of the analysis, we used pseudo-observed simulated data sets to quantify the type I error rate (risk to exclude the focal scenario when it is the true one) and the type II error rate (risk to select the focal scenario when it is false) 49 . All simulations and ABC analyses were carried out in DIYABC (v.2) sofware 52 .

Results
No evidence of linkage disequilibrium was detected for all 14 microsatellite loci. Four loci were excluded from the analysis because they had more than 20% of null alleles and were not at Hardy-Weinberg equilibrium (i.e., cc7b, cc16, cc65, cc6b). With the remaining ten microsatellites loci we found that the native population from Yuquerí (Argentina) had on average the highest number of alleles (3.8), while the lowest average number of alleles (1.6) was observed in one of the most recent invasive population (Lee Co., USA) ( Table 1). Allelic richness (AR) was higher in native range (2.016-2.185) than invaded populations (1.453-1.992; Table 1). A deficit of heterozygotes was detected in one of the native populations (Ayuí), and in three of the Caribbean invaded populations (Antigua, Nevis and Santiago; Table 1). On average, lower levels of heterozygosity were recorded for invaded than native populations (Table 1).
Pairwise comparisons indicated significant differentiation between the native Argentinian populations and all invaded populations in the Caribbean and Florida (mean F ST = 0.229, range = 0.128-0.49). Within the invaded region, genetic differentiation ranged from almost no differences to rather high levels of genetic differentiation ( Table 4). The average level of differentiation within the invaded region was F ST = 0.185. In general, early invaded populations (Nevis and Antigua) presented lower levels of differentiation with the rest of the populations in the region than recently invaded populations. The Jamaican population (Palisadores) presented the greatest level of genetic differentiation within the Caribbean (mean F ST = 0.310; Table 4), whereas Florida populations had a low level of differentiation with population La Romana (Dominican Republic) and low to moderate differentiation with two of the closest Cuban populations (Pinar del Río and Trinidad; Table 4).
Clustering genetic analyses using STRUCTURE on the whole dataset revealed the presence of two main groups (K = 2, with the delta K method) distinguishing native from all invasive populations (Fig. 3B), and tend to confirm a single introduction event. When analyzing only the group of invasive populations, K = 2 is also selected (Fig. 3A). While individuals from Puerto Rico, Palisadores and Pinar del Rio belong to a well-defined cluster, the rest of the populations showed an admixed pattern. Analyses of isolation by distance and connectivity including only the invading populations indicated a low but almost significant level of isolation by distance (r = 0.165; P = 0.05). On the other hand, the isolation by resistance hypothesis based on hurricane and wind currents was significantly supported (r = 0.435; P = 0.008).
Results from Approximate Bayesian Analyses (ABC) support quiet well the hypothesis described in the third scenario (Fig. 2), due to (1) a higher consistency between simulated and sampled data (i.e. the highest posterior probability; P = 0.7798), (2) non-overlapping confidence intervals, and (3) low Type I and mean Type II errors (0.038 and 0.000; respectively). This dispersal scenario considers that the invasion to Florida (represented by the Highlands population) was funded by insects belonging to at least two sources of dispersal, one from Dominican Republic and the other from Cuba (represented in the analysis by the Santiago population; Fig. 2). The posterior probability estimation of parameters of model three assigned a median admixture rate of R = 0.4.  Table 3. Parameters used for data simulation in the three competing scenarios using ABC analyses.

Discussion
During the 1981-1991 decade, North American customs authorities intercepted Opuntia plants infested with C. cactorum. Seizures occurred in plant shipments sent from the Dominican Republic to Miami and in luggage to Dallas Texas International Airport 53 . Commercial trade of ornamental cacti from the Dominican Republic to the United States of America coupled with the constant influx of American tourists to the Caribbean islands lead to the proposal that human-mediated dispersal has been the major agent of migration to the continent. Similarly, the comparison of the mitochondrial haplotypes from eastern populations of C. cactorum of Florida and those of Dominican Republic supports this idea 21,31,32 . Despite other potential mechanisms of dispersal as island-hopping and tropical storms or hurricanes 30 , the present study support the genetic similarity and likelihood of migration between La Romana (Dominican Republic) and Highlands (Florida) populations. However, the distribution of genetic variation and ABC analyses of invasion routes within the Caribbean region also suggest the presence of other nonexclusive routes of dispersal.
In the Caribbean and North Atlantic region, hurricanes are a ubiquitous temporary phenomenon with consistent wind current patterns. Studies of genetic variation have shown that hurricanes and marine currents influence the patterns of dispersal and constitute a source of variation in ecological and demographic processes 54 . Wind has been identified as one of the most important factors promoting the spread and long-distance dispersal   www.nature.com/scientificreports www.nature.com/scientificreports/ of invasive insect species 7,55 and to enhance distance and/or speed of dispersion as in the case of the monarch butterfly 56 , the cattle egret 57 or the mealy bug, Oracella aculata 58 . Also, the grasshopper Eumetopina flavipes, vector of the sugarcane virus Ramu, entered Australia from New Zealand favored by wind currents 59,60 . Wind is also responsible for long-distance dispersal of the wasp Megastigmus schimitscheki in France 7 . Although previously suggested, the role of hurricanes on local dispersal of C. cactorum within the Caribbean was not examined until recently 32 . Our previous results using mitochondrial COI 32 coupled with the findings of this study using nuclear microsatellites, are consistent and provide new pieces of evidence supporting the role of hurricanes on the spread of C. cactorum toward mainland North America. In comparison with the isolation by distance model, the isolation by resistance migration approach was a much better predictor of the current pattern of genetic variation of C. cactorum in the Caribbean and Florida. This finding suggests that hurricanes were one of the main drivers of dispersal from the Caribbean to mainland North America as previously shown using a more conserved mitochondrial marker 32 . In addition, ABC analyses of the invasion routes supported our hypothesis that moths entered Florida both through commercial transportation from the Dominican Republic and from Cuba as a consequence of climatic events (hurricanes); whereas trade with Cuba was suspended as a result of the embargo. Moreover, this finding suggests that hurricanes may have also dispersed C. cactorum from its initial introduction in Nevis, Antigua and Montserrat to other islands onto the Caribbean region. Historical records of hurricanes in the region (https://www.nhc.noaa.gov/) further support this statement. After the initial human-mediated introduction of C. cactorum, three hurricanes connected the Lesser Antilles and the Dominican Republic between 1963 and 1967, four impacted the south of Cuba and the Dominican Republic between 1975 and 1980, and two passed from Cuba and impacted the USA in 1985 and 1987. Nevertheless, the presence of the moth in the lower Bahamas in 1983 61 , and the possibility of island hopping, as recorded in the Hawaiian archipelago, can also represent a possible route of dispersal towards Florida. During the last decades at least three hurricanes impacted the Bahamas before reaching the Atlantic coast of Florida (Mitch in 1988; Andrew in 1992, and Katrina in 2005). Hence a more extended survey of populations in this area of the Caribbean will provide further insight on the dispersion routes across the region.
The process of invasion is generally associated with the loss of genetic variation due to successive demographic bottlenecks as the invasion of new locations proceeds 62 . This process can be counteracted if repeated invasion events occur in a given location, thus increasing the amount of genetic variation to equivalent, or larger levels of variation, than that observed on the native region 17,63 . Our analyses revealed that populations of C. cactorum from the Caribbean and Florida, sustained less genetic variation than in the native region, likely as a consequence of genetic drift. Given that present populations of C. cactorum outside the Caribbean region are located more that 8,000 km away in other biogeographical areas and continents, the genetic structure results support previous analyses using mitochondrial markers suggesting that the insect entered the Caribbean probably once in 1957 21 . Similar to the Lepidopteran tomato pest, Tuta absoluta, in Africa and the Mediterranean Basin 62 , the occurrence of just one introduction of the cactus moth was enough to allow its spread throughout the Caribbean despite the lower genetic variability than the source population from South America. Unlike the tomato pest and other invasive insect species 1,9 significant genetic differentiation was detected in the invaded region suggesting again a possible effect of drift. A more complete survey of invaded populations in Australia, South Africa and in other continents will help to decipher whether the magnitude of reduction in genetic variation will affect the adaptive potential of the moth into other continental areas with different climatic conditions as those of the Caribbean and Florida.
The combined effect of hurricanes and genetic drift reducing genetic variation through space will likely affect the evolutionary potential of the moth during future introductions to the continent. Hence, monitoring programs of climatic events like hurricanes in the region and environmental genomics surveys will add new insight to better understand possible barriers to the expansion of the moth into continental areas.