Elucidating the invasion history of introduced bullfrogs in New Mexico using population genetic approaches

Sample collection

We obtained bullfrog tissue samples from museum loans and recent field collections. Field collections were completed under appropriate local and state permits: New Mexico Department of Game and Fish 3734, Rio Grande Nature Center 2023-010, and Bosque del Apache National Wildlife Refuge 2023-012. Samples archived at the Museum of Southwestern Biology (hereafter, MSB) at the University of New Mexico (UNM) were collected across the various river basins in northeastern and southwestern New Mexico between 2007–2022. Additional sampling was completed by us in the summer of 2023 within the Middle Rio Grande Basin and Lower Colorado River Basin. Specimen information is readily available for public access on the Arctos museum database (Table S1). Together the sampling consists of bullfrogs from Albuquerque (ABQ, n = 18), Socorro (SC, n = 3), Bosque del Apache (BDA, n = 4), Upper Gila River (UGR, n = 33), Cliff (CF, n = 4), Rodeo (RO, n = 14), and Mora River (MR, n = 1). We aimed to collect individuals from different ponds and along rivers approximately 500–1,000 m apart at sites to minimize collecting full sibling individuals. Individuals (various life stages: larval, juvenile, and adult) were caught by hand or dip net and transported to MSB where individuals were euthanized following approved University of New Mexico Institutional Animal Care and Use Committee Protocols (Protocols 20-201006-MC and 23-201375-MC). Specifically, we applied 20% benzocaine to the ventral side of the frog or tadpole until complete unresponsiveness was found using reflex tests, and then removed the heart and other tissues. Dissection tools were flame sterilized between frogs, and tissues were immediately flash frozen in liquid nitrogen and stored frozen until DNA extraction.

DNA extraction and sequencing

We extracted genetic DNA from liver, muscle, toe, and tail tissues using the E.Z.N.A. Tissue DNA Kit (Omega Bio-Tek, Norcross, GA, USA) following the manufacturer’s protocols. We used established polymerase chain reaction (PCR) protocols to amplify the mitochondrial gene, cytochrome b (cytb). Specifically, we amplified the 1,047-base pair (bp) gene region using the HERP328_cytb and HERP329_cytb forward and reverse primers (Yuan et al., 2016). This gene region has previously been used by other studies to investigate the invasion history of bullfrogs (Ficetola, Bonin & Miaud, 2008; Funk et al., 2011; Kamath, Sepulveda & Layhee, 2016) and to conduct phylogenetic analyses across the native bullfrog range (Austin, Lougheed & Boag, 2004).

We used 25 µL reaction volumes, each containing 30 ng of template DNA, one µL each of 10 mM forward and reverse primers, one µL of 25 mM MgCl2, one µL of 10 mM deoxynucleotide triphosphates (dNTPs), five µL of clear 5x GoTaq flexi buffer, 0.125 µL of GoTaq polymerase (Promega Corp., Madison, WI), and molecular-grade water. A negative control consisting of molecular grade water in place of template DNA was included with each batch of reactions. We ran reactions on a Bio-Rad T100 thermal cycler (Bio-Rad Laboratories, Hercules, CA) with the following thermal profile: an initial denaturation at 94 °C for 6 mins, 35 cycles of 94 °C for 20 secs, 52 °C for 30 secs, and 72 °C for 30 secs, and a final extension at 72 °C for 5 mins. The PCR products were visualized using a 1% agarose gel to confirm amplification of the gene region. Finally, PCR products were cleaned using Exo-SAP-IT (Bell, 2008) and sent for Sanger sequencing at Psomagen, Inc. (Rockville, MD).

Haplotype comparisons

