Propagule pressure and hunting pressure jointly determine genetic evolution in insular populations of a global frog invader

Islands are often considered to be more susceptible to biological invasions and to suffer greater impacts from invaders than mainland areas, and this difference is generally attributed to differences in species introductions, ecological factors or human activities between islands and mainland areas. Genetic variation, as a good estimate of evolutionary potential, can influence the invasion process and impacts of alien species. However, few studies have compared the genetic diversity of alien species between islands and a corresponding mainland. Here, we examined the genetic variation and differentiation in feral populations (30 sampled individuals/population) of a globally invasive species (the American bullfrog, Lithobates catesbeianus) that was extensively farmed on 14 islands in the Zhoushan Archipelago of China and in three nearby regions on the mainland. We quantified the relative importance of propagule pressure and hunting pressures on the genetic variation of bullfrog populations and found that insular populations have greater genetic variation than their mainland counterparts. Although genetic differentiation between the populations was observed, no evidence of recent bottlenecks or population expansion in any of the tested population was found. Our results suggest that the propagule pressures of bullfrogs escaping from farms, multiple releases and hunting pressure influence the genetic variation among bullfrog populations. These results might have important implications for understanding the establishment and evolution of alien species on islands and for the management of invasive species.

Islands are often considered to be more susceptible to biological invasions and to suffer greater negative impacts from invaders than the mainland 1,2 . This difference is generally attributed to the facts that islands often harbor more alien species, that islands have lower biological resistance (e.g., lower competition pressures or predation pressures) to invaders than the mainland areas; or that islands experience more human disturbances than the mainland areas [3][4][5][6][7][8][9] , and/or that islands exhibit a lower hunting pressure 10 . Genetic variation, a good estimate of population fitness and adaptive potential, may influence the invasion process and impacts of alien species 11,12 . Greater genetic diversity is hypothesized to promote rapid evolutionary responses to selection pressures in novel environments and may therefore facilitate invasion success and range expansion 13,14 . However, few studies have compared the genetic diversity and structures of alien species between islands and the mainland 15,16 .
Identifying the factors that influence the genetic variation of alien species on islands is not only important for understanding their successful establishment and evolution but is also crucial for the prevention and eradication of such species. The introduction history (e.g., propagule pressure and residence time), island environment and human activities (e.g., hunting pressure) are likely to influence the genetic variation of insular introduced

