Genetic structure of American bullfrog populations in Brazil

Non-native species are a major problem affecting numerous biomes around the globe. Information on their population genetics is crucial for understanding their invasion history and dynamics. We evaluated the population structure of the non-native American bullfrog, Aquarana catesbeiana, in Brazil on the basis of 324 samples collected from feral and captive groups at 38 sites in seven of the nine states where feral populations occur. We genotyped all samples using previously developed, highly polymorphic microsatellite loci and performed a discriminant analysis of principal components together with Jost’s D index to quantify pairwise differentiation between populations. We then amplified 1,047 base pairs of the mitochondrial cytochrome b (cytb) gene from the most divergent samples from each genetic population and calculated their pairwise differences. Both the microsatellite and cytb data indicated that bullfrogs comprise two populations. Population grouping 1 is widespread and possesses two cytb haplotypes. Population grouping 2 is restricted to only one state and possesses only one of the haplotypes from Population grouping 1. We show that there were two imports of bullfrogs to Brazil and that there is low genetic exchange between population groupings. Also, we find that there is no genetic divergence among feral and captive populations suggesting continuous releases. The limited genetic variability present in the country is associated to the small number of introductions and founders. Feral bullfrogs are highly associated to leaks from farms, and control measures should focus on preventing escapes using other resources than genetics, as feral and captive populations do not differ.