We edited and assembled forward and reverse sequences within each individual using Geneious version 2025.03 (Biomatters Ltd, Aukland, NZ). Raw forward and reverse sequences were trimmed to remove primers. Consensus sequences were then aligned using the Muscle alignment software implemented in Geneious (Kearse et al., 2012) with a reference bullfrog sequence from NCBI GenBank (AF205089.1). We collapsed the alignment of New Mexico sequences to unique haplotypes using the FaBox web tool (Villesen, 2007). We then used the standard nucleotide NCBI Basic Local Alignment Search Tool (BLAST; Altschul et al., 1990) to compare each New Mexico haplotype to publicly available sequences and identify previously sequenced individuals with identical haplotypes. Where possible, we determined the locations of those sequences from referenced literature or museum catalog numbers. We considered identical haplotypes based on both full-length sequences and shorter sequences (408 bp) available from previous studies.

Population genetic diversity

We used the R package pegas (Paradis, 2010) to calculate genetic diversity metrics including haplotype diversity (H_d) and nucleotide diversity (π), using the functions hap.div() and nuc.div(), respectively. We used both full-length and shorter sequences (408 bp) for these analyses. Diversity metrics were calculated using default settings for each New Mexico site except Mora River (MR), which was excluded because it consisted of a single individual. We compiled H_d and π metrics previously reported in tables and supplemental materials from the following studies to make comparisons between invasive and native populations: Ficetola, Bonin & Miaud, 2008; Funk et al., 2011; Kamath, Sepulveda & Layhee, 2016; LaFond et al., 2022. We computed the same genetic diversity metrics (H_d and π) for the native populations sampled in Austin, Lougheed & Boag (2004) using Arlequin 3.5 (Excoffier & Lischer, 2010) and modified supplemental files from Kamath, Sepulveda & Layhee (2016). These data are summarized by site in Table S2. We ran a Shapiro normality test to assess whether the diversity metrics followed a normal distribution. Neither index met that assumption; therefore, we decided to use Mann–Whitney U tests (Hart, 2001) via the wilcox.test() function to determine whether H_d or π differed between invasive and native populations. Effect sizes were then calculated with the rstatix package (Kassambara, 2023).

Haplotype network and phylogenetic reconstruction

We assessed phylogenetic relationships of New Mexico haplotypes in comparison to native and other introduced haplotypes using sequences from GenBank. We selected one representative of each haplotype and included sequences that were longer than 600 bp (n = 35). The sequence alignment (full length = 949 bp) included both native (Austin, Lougheed & Boag, 2004; LaFond et al., 2022; MacGuigan et al., 2022) and invasive (Funk et al., 2011; Kamath, Sepulveda & Layhee, 2016; LaFond et al., 2022; Rodgers et al., 2023) bullfrog sequences.

We created a haplotype network in POPART v.1.7 using the Minimum Spanning Network (MSN) function for cytb sequences including one representative of each haplotype (haplotypes identical with New Mexico were combined) and included sequences longer than 600 bp (n = 27). We also used a maximum-likelihood (ML)-based approach to construct phylogenetic relationships among haplotypes using IQ-TREE version 2.3.6 (Nguyen et al., 2015). We used ModelFinder (Kalyaanamoorthy et al., 2017) to select an appropriate substitution model and UFBoot (Minh, Nguyen & von Haeseler, 2013) with 1,000 ultrafast bootstrap replicates to assess branch support. Nodes were considered well-supported with ML bootstrap support > 70% and we collapsed poorly supported nodes using TreeGraph 2 (Stöver & Müller, 2010). The outgroups used were Rana heckscheri and Rana Clamitans (GenBank Accessions AY083299 and AY083281).

Evaluating introduction history

To identify the minimum number of introductions of bullfrogs to New Mexico, we estimated pairwise ΦST between each New Mexico population pair with significance evaluated based on 10,000 permutations. Analyses were conducted in Arlequin 3.5 (Excoffier & Lischer, 2010) using both full length and shorter sequences (408 bp). We applied a sequential Bonferroni correction (Holm, 1979) to correct for multiple pairwise tests (n = 15). Significant ΦST values indicate that two populations are genetically distinct in relation to haplogroups, which are interpreted as evidence of independent introductions (Ficetola, Bonin & Miaud, 2008; Funk et al., 2011). In contrast, nonsignificant ΦST values indicate a lack of genetic distinctiveness and therefore no support for independent introductions of those populations from separate sources (Waples & Gaggiotti, 2006).