Number of bullfrogs raised in an enclosure, number of farms, residence time and island area.
We obtained information on the owners of the bullfrog farms from the Fisheries Bureau for the counties where the sampling sites were located. According to a regulation in China, a permit issued from the bureau is needed before a bullfrog farm can be operated. We then visited each owner to collect information on the number of bullfrogs raised on the farm and the periods during which bullfrogs were raised using a questionnaire survey 10 . The location of each farm was determined by GPS (geographic positioning system; Magellan eXplorist 210, Santa Clara, CA, USA). Bullfrog farms are generally located near a permanent still-water body or rice fields near permanent still-water bodies (less than 1 km from a water body) because they require a water supply. We counted the number of bullfrog farms within a plot (with a radius of 1 km centered at a sampling site) for each of the sampled sites in a region or island and summed the number of bullfrogs raised on the farms. For the sampled sites where a bullfrog farm was located, the residence time was considered as the time since the farm began to raise bullfrogs. We obtained information on residence time in the sampled sites without a bullfrog farm by using questionnaire surveys (Table S1, Supporting Information) 35 . We usually interviewed two or three resident individuals living around each of these sampled water bodies. The residence time was based on the time since these individuals first saw juvenile or adult bullfrogs or heard bullfrog calls. When several interviewers gave different answers for the residence time for a bullfrog-invaded site, the average value (year) of their answers was used. The residence time for all the sampled sites in a region or island is defined as the longest reported residence time.
Hunting pressure. We obtained information on the hunting pressure in invasive regions through the use of questionnaire surveys and field monitoring in 2011-2018 10 . For the questionnaire surveys, we interviewed three resident individuals living in the invasive region and asked "is there frog hunting activity at night. " The answers were usually "frequent hunting", "occasional hunting" or "no hunting", and the answers provided by the three residents in each invasive region were always consistent. For field monitoring, we recorded any human hunting activity in the invasive regions during three consecutive nights. Hunters usually used electric torches to find bullfrogs and caught them by hand or with a fishing net during the night. Therefore, we could find hunters by torchlight during the night. "Frequent hunting" was associated with higher hunting pressure. We classified hunting pressure as a dichotomous variable with the values "frequent hunting" and "no hunting". An invasive region was considered to "frequent hunting" if the questionnaire surveys yielded the answer "frequent hunting" or hunting activity was detected in two of the three nights monitored, and all other cases were defined as no (or occasional) hunting.
Sampling, DNA extraction and multilocus microsatellite genotyping. We randomly captured 30 postmetamorphic bullfrogs in invaded water bodies manually or with long-handled nets between May and October in 2011-2018. The distal third of the third toe on the right hind-foot from each frog was clipped. We then released the frog at its capture site. Each tissue sample from the toe clips was preserved separately in 95% ethanol in a 2-ml screw-cap microcentrifuge tube and stored at −20 °C in the laboratory. DNA was extracted following a published procedure 38,39 . In brief, the tips of bullfrogs' second hind toes (approximately 3 mg) were clipped and placed in a 2-ml centrifuge tube with 100 μl of lysis buffer containing 0.01 M NaCl, 0.1 M EDTA, 1 mg/ml proteinase K, 0.01 M Tris-HCl (pH 8.0) and 0.5% Nonidet P-40. The tube was vortexed for 1 min at ambient temperature, centrifuged to recover all material from the bottom of the tube, and incubated at 50 °C for 120 min and at 95 °C for 20 min. The tube was then centrifuged at 12,000 rpm for 3 min in a cold temperatures, and the extract was diluted to one tenth of its original concentration for polymerase chain reaction (PCR) amplification. The individuals were genotyped using nine nuclear microsatellite loci (BF1, BFD11 and GenBank accessions AY323928-AY323934). Sequence information for the primers and the amplification conditions used for these loci were based on previously published data 40,41 (Table S2, Supporting Information). All the primers were tagged with 5′-fluorescein bases (TAMRA, FAM or HEX). PCR reactions were performed as previously described 41,42 . The amplification conditions consisted of an initial denaturation at 94 °C for 3 min followed by 35 cycles of 10 s at 94 °C, 30 s at the annealing temperature (Table S2), and 30 s at 72 °C and a final 10-min extension at 72 °C. The PCR products were then separated by 2% agarose gel electrophoresis. Following amplification, the PCR products were resolved using an ABI PRISM 377 DNA Sequencer (Applied Biosystems) to resolve these PCR products. GENESCAN version 3.7 (Applied Biosystems) and GeneMarker version 1.71 (SoftGenetics) were used to score the microsatellite fragments.