www.nature.com/scientificreports/ population dynamics 35 , but a study focused on its known and potential impacts indicates that the presence of bullfrogs can cause changes in activity and habitat use preferences of native species 36 .
Although there is very limited knowledge about the possible impacts of bullfrogs over native species in Brazil, this species is widespread in the south and southeast regions of the country, with limited distribution in the northeast and north 37,38 . There are over 150 known breeding facilities in the country, that produce around 400 tons/year, moving around 1.9 million USD 38 . According to published accounts, bullfrogs were first imported to Brazil to produce meat for human consumption in 1935 by a Canadian technician named Tom Cyrril Harrison who brought either 300 individuals 39 or pairs 40 from an unknown locality to a Brazilian government breeding facility in the city of Rio de Janeiro 41,42 . Hundreds of tadpoles were subsequently sent to new breeding facilities around the country in an agricultural program encouraged by federal and state governments 41,42 . The state-based characteristics of this governmental program facilitated the appearance of feral bullfrog populations in different political regions of the country, as today feral bullfrogs are widespread in the south and southeast regions of Brazil, with additional localities in the central-west, northeast, and north 41 . Escape and release of bullfrogs from farms are spatially correlated with feral bullfrog populations [43][44][45] , being the most likely cause considering the disconnected distribution of these populations in Brazil, although the dispersive potential of bullfrogs is not well documented. Although advances in technology and increased market demand have enabled bullfrog farming to expand greatly over the last 10-15 years 46 , many of the earlier frog farms failed in the 1990s 38 , resulting in massive releases of bullfrogs when businesses closed (Valdir Alves, owner of Pedrinhas bullfrog farm, personal communication). In the 1970s, an additional 20 pairs of adult bullfrogs were imported to Brazil 28 from the University of Michigan to São José do Rio Preto in São Paulo state by Luiz Dino Vizotto (C. M. Ferreira, personal communication). Little information about these individuals is available, and the outcome of this introduction is not confirmed nor discussed in the literature. The same is true for small-scale occasional introductions supported by individual breeders, without any official registers, such as a supposed bullfrog batch introduced from Mexico in the early 2000's 38 .
If all the introduction events were successful, different lineages in the country might be going through an introgressive hybridization process due to the common practice among farmers of purchasing breeding stock from multiple farms in different states with the goal of increasing the genetic diversity of the frogs bred in their facilities (Romar's bullfrog farm owner, personal communication) 38 . If only one event was successful, followed by strong selective pressures of breeding facilities, the genetic diversity of bullfrog populations in Brazil should be quite limited. The limited information about introduction history of the bullfrog in Brazil raises the possibility of different genetic scenarios. At least one population was successfully introduced in the country, and other lineages could have also contributed genetically to the introduced populations.
Information on introduced bullfrog population genetics is crucial to understand invasion history, structure and dynamics of gene pool exchange between populations 47,48 , and can be helpful to develop management and control programs 49,50 . Efforts to understand the genetic structure of introduced bullfrog populations have been undertaken in several regions, such as Europe 51 , China 43 and western USA 48,52 . Unlike the invasion history in Brazil, Europe had several events of introduction from different origins 51 , and the western region of the USA is continually flooded with individuals from the native range 52 . China also had different events of introduction, although little genetic diversity was found 43 . These examples can serve as a baseline for comparison in new genetic research, and give indications about the invasiveness of a non-native population.
Because multiple introduction events can create introduced populations that present more diverse gene pools 52 , we expect introduced bullfrog population in Brazil to be less diverse than other non-native populations previously studied in different regions of the globe, as fewer introduction events were reported in Brazil in relation to other areas. Considering that more unreported introduction events could have happened, our main goal in this study was to assess the genetic structure of Brazilian bullfrog populations and possible gene flow between them. To achieve our goal, we focused on four main questions: (1) Are captive and feral populations genetically different? (2) Does bullfrog genetic structure corroborate the hypothesis that there were two or more introduction events in Brazil? (3) Is there gene flow between populations? (4) Do Brazilian bullfrog populations present lower genetic diversity than non-native populations from other regions?

Materials and methods
We obtained 324 bullfrog skin, liver, and/or muscle tissue samples, of which 128 were purchased from 11 farms and 194 were collected from feral specimens or sampled from museum or private collections ( Fig. 1; Table SI). Specimen collection was conducted in accordance to the guidelines and authorized by the Instituto Chico Mendes de Conservação da Biodiversidade (SISBio 56772-1). The sampling methodology also followed the ARRIVE guidelines 53 . Farm animals were not sacrificed for this research. Tissues were stored in 99% ethanol and kept at − 20 °C until DNA extraction, which we performed with the DNeasy Blood and Tissue extraction kit (Qiagen, Valencia, CA, USA) following the manufacture's guidelines.

Microsatellite analyses.
We amplified seven nuclear microsatellite loci using the library developed by Austin et al. 54 . Multiplex PCR was performed on a Veriti™ thermal cycler (Applied Biosystems) with a thermal profile consisting of 95 °C for 7 min followed by 10  www.nature.com/scientificreports/ provider to avoid bias due to calibration differences. We scored results using Gene Marker v. 2.6.3 (SoftGenetics) and tested for the presence of null alleles, allele dropout, and stuttering using Microchecker 55 . We also genotyped 91 samples (24.1%) a second time to evaluate the percentage of homozygotes that were actually null alleles. We included control samples in each procedure. Given that published sources report a total of 340 or 640 bullfrogs being imported to Brazil (300-600 in 1935 and another 40 in the 1970s 28,39,40 , we anticipated two possible scenarios: (1) extremely low genetic divergence and diversity among bullfrog populations in Brazil. This could have been caused by limited genetic diversity in the founding event of Brazilian populations, together with strong selective pressure in breeding facilities; 2. Increased genetic diversity caused by founding populations that were prevenient from different native populations. Considering that each scenario would require different analytical methods, we chose for assessing the genetic structure using two different analytical toolkits. A flowchart with the analytical decision-making is available in Fig. 2.
For scenario 1, we performed a discriminant analysis of principal components (DAPC) 56 in the adegenet R package 57 to assess current genetic structure. This analysis includes all bullfrog samples in the country grouped by state, in recognition of the state-based government programs regulating bullfrog farming. To determine if there is any genetic differentiation between groups defined by feral or captive origin from the same locations-which could be used to inform law enforcement in case of releases-we performed one DAPC only for states in which both groups were well sampled (Rio de Janeiro and São Paulo). www.nature.com/scientificreports/ We used Jost's D index 58 to quantify pairwise genetic differentiation among state populations, which was calculated with the mmod R package 59 . Next, we tested the significance of differentiation between pairs of populations using the DEMEtics R package 60 with 10,000 permutations. We applied the Benjamini-Hochberg Procedure 61 to control for false discovery rates (FDR) and avoid type I error 62 using the p.adjust function from stats package 63 . We used the gstudio package 64 to calculate the following diversity indices for each genetic population: expected (He) and observed (Ho) heterozygosity, number of alleles (A) per locus, effective number of alleles (Ae) per locus, and size-corrected Wright's inbreeding coefficient (Fis). Through a permutation procedure with 100 batches of 1,000 iterations, we checked for Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) in each state population by performing a probability test in Genepop ver. 3.4 65 . To avoid type I error, P values were corrected following the Benjamini-Hochberg Procedure.
Considering scenario two, we performed a principal coordinate analysis (PCoA) to visualize genetic differentiation using the R package Adegenet 57 . We estimated admixture using STRU CTU RE software 66 . This analysis was performed using nine runs for every K (from K = 1 to K = 9) using state populations in the locprior model and a burn-in period of 10.000 steps and 20.000 Markov chain Monte Carlo (MCMC) repetitions. The best fitted K number of populations was calculated using STRU CTU RE HARVESTER 67 , and interpreted through the Evanno method 68 .
We also calculated the differentiation among and within state populations and feral and captive populations through an analysis of molecular variance (AMOVA). We used Arlequin software 69 to perform these analyses.

Cytochrome b analyses.
On the basis of the microsatellite results, we selected the 18 most divergent samples from each genetic population (Table 1) 72 for sequence editing and contig formation of the cytb sequences based on the chromatograms obtained from the automated sequencer. All samples were sequenced in both directions to check for potential errors. The sequences were aligned with the MAFFT 73 plugin in Geneious 10.2.3 (Biomatters) with the G-INS-I strategy. We trimmed the sequence alignment to a length of 937 bp. We used Arlequin 3.5 69 for most genetic analysis. We calculated the haplotype diversity (Hd) and nucleotide diversity (π) indices and report values ± SD for each genetic population. We also calculated the pairwise differences between the introduced populations using the θST index and evaluated its significance by performing 10,000 permutations. Using the p.adjust function from stats package 63 , we performed the Benjamini-Hochberg Procedure 61 to control for FDR.

Results
The presence of null alleles was indicated for all loci, whereas allele dropout and stuttering measures were not significant. Genotyping error was estimated to be less than 2%. The state-based total-data DAPC (Fig. 3) showed extensive overlap between most sampling units (defined as Population grouping 1) except individuals from Minas Gerais (Population grouping 2).
Pairwise multilocus Jost's D between Population grouping 1 and Population grouping 2 was significant even after FDR correction (Jost's D = 0.25638, P = 0.001). Most pairwise comparisons indicated significant LD (Table SII) and deviation from HWE (Table 2), even after FDR correction. The mean number of alleles per locus  www.nature.com/scientificreports/ (A) within each population was 4.28-9.57, while the effective number of alleles was only 2.86-4.05 due to the strong dominance of one or two alleles in most loci (dominant alleles accounted for 20-50% of the total at each locus). The observed heterozygosity did not match expected heterozygosity and both population groupings showed positive inbreeding values ( Table 2). The PCoA analysis did not present any patterns of genetic differentiation (Fig. 4. Conversely, the clustering analysis performed by STRU CTU RE found K = 3 to the best fit for genetic populations (Fig. 5), supported by the Evanno method (Supplemental material, Fig. SII). The state-based AMOVA (Table 3) indicated 77.49% of the genetic variation should be within individuals, and only 2.98% among groups (states). Similar results were found for the AMOVA between feral and captive specimens (84.39% and 0.07%, respectively). The DAPC that tested the differentiation between all feral and captive groups (Supplemental material, Fig. S1) presented almost complete overlap of sample units, supporting the results found in the AMOVA and leading to the conclusion that feral and captive animals are not significantly different.
We identified only two cytb haplotypes in Population grouping 1 haplotypes (Haplotypes A and B) that differ from each other in 14 base pairs (π = 0.005324 ± 0.003); Haplotype A was more frequent than Haplotype B, occurring in 80% of specimens (Hd = 0.3556 ± 0.159). Population grouping 2 contains Haplotype B exclusively.  www.nature.com/scientificreports/ The population differentiation test between the two population groupings was significant after FDR correction (pairwise θST = 0.75535, P = 0.00147).

Discussion
All bullfrogs in Brazil descend from a relatively small number of founders introduced at most 85 years ago, in at least two different events 28,39,40 . The introduction history of this species resulted in high rates of homozygosity, high inbreeding coefficient values, low allelic richness, and consequent overestimations of LD and HW disequilibrium. These characteristics are predictable, as they are associated with the founder effect and inbreeding in these populations following their importation to Brazil and subsequent intensive reproduction and presumed artificial selection. Populations with low genetic variation, as expected in cases of non-native populations with few founders or populations that have undergone major reductions 73 , present special conditions that require appropriate methods. As such, we used DAPC to analyze genetic structure because it makes no assumptions about HWE, LD, or gene flow, like more commonly used analyses such as STRU CTU RE 66 .
Similarly, although null alleles are known to inflate measures of genetic differentiation and create false homozygotes 74,75 , the common practice of discarding loci inferred to possess null alleles [e.g. 50,77,78 is contraindicated in cases such as the introduced bullfrog in Brazil. Specifically, Microchecker indicates the presence of null alleles on the basis of excess homozygotes being evenly distributed across homozygote classes 55 , which is precisely the situation expected to occur in recent invasions subjected to high levels of inbreeding [e.g., 9 . The indication of null alleles in all genes reinforces this point, as all genes underwent the same biological process and none of them are in HWE. Nevertheless, the double genotyping of almost a quarter of all microsatellite data showed that the presence of null alleles is not significant. Thus, we included all loci in our analyses.
Considering that other non-reported introductions could have happened, increasing the gene pool and contradicting our first methodological toolkit assumptions, we reanalyzed our dataset with a more traditional approach. This included PCoA and STRU CTU RE, to visualize genetic differentiation and admixture, respectively; and AMOVA, to teste differentiation among populations.
Both methodological toolkits presented similar results, with low differentiation among population groupings. Our state-based analysis of microsatellite data corroborates published accounts that there have been only two importations of bullfrogs to Brazil, a well-known importation of 300-600 individuals in 1935 and a lesser-known importation of 20 pairs in the 1970s 28 . Specifically, although we observed that samples from most areas overlap extensively by occupying the same multivariate space in the DAPC (Population grouping 1), samples from Minas Gerais did not (Population grouping 2). Similarly, the best number of Ks indicated for the STRU CTU RE analysis was three. Although there was a difference in the interpreted number of populations between the two used approaches, it could have happened due to biases in both of them. The DAPC approach is subjected to a certain amount of subjectiveness, as different populations are interpreted by their distribution in the multivariate space. In this case, the population from Parana, for instance, has a distribution that does not completely overlap with the rest of samples named here as Population grouping 1. This could be the third K represented in the STRU CTU RE approach. At the same time, as mentioned before, the assumptions about HWE, LD, or gene flow in the STRU CTU RE analysis could also be responsible for inflating the real number of Ks.
The results of the cytb analysis corroborate the DAPC results and further suggest that either founders of Population grouping 1 contained both mtDNA haplotypes, or some degree of introgression of the haplotype from Minas Gerais (Population grouping 2) into Population grouping 1 has occurred but not the other way around. This might be correlated with the decrease of frog farmers in the state of Minas Gerais, previously known as one of the biggest producers in the country. Farmers from other states, like São Paulo where some haplotypes from population grouping 2 were found, made efforts to enhance the genetic diversity of their breeding population by importing animals from other states while breeders from Minas Gerais closed or reduced their facilities, preventing Population grouping 2 from receiving migrants from Population grouping 1. We therefore conclude that populations from the same genetic units are probably exchanging migrants, since there is no signal of genetic differentiation between them, what agrees with the information that farmers keep exchanging animals. This practice is unlikely to be effective, given the low genetic variation of bullfrogs in Brazil, and might also have contributed to preventing the genetic differentiation of introduced populations. That is, by exchanging "migrants" with each other and replenishing feral populations with constant leaks 44 , breeding facilities effectively maintain bullfrog populations across the country in an infinite island model. This process was documented in a recent work 32 , where authors investigated the historical and present bullfrog trade among Brazilian states. In their work, the states of São Paulo, Rio de Janeiro and Minas Gerais seem to be central to this market. They did not contrast the trade of living bullfrogs from the meat trade, what might explain why Minas Gerais still has differentiated genetic traits even though it is deeply entangled in the commercial net. www.nature.com/scientificreports/ The lack of differentiation among populations prevents the use of genetic resources to diagnose and control leaks from breeding facilities. As shown in both DAPC and AMOVA, captive and feral specimens do not differentiate. This result can be explained not only by biological processes, like the low number of founders and intense inbreeding following importation, but also by the common practice among farm managers of purchasing breeding stock from multiple farms in different states with the goal of increasing the genetic diversity among the frogs bred in their facilities.
The practices of bullfrog farmers in Brazil tend to generate an undifferentiated population in the country over time. The same process seems to have occurred in other countries where populations were analyzed with the same genetic markers used in this study. Similar patterns were observed in China 43 -where also only two haplotypes were found among over 500 samples-and in some degree in Europe (5 haplotypes among nearly 400 samples from 8 countries) 47 , although European countries have more populations and more genetic diversity, due to the relatively high number of introduction events that happened in the continent (at least 25) 51 . The process of transforming introduced populations into a single population due to exchange of migrants may obscure population and invasion genetics inference, as the genetic signal is lost making it nearly impossible to clarify the history of events that followed the introductions. This also brings consequences to management efforts as important information like the origin and the dynamics of the introduced populations are let largely unanswered.
The number of introduction events seems to be directly related to the degree of genetic diversity in the region, and the bullfrog populations in Brazil exhibit the lowest number of mtDNA haplotypes of all studied non-native populations of this species examined so far 43,47,52 . Although the genetic diversity of bullfrogs in Brazil seems to be low, especially when compared to populations in its native range (42 haplotypes were found in the United States 79 ), farming does not seem to suffer any impacts.
The exceptionally reduced genetic variation found in the population from Minas Gerais might explain the observation of morphological anomalies reported by Ferrante et al. 35 . Conversely, we observed a small percentage of anomalous animals when visiting farms, and farmers did not report any negative effects in their business. Research focused on other invasive amphibian species, such as the cane toad (Rhinella marina + R. horribilis species complex), show that reduction in genetic diversity after introduction does not seem to affect ecologically relevant traits 80 . This matter should be accounted when evaluating the potential of bullfrog invasiveness. Feral bullfrogs are highly associated to leaks from farms, and control measures should focus on preventing escapes using other resources than genetics, as feral and captive populations do not differ. The limited genetic variability present in the country is associated to the small number of introductions and founders.

Data availability
Available as Supplementary Data 1 and 2. www.nature.com/scientificreports/