We assessed possible source populations for New Mexico bullfrog populations using evidence from the following two approaches. First, we inferred potential source populations based on the geographic distribution of identical haplotypes within the native and invasive ranges outside of New Mexico. Second, we used an analysis of molecular variance (AMOVA) to yield statistical evidence of a potential source population (Excoffier, Smouse & Quattro, 1992). This analysis estimates the amount of genetic variation between groups of populations and requires subdivision of the potential source population range into meaningful biological areas. We grouped the native range into four meaningful areas (Midwest, Overlap, East, and Northeast) according to the nested clade analysis previously published by Austin, Lougheed & Boag (2004) and following the same framework as previous studies of other invasive regions (Ficetola, Bonin & Miaud, 2008; Funk et al., 2011; Kamath, Sepulveda & Layhee, 2016) (Fig. 1C). The Midwest region in our analysis corresponds to the West region in prior studies. We then divided the other invasive populations into two regions, Southwest and Northwest, based on their geography.

We ran the AMOVA analyses in Arlequin 3.5 using modified input files from Kamath, Sepulveda & Layhee (2016). We used all available sequences from Austin, Lougheed & Boag (2004), Funk et al. (2011), Kamath, Sepulveda & Layhee (2016), and LaFond et al. (2022) that were assigned to haplotype and site. The East region included 20 sites with 185 sequences, the Northeast included 15 sites with 68 sequences, the Overlap included 16 sites with 98 sequences, the Midwest included six sites with 33 sequences, the Northwest included 15 sites with 506 sequences, and the Southwest included seven sites with 130 sequences. Each of these regions was compared to New Mexico with six sites and 76 sequences (excluding Mora, which had a single sequence). Separate analyses were conducted to assess the amount of molecular variance attributed to differences between each region (Midwest, Overlap, East, Northeast, Southwest, or Northwest) and New Mexico. Significant among-group variation suggests that a region is unlikely to be a source population. We used 10,000 permutations to determine significance.

Haplotype comparisons

We identified eight cytb haplotypes in bullfrogs across different sites in New Mexico (haplotypes NM1-NM8; Fig. 1A). Five haplotypes were restricted to a single site (NM4-NM7 in ABQ only; NM8 in CF only), while three haplotypes (NM1-NM3) were found at more than one site in NM. When we compared the full-length sequences of NM haplotypes to previous sequences, three were identical to haplotypes in the invasive populations in Montana and Wyoming (NM1 = MT1, NM3 = WY8, NM7 = MT2) and two haplotypes were identical to invasive bullfrogs sequenced in Arizona and California (NM2 = Rcat111, NM8 = Rcat162; Fig. 2). We also found two NM haplotypes that were identical to those from the native bullfrog range, one in the Overlap region in Mississippi (NM3 = RcatF = Rcat624) and one in the Northeast region in Maryland (NM6 = Rcat1415).

When the shorter, 408 bp sequences were compared, NM1 was identical to bullfrog sequences from the invasive range including California, Arizona, and Oregon; two from the native Overlap region in Ohio and Illinois; and one from Canada (Fig. 2). The short version of NM2 was identical to bullfrog sequences from Arizona and Oregon. The short NM3 sequence was identical to that of bullfrogs from the invasive range in Arizona, Oregon, and California, as well as within the Northeastern and Overlap native zones, including Mississippi, Tennessee, and Alabama. When only the short sequence is compared, NM4 was identical to NM1 and its matches in the invasive and native range. The short sequence of NM6 was identical to those in the Northeast, East, and Overlap native zones including states such as Alabama, Florida, Georgia, Virginia, and Maryland, and to sequences from the invasive range in the Northwest region in Oregon. The short version of NM7 was identical to sequences from the Midwest region of the native range, including Kansas, and from the invasive Northwest region, including Oregon and Montana. The short NM8 was identical to bullfrogs from the invasive range in California and Arizona. For both the full-length and short sequences, NM5 was not identical to any other haplotypes.