Statistical analysis.
We applied MICRO-CHECKER 2.2.3 42 to quantify the scoring errors resulting from factors such as large allele dropout or stuttering and null alleles. GENEPOP version 4.0 43 was employed to test the linkage disequilibrium and Hardy-Weinberg equilibrium. We used Bonferroni corrections to test multiple comparisons 44 . Because the sampled bullfrogs within entire invasive regions or islands were collected, there might be some genetic structure among the sampled bullfrogs. We therefore examined the presence of genetic structure for the sampled bullfrogs within each sampled region or island using STRUCTURE version 2.3.3 45 based on the Bayesian clustering approach. This analysis revealed no genetic structure for the sampled bullfrogs for all sampled regions and islands (≥2 sampling sites) (Fig. S1, Supporting Information). Consequently, we treated the sampled bullfrogs within each region or island as a population.
We applied GenAlEx 6.5 46 to quantify the expected heterozygosity (He) and the observed heterozygosity (Ho), and the mean number of alleles (Na) for each region or island. Asymmetric migration rates were estimated using the MIGRATE-N 3.2.7 program 47 .
We examined the difference in the mean genetic diversity and allelic richness for bullfrogs between introduced regions on the mainland and the islands using a one-way ANOVA. We tested the differences in residence time, the number of bullfrogs raised and the number of farms between introduced regions on the mainland and the islands using a Wilcoxon rank sum test.
We determined the genetic divergence between regions or islands using Dest in the R package DEMEtics V0.8.0 48 , which is an unbiased estimator of genetic divergence for multiple loci 49 . The significance of the Dest results was assessed by bootstrapping with 1,000 permutations. In addition, we examined the relationship between the pairwise Dest and relative residence time using a Mantel test implemented in ARLEQUIN3.5 with 100,000 permutations 50 .
We tested for population bottlenecks using the program BOTTLENECK 51,52 . This software measures heterozygosity excess with respect to the mutation-drift equilibrium characteristic of population-size reductions for each population to infer recent population bottlenecks. We applied the one-tailed Wilcoxon signed rank test to determine the significance of heterozygosity excess. Two-phase mutation models were used: 90% single-step mutations and 10% multiple-step mutations with 10,000 replicates. To assess population expansion, we performed an intralocus k-test and interlocus g-test 53,54 using KGTESTS Excel Macro 53 . The intralocus k-test estimates the P value of k using the one-tailed binomial distribution. The significance of the g value for expanding populations was set at α = 0.05 54 .
Island area and the number of bullfrogs raised were log-transformed to improve linearity. Because several predictors were collinear (Fig. S2, Supporting Information), we used multimodel inference based on information theory 55 to evaluate the relative importance of the introduction history and hunting pressures on the expected heterozygosity. The full model was a multiple linear regression model with He as the response variable and with hunting pressure (frequent hunting = 1, no or occasional hunting = 0), residence time, number of bullfrogs raised and number of farms as predictors for island populations. All the models that included all possible combinations of the four (2 4 −1 = 15 models) variables were ranked. We compared the alternative models using the second-order Akaike information criterion (AICc) 55 and calculated the parameter estimates and their variances based on Akaike weights. Those models that were within two AICc units of the highest-ranked models (i.e., ΔAIC ≤2) 55 were reported. We conducted these analyses using the lme4 and MuMIn packages in R version 2.15.2 56 .

Results
Genetic diversity in insular and mainland populations. The number of sampled sites on an island or in a region ranged from one to four bodies of water and there was no significant difference between the island and mainland sampling points (Mann-Whitney U test: Z = 0.066, P = 0.948). The number of bullfrogs sampled in each site ranged from seven to 30 individuals in an island or from five to 30 individuals in a region (Table S3,  Supporting Information). MICRO-CHECKER did not detect the presence of null alleles or scoring errors. After a Bonferroni correction, we found no patterns of linkage disequilibrium among the loci and populations (P > 0.05). We also found that none of the populations or loci showed significant departure from Hardy-Weinberg equilibrium after Bonferroni corrections (P > 0.05). The loci were polymorphic within all populations.
The expected heterozygosity (He) and observed heterozygosity (Ho) for mainland populations ranged from 0.52 to 0.56 and 0.56 to 0.67 (Table 1), respectively, and for insular populations from 0.49 to 0.73 and 0.5 to 0.76, respectively. He differed from Ho across sampled regions or islands (Paired sample t-test, P < 0.001). The mean number of alleles (Na) and effective population size (θ) for mainland populations ranged from 3.78 to 4.33 and 0.0022 to 0.0145, respectively, and for island populations from 4.22 to 6.78 and 0.0022 to 0.0971, respectively. Overall, the genetic diversity (He) and allelic richness (Na) for island populations were higher than those for mainland populations (one-way ANOVA, F = 6.482, d.f. = 1, P = 0.022 for He; F = 6.604, P = 0.027 for Na).
Genetic differentiation, gene flow and genetic bottlenecks among insular and mainland populations. The pairwise population Dest differences ranged from 0.037 to 0.727 (significant α = 0.05 after Benjamini-Hochberg correction, P < 0.001 for all comparisons, Table S4, Supporting Information). These results indicated genetic differentiation between the population pairs.
Low-immigration-rate 95% CIs calculated by MIGRATE-N demonstrated that differences among populations did not differ significantly from zero, indicating that gene flow between these populations is likely lacking (Table S5, Supporting Information). We did not detect evidence for recent bottlenecks (one-tailed Wilcoxon signed rank test, P > 0.01) or population expansion in any of the populations (P > 0.05 Table S6, Supporting Information).