Population genetic diversity

The Albuquerque and Cliff populations exhibited the highest genetic diversity among the New Mexico populations (ABQ: H_d = 0.830, π = 0.009; CF: H_d = 0.833, π = 0.0034; Table 1; Fig. 1A). We observed the lowest genetic diversity in the Upper Gila River population (UGR: H_d = 0.435, n = 0.0009) and in the Rodeo (RO) site, which had a single haplotype sequenced from 14 individuals. Overall, the mean haplotype diversity was moderately high for New Mexico, but the mean nucleotide diversity was low (Table 1).

Mann–Whitney U-tests revealed significant differences in both haplotype diversity (W = 1511, P = 0.0004, r = 0.366) and nucleotide diversity (W = 1349.5, P = 0.0234, r = 0.236) among invasive and native bullfrog populations (Fig. 3). Both haplotype and nucleotide diversity were significantly lower in introduced populations compared to populations sampled from the native regions.

Haplotype network and phylogenetic relationships

The haplotype network showed the cytb haplotypes of New Mexico bullfrog populations were grouped within the two haplogroups separated by 11 mutations (Fig. 1B). Most haplotypes (NM1, NM2, NM3, NM4, NM7, NM8) were within Haplogroup 1, while only two haplotypes (NM5, NM6) were within Haplogroup 2. Most haplotypes were only separated from other closely related haplotypes by 1–3 mutations. The best-fit substitution model for the phylogenetic analysis based on ModelFinder was K3Pu+F+I, which includes three substitution types, empirical base frequencies, and a proportion of invariable sites. Similar to the haplotype network, the maximum likelihood analysis inferred two cytb haplogroups with moderately high support (bootstrap value = 82, Haplogroup 1; bootstrap value = 80, Haplogroup 2, Fig. S1). These two haplogroups are broadly distributed across the native range (Fig. 1C) and have a large area of overlap. There is apparent clustering of New Mexico haplotypes within the network and phylogeny that could indicate a common source for NM1 and NM4, and for NM5 and NM6 (Fig. 1B; Fig. S1).

Evaluating introduction history

Population genetic differentiation, measured as pairwise ΦST, between NM populations was low and not significant for most population pairs after sequential Bonferroni correction (Table 2). Two population pairs, ABQ-UGR and BDA-RO, were significantly differentiated after correction for multiple tests. However, based on chains of non-significant tests with other populations, the results suggest a minimum of one introduction to New Mexico.

The AMOVA revealed significant among-group variation between New Mexico bullfrog populations and three regions of the native range (Overlap: 39.8%, P= 0.0029, Northeast: 78.1%, P < 0.0001, East: 72.5%, P < 0.0001; Table 3), indicating these regions are unlikely source populations. Lower, marginally significant among-group variation was observed between New Mexico and the Midwest region within the native range (22%, P= 0.0503). The AMOVA also showed no among-group variation between New Mexico bullfrog populations and both regions in the invasive range (Northwest: 0%, P = 0.6411, Southwest: 0%, P = 0.5124).

Genetic diversity in introduced and native populations

In this initial characterization of the genetic structure of New Mexico bullfrog populations, we found intriguing variation in genetic diversity among sites. Specifically, we found lower genetic diversity at more isolated sites (two well-sampled sites, RO and UGR, had only one and two haplotypes, respectively), while genetic diversity was higher in the site with increased human activity (ABQ had six haplotypes). These haplotypes occur within two different cytb groups found across the native bullfrog range, and ABQ is the only site we sampled that included haplotypes from both haplogroups. We note that the phylogenetic support for these groups was moderate (bootstrap values 80–82), which may be attributed to the relatively short gene fragment used. This pattern could be attributed to increased opportunities in urbanized areas for invasive species to move to different populations and increase genetic diversity (Santana Marques et al., 2020). Six haplotypes found across New Mexico were identical to haplotypes dispersed across the native and other introduced ranges, and one more was identical when the shorter 408bp sequence was compared. This similarity is consistent with a recent introduction, presumably within the last century (e.g., NMDGF, 1921). The novel haplotypes we described could represent new mutations during that time or rare variants that would be revealed with more intense population sampling from other regions.

When comparing diversity metrics to populations from the native range, we found that invasive populations, including New Mexico, the Northwestern and Southwestern United States, and Europe, had significantly lower genetic diversity. This result corroborates previous studies that indicate successful bullfrog establishment is not constrained by lower genetic diversity (Ficetola, Bonin & Miaud, 2008; Funk et al., 2011; Bai et al., 2012; Kamath, Sepulveda & Layhee, 2016). Invasions are usually associated with genetic bottlenecks, which would be expected to limit the genetic diversity and adaptive potential of invasive populations (Ficetola, Bonin & Miaud, 2008; Estoup et al., 2016; Bors et al., 2019). Many studies have shown that multiple introductions from various regions could potentially introduce enough genetic diversity for an invasive species to thrive (Dlugosch & Parker, 2008; Winkler et al., 2019). For example, multiple introductions of Cuban brown anoles to Florida introduced higher genetic diversity compared to native counterparts that could be contributing to the success of invasive species (Kolbe et al., 2004). In bullfrogs, studies of other introduced populations in Montana (Kamath, Sepulveda & Layhee, 2016) and Europe (Ficetola, Bonin & Miaud, 2008) concluded there were multiple introductions, but with low genetic diversity.

In contrast, our data is consistent with one introduction of bullfrogs to New Mexico. Interestingly, New Mexico has several haplotypes that are identical to the native range, suggesting that a single invasive founding population could be from multiple sources that may have yielded enough genetic variation for the success and dispersal of bullfrogs across the state (Lavergne & Molofsky, 2007; Schrieber & Lachmuth, 2017). Many other factors could contribute to the success of invasive bullfrog populations despite low genetic variation. Their opportunistic predation habits, high mobility, and larger body size may confer advantages compared to native species (Snow & Witmer, 2010; Jancowski & Orchard, 2013). Furthermore, faster maturation time and presumed higher fecundity in introduced populations compared to native populations could provide a competitive advantage and maintain high population sizes (Urbina et al., 2020). Phenotypic plasticity may also be contributing to the success of bullfrogs in New Mexico enabling survival in new environmental conditions (Hagenblad et al., 2015). However, we caution that the inference of a single introduction is based on a single gene and ΦST estimates with large variances that may not be supported in future studies with more comprehensive genetic data, including additional genetic markers (Keller & Taylor, 2008; Sturk-Andreaggi et al., 2022). Sampling intensity has also been connected to the accuracy of ΦST estimates, therefore affecting the approximation of the number of introductions and gene flow (Dlugosch & Parker, 2008).

Identifying potential sources and connectivity of introduced populations

Bullfrogs have been linked to various source populations and modes of introduction in previous studies. Many populations in China were attributed to frogs escaping bullfrog farms and food markets that were first imported from Cuba and Japan (Xuan, Yiming & McGarrity, 2010; Bai et al., 2012). In Montana and Wyoming, the native source population inferred was the Midwest, which may also be linked to bullfrog farming (Kamath, Sepulveda & Layhee, 2016), but other vectors such as the exotic pet trade and bait market for fishing, were also considered (Strecker, Campbell & Olden, 2011; Sepulveda et al., 2015). In Oregon, the native source population inferred was the Overlap region, which was attributed to a possible connection with bullfrog farms (Funk et al., 2011). In our study, we did not find conclusive evidence for a native source population, though the native group most similar to New Mexico, based on the AMOVA, was the Midwest, as in Kamath, Sepulveda & Layhee (2016). In addition, the AMOVA comparison found no significant differences between New Mexico and the invasive Northwest and Southwest regions, indicating a possible link between these populations. This finding could indicate a common source population for these introductions or a recent movement of populations from one invasive area to another.