Correlation between factors and genetic diversity in insular populations. The number of bullfrogs
raised on farms and the number of farms around sampled water bodies ranged from 1,300 to 533,330 bullfrogs and one to four farms on each island and from 8,300 to 12,200 bullfrogs and only one farm in each region of the mainland ( island was between 19 and 23 years, and that for a region of the mainland was between 15 and 27 years (Table S7, Supporting Information). The number of bullfrogs raised on islands was higher than that in invaded regions of the mainland (Mann-Whitney U test, z = −2.269, P = 0.023), but there was no difference in the residence time between islands and regions of the mainland (z = −1.821, P = 0.069 for the number of farms; z = −1.204, P = 0.229 for the residence time).
The three most highly supported models (i.e., ΔAICc ≤2) for insular populations contained three variables: the number of bullfrogs raised, the residence time and hunting pressure ( Table 2). These predictors explained a proportion of the variation in the models (R 2 ) ranging from 0.8015 to 0.8691 for insular populations (Table 2). However, there was moderate model selection uncertainty across models (Wi = 0.206 to 0.304) ( Table 2).
The model averaging analysis also showed that the number of bullfrogs raised, the residence time and hunting pressure had high relative importance values (0.66-0.85) for expected heterozygosity in insular populations ( Table 3). The model-averaged 95% confidence intervals for these variables also excluded zero. The expected heterozygosity (positive parameter) increased with increases in the number of bullfrogs raised and in the residence time for insular populations, and the expected heterozygosity (negative parameter) decreased with increases in hunting pressure for insular populations. Figure 2 shows the relationship between expected heterozygosity and these predictors. The number of bullfrog farms had lower relative importance values (0.16) for insular populations, and the model-averaged 95% confidence intervals for these variables also included zero. Similar relationships between neutral genetic variation and ecological factors were found (Tables S8-10).

Discussion
We detected a statistically significant pattern in which greater genetic diversity and allelic richness in feral populations of bullfrogs on the islands of the Zhoushan Archipelago than on mainland China, although some island populations (such as Jintang, Daxie and Sijiao) have genetic diversity similar to mainland populations. A feral population would have greater genetic variation if there had previously been a larger number of bullfrogs raised on farms in the invaded region or island or if bullfrogs invaded early or invaded islands. These results suggested that both the introduction history and hunting pressure likely influenced the genetic diversity in bullfrog populations.
Several studies have compared the genetic variation in introduced populations between islands and the mainland. For example, Lade et al. 16 detected lower genetic variation in one insular population of Vulpes vulpes on Phillip Island compared to three populations on the Australian continent 16 . Similarly, higher genetic diversity was measured in eight populations of Neovison vison in Scotland (treated as a mainland) than one Outer Hebridean island population due to the admixture of introductions in the Scotland populations 15 . Therefore, the higher genetic diversity of bullfrog populations on the islands of the Zhoushan Archipelago than on mainland China is correlated to larger numbers of bullfrogs being raised on the islands than on the mainland and to lower hunting pressure in the islands. The propagule pressure of introductions plays an important role in determining the genetic variation of alien species 11,12,17,57,58 . The large number of bullfrogs raised may increase the number of bullfrogs escaping from the enclosure given that individual bullfrogs in the enclosures have the same chance of escaping, which would increase the input of more variation from escaped bullfrogs and therefore have positive effects on the genetic diversity of feral populations. In contrast to another study on invasive salmonids 20 , we found a positive effect of residence time on genetic variation in bullfrog populations, indicating that there might have been multiple releases at the sampled sites following the original introduction by escaped bullfrogs from farms. Bullfrog releases by both religious communities and individual persons are actively observed in China, where our sampled sites are located 33,36 . Such releases might have increased propagule pressure on the studied water bodies. In model averaging, hunting pressure might represent the main factor that influences the survival of bullfrogs escaping from farms on the islands. Human hunting pressures on bullfrogs were lower on the Zhoushan Archipelago than in the Beilun region of the mainland likely due to smaller human populations and lower demand for frogs as food throughout much of the archipelago 10 . Lower hunting pressures would be more likely to reduce the mortality of bullfrog founders from farms on islands and could therefore could result in higher genetic diversity in insular populations than in their mainland counterparts.
Studies have found that the admixture of multiple introductions from different source populations in a native range can increase the genetic variation in the invading populations 18,19 . The higher genetic diversity in insular bullfrog populations is unlikely due to such an admixture for three reasons. First, seawater might be a significant barrier to the dispersal of bullfrogs among islands due to the species' intolerance of seawater 30 . As a result, gene flow between island bullfrog populations might have been limited among the islands of the Zhoushan Archipelago. Second, a lack of genetic structure within all bullfrog populations indicates that multiple introductions from different source populations in the native range are unlikely to have occurred and contributed to genetic structure in each of these invasive populations. Third, this result is in line with a recent study based on the mitochondrial cytochrome b gene 34 , which suggests that a single introduction from a local population in the native range might have driven the limited haplotype diversity (only two haplotypes) in invading bullfrog populations in China 34,58 .
In contrast to studies that show population expansion without recent bottlenecks in invasive species 41,57,59 , we found no evidence for either recent bottlenecks or population expansion in either insular or mainland bullfrog populations (Table S6, Supporting Information), suggesting that the effective population size of bullfrog populations are stable after the establishment. Bullfrogs are characterized by their high fecundity (approximately 10,000 eggs/clutch) 60,61 , a potential that can result in high population densities and rapid population expansion. Feral bullfrog populations often meet adverse conditions in novel environments in the invaded range and cannot breed successfully in some years 28,35 , which might restrict this potential. Such demographic growth after introductions may result in the colonization of bullfrog populations with a minimal bottleneck effect or expansion effect on genetic diversity.
Both high and low levels of genetic differentiation have been observed in invasive species, including those with highly structured populations, such as Sturnus vulgaris 62 , Neovison vison 63 , Metrioptera roeselii 64 and Drosophila suzukii 65 , and populations with a low structure, such as Harmonia axyridis 66 and Galapaganus howdenae howdenae 57 ). We detected genetic differentiation in both insular population and mainland population pairs, suggesting the existence of genetic structuring exists among bullfrog populations. However, the differentiation (Dest) between pairs was not related to differences in residence time (Mantel test, r = 0.131, P = 0.088), suggesting that there were no differences in genetic differentiation between old and young bullfrog populations. Three factors  may be responsible for such differentiation. First, gene flow among mainland and insular populations might be very low (Table S5, Supporting Information), which would facilitate genetic differentiation among populations 67 . Similar to most amphibians 68,69 , bullfrogs are intolerant of seawater 30 , and sea barriers might isolate gene flow among insular populations. Second, large differences in the number of bullfrogs raised on farms (Table S7, Supporting Information) might lead to genetic differentiation among feral populations that originated from escaping founders. Finally, multiple releases might have an effect on the differentiation when such releases involve founders that are different from those of the first introduction 20 , which can increase or restore genetic diversity in invading populations. The results from this study might have important implications for understanding the establishment and evolution of alien species and for the management of invasive populations on islands. Enhanced genetic variation may increase the potential for adaptive evolution, which could lead to rapid evolution by natural selection in novel environments 12,19 . Although evidence exists for successful invasions that display reductions in genetic variation 70,71 , increasing numbers of studies show that higher genetic diversity facilitates invasion success 13,14,72,73 . Higher genetic variation in insular invasive populations than in their mainland counterparts would make islands more prone to biological invasions. The available data for bullfrogs are limited but consistent with this conjecture. Bullfrogs are more likely to invade from farms located on the Zhoushan islands than from farms in the nearby the mainland 10 . Future work is required to quantify the relative contributions of ecological factors and genetic factors to the differences in the invasion success and evolution of bullfrog populations between islands and the mainland.