While bullfrog farms may have contributed to introductions in the Northwest (Funk et al., 2011; Kamath, Sepulveda & Layhee, 2016), they are not a current threat for continued introduction because of their closures across North America (Helfrich, Neves & Parkhurst, 2009), nor are they likely to explain introductions to the Southwest. Historically, the Southwest did not have established bullfrog farms, but rather a record of intentional release and cultivation of bullfrogs for recreational hunting and fish stockings. In Arizona, bullfrogs were introduced for recreational hunting in the 1950s near Tucson and were stocked in rivers and lakes within a 100-mile radius (Brennan & Holycross, 2006). In California, bullfrogs were first introduced in 1912 for recreational hunting by a Louisiana dealer to the Amargosa River and were then moved to other aquatic environments in the Mojave Desert (Washington Department of Fish and Wildlife, 1995). There was also a separate introduction in California that involved a different population of bullfrogs from a San Francisco frog merchant in 1914, which was introduced into a nearby reservoir (Washington Department of Fish and Wildlife, 1995). In New Mexico, there was a similar history where the government intentionally released bullfrogs to the Rio Grande from a single population of bullfrogs from the Kansas State Fish Hatchery in 1921 (NMDGF, 1921).

The genetic similarities we identified among populations in New Mexico and the Southwest region, including California and Arizona, suggest a link between these populations. Arizona is geographically proximate to New Mexico and is connected through the Colorado River basin, but migration of bullfrogs over longer distances is less likely (Sepulveda & Layhee, 2015). Human-mediated dispersal is therefore a more probable explanation for the introduction and spread of bullfrogs across New Mexico. The identical haplotypes and lack of among-group variation between New Mexico and the Southwest region suggest another invasive population could have been a direct source. While our data is consistent with one single introduction to New Mexico, the current bullfrog range in the state far exceeds the Rio Grande, where bullfrogs were initially introduced, and the northeastern edge of the state, which may be part of the native range. Secondary spread after initial introduction is likely (Vander Zanden & Olden, 2008), and there is now greater concern to investigate secondary spread to improve management of New Mexico populations. Management strategies such as implementing laws to stop the import of invasive species into a region (Reed et al., 2023) and eradicating current populations (Hossack et al., 2023) have been successful in halting continued introduction and secondary spread in some states in the U.S. Current eradication efforts in New Mexico are focused on specific areas of overlap with vulnerable native species, but more concerted removal efforts and monitoring are critical to understanding the impacts of invasive bullfrogs across the state.

Limitations and future directions

This study applied genetic approaches to provide a better understanding of the invasion history of an immensely successful invader. We provide key insights about the genetic composition of invasive bullfrog populations in the western U.S., including the first accounts of the invasion history in New Mexico. However, we recognize the limitations of using a single mtDNA gene region to uncover fine-scale dispersal patterns and gene flow (Hurst & Jiggins, 2005). Given the lower effective population size of mtDNA compared to nuclear DNA, this can affect inferences about genetic bottlenecks in population genetic studies (Ferreira & Rodriguez, 2024). Using additional genetic markers, including genome-wide loci, would improve resolution and enable the study of genetic adaptation (Sturk-Andreaggi et al., 2022). We also acknowledge that our sampling does not encompass all bullfrog populations, and varying sample sizes, especially low sample sizes in close geographic proximity, can affect measures of genetic diversity and haplotype estimates (Rosenberger et al., 2021). Increasing sample sizes and expanding sampling to include other populations would be helpful in future work. Our results provide insights into the connectivity of invasive populations and the importance of considering secondary spread to inform management strategies. Understanding the history of invasive species is critical in managing existing populations, and future work should continue monitoring population growth to avoid further declines of native species.