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

TL;DR

This study uses genetic analysis to trace the invasion history of bullfrogs in New Mexico, revealing patterns of genetic diversity and likely sources of introduction.

Contribution

The study provides new insights into the genetic structure and invasion pathways of bullfrogs in New Mexico using population genetic approaches.

Findings

Bullfrog populations in New Mexico show moderate genetic diversity with higher diversity near human activity.

Introduced populations have significantly lower genetic diversity than native populations, suggesting recent colonization.

Analysis suggests bullfrogs in New Mexico likely originated from a single source population.

Abstract

The American Bullfrog, Rana (Aquarana) catesbeiana, is an extremely successful invader that has spread globally in the last century, impacting vulnerable ecosystems. In the western U.S., bullfrogs were intentionally introduced in the early 1900s, but little is known about their subsequent colonization into the state of New Mexico. We evaluated a single mitochondrial gene region, cytochrome b, using population genetic approaches to investigate the invasion history of bullfrogs across their introduced range. Specifically, our objectives were to (1) assess the level of genetic diversity and identify haplotypes within bullfrog populations in New Mexico, (2) compare the genetic diversity of native and invasive bullfrog populations, (3) infer the number of introductions into New Mexico, and (4) identify potential native source populations. Using haplotype and nucleotide diversity estimates, we found moderate genetic variation within New Mexico (Hd = 0.648, π = 0.0036) with higher diversity at sites with increased human activity. However, there was significantly lower genetic diversity in introduced populations compared to native populations, consistent with expectations of recent colonization. Based on haplotype diversity estimates and BLAST results, we found a total of eight haplotypes across New Mexico, of which six haplotypes were found across the native and other introduced ranges. Pairwise ΦST revealed minimal differentiation between New Mexico sites, consistent with introduction from a single source population. Lastly, the analysis of molecular variance (AMOVA) conveyed that bullfrogs from the Northwest (Northwest: 0%, P = 0.6411) and Southwest (Southwest: 0%, P = 0.5124) invasive regions showed no significant differences compared to New Mexico populations, suggesting either recent connectivity or similar origins. This study reinforces the importance of managing the movement of invasive species and demonstrates how evaluating the genetic composition of an invasive species can reveal key points of its invasion history.

Figures

Haplotype and haplogroup distribution of bullfrogs sampled across the United States.(A) Cytochrome b haplotype proportions by sampling sites across New Mexico. Different haplotypes are represented by colors. (B) Haplotype network for cytochrome b haplotypes across the native and invasive range. Color indicates region where haplotype is found. Size of pie charts is proportional to haplotype freqency. Dashes on lines between haplotypes represents number of mutations. (C) Distributions of estimated bullfrog haplogroups according to phylogenetic analyses and haplotype network across invasive and native ranges. Silhouette of bullfrog was obtained from phylopic.org. U.S. County-level ranges were obtained from The International Union for Conservation of Nature’s Red List (IUCN, 2022).

Distribution of cytochrome b haplotypes found in New Mexico bullfrogs across native and invasive ranges.Colors represent identical haplotypes to invasive populations found in New Mexico. Other haplotypes not identical to those in New Mexico are shown in grey. Both invasive and native ranges are shown on the map as shaded regions and were originally obtained from The International Union for Conservation of Nature’s Red List (IUCN, 2022) and adjusted according to new and existing museum records. The native range was divided into four meaningful biological areas, visualized by the dashed lines, based on the nested clade analysis published in Austin, Lougheed & Boag (2004). Additional sequences from Funk et al. (2010), Kamath, Sepulveda & Layhee, 2016, Lafond et al. (2022), MacGuigan et al. (2022), Rodgers et al. (2023) were assigned to these regions according to geography (Northeast, East, Overlap, Midwest). The invasive range was split into two regions (Northwest and Southwest) based on geographical boundaries.

Distribution of genetic diversity metrics for native and invasive populations.Distribution of haplotype diversity (A) and nucleotide diversity (B) of native and invasive populations. Mann–Whitney U test statistics (W) and significance values (P) are provided for each comparison.

Full text

Introduction

The American Bullfrog, Rana (Aquarana) catesbeiana, is a successful invader that was initially introduced to non-native regions both intentionally, for recreational hunting as a food source (Brennan & Holycross, 2006), and unintentionally, through release from bullfrog aquaculture (Jennings & Hayes, 1985) and fish stockings (McAuliffe, 1978). Since then, American Bullfrogs (hereafter, bullfrogs) have expanded from their native range in the eastern United States of America (hereafter, U.S.) to the western U.S. and to more than 40 countries (Lowe et al., 2019). Listed in the IUCN’s Top 100 Most-Invasive Species, bullfrogs have been linked to species’ declines across the globe because of disease spread (Urbina et al., 2018), direct predation (Jancowski & Orchard, 2013), and competition for resources with native species (Snow & Witmer, 2010). This impact is especially prominent where native populations are already under high environmental stress due to a changing climate, including limited water availability and rising temperatures (Prăvălie, 2016).

The Southwestern U.S. is a climate-stressed region in which the interactions between invasive and native species have been relatively understudied (Archer & Predick, 2008; Dettinger, Udall & Georgakakos, 2015; Dueñas et al., 2018). Bullfrogs exacerbate environmental stressors by outcompeting native amphibians for limited aquatic habitat and preying on vulnerable populations (Hayes & Jennings, 1986; Hossack et al., 2023). For example, bullfrog establishment in the Southwestern U.S. has been directly linked to declines of the federally endangered Chiricahua Leopard Frog (Rana chiricahuensis) and Sonoran Tiger Salamander (Ambystoma mavortium stebbinsi) (United States Fish and Wildlife Service, 2002; United States Fish and Wildlife Service, 2007). In extreme environments such as the deserts of the American Southwest, where native species are already pushed to their physiological limits, the presence of a dominant invasive predator may further accelerate biodiversity loss (Hossack et al., 2023).

The invasion history of bullfrogs and their potential impacts on native ecosystems have been studied for over a decade, but little is known about the introduction of bullfrogs in the southwestern state of New Mexico. Bullfrogs are considered invasive across most of the state, though there has been considerable debate about how far west the native range spans. In New Mexico, bullfrogs are not native to the Rio Grande and the southwestern area of the state, but the northeastern edge of the state may be considered part of the native range (Degenhardt, Painter & Price, 1996). The first potential introduction of bullfrogs in New Mexico was described in a 1921 Game and Fish report (NMDGF, 1921), in which the State Game and Fish Commission instructed the Game Warden to secure bullfrogs from the Kansas State Fish Hatchery for stocking of the Rio Grande. Therefore, bullfrogs may have been introduced from the native range intentionally by government agencies for recreational hunting at least once. Records from annual reports between the years 1963–1968 and 1970–1971 (NMDGF, 1963; NMDGF, 1964; NMDGF, 1965; NMDGF, 1966; NMDGF, 1967; NMDGF, 1968; NMDGF, 1970; NMDGF, 1971) state that bullfrogs were acquired, kept on hand, and disposed of, which could indicate that more bullfrogs were introduced to New Mexico, but the release locations and source populations are unknown. Based on museum records, bullfrog populations far exceed the initial Rio Grande valley in New Mexico and have spread to disconnected aquatic habitats likely due to human-mediated dispersal, as in Kamath, Sepulveda & Layhee (2016) in Montana. To further investigate these unknowns, population genetics can be used to give first descriptions of the invasion history of bullfrogs into New Mexico and other introduced regions.

Population genetic analyses have successfully furthered our understanding of invasion biology across numerous taxa (e.g., Baker & Stebbins, 1965; Lu et al., 2022; Thompson et al., 2024). For example, assessing haplotype frequency patterns across geographic space has proven helpful in identifying potential source populations and the number of introductions of zebra mussels into North America (May et al., 2006). Similar approaches have been applied to investigate global bullfrog introductions, demonstrating that source populations for bullfrogs in various European countries included midwestern and eastern native populations in the U.S. (Ficetola, Bonin & Miaud, 2008). Within the northwestern U.S., multiple introduction events from various native source regions were inferred and have likely contributed to the invasion success of bullfrogs in Montana (Kamath, Sepulveda & Layhee, 2016). Comparing the genetic diversity of invasive populations to their native counterparts provides insight into how genetic variation influences invasion success (Li et al., 2022). For example, genetic diversity can impact the success of invasive species by either limiting their ability to adapt to novel environmental stressors and putting them at a greater risk of extirpation when variation is low, or by enabling adaptive potential when variation is high (Kolbe et al., 2004).

Many invasive populations exhibit lower genetic variation compared to their native counterparts because of genetic bottlenecks associated with colonization (Estoup et al., 2016). Though this decreased variation should negatively impact their ability to colonize and adapt to new environments (Bouzat, 2010), many invasive populations can thrive. This scenario, known as the genetic paradox of invasions (Allendorf & Lundquist, 2003), has garnered sustained attention in recent years, with focused research into other mechanisms responsible for invasion success. For example, aspects such as phenotypic plasticity (Hagenblad et al., 2015) via transposable elements (Stapley, Santure & Dennis, 2015), host traits such as timing of reproduction (Urbina et al., 2020), or behavior related to recognizing predators (Garcia et al., 2012), have all been linked to invasion success despite low genetic diversity. One crucial factor in invasion success is the number of introductions, which can increase the genetic pool and promote successful establishment (Kolbe et al., 2004; Lavergne & Molofsky, 2007). Bullfrogs demonstrate low genetic diversity in invasive populations with a variable number of introductions and source populations across North America (Funk et al., 2011; Kamath, Sepulveda & Layhee, 2016) and in other global studies (e.g., China: Zhang et al., 2024), but they have been consistently successful in invading novel environments. These contrasting patterns highlight the complexity of identifying determinants of invasion success, emphasizing the need for further research to understand how genetic diversity, introduction history, and ecological context interact to shape establishment and spread.

The goal of this study was to use population genetics and phylogeographic approaches to uncover the invasion history of bullfrogs across their introduced range in the western U.S. with a focus on New Mexico. Specifically, we sequenced cytochrome b (cytb) in the mitochondrial genome of 77 individuals, sampled from multiple waterways across New Mexico. We then compiled and analyzed data from prior studies in other invasive and native regions and combined these with our new data to (1) assess the level of genetic diversity and identify haplotypes within bullfrog populations in New Mexico, (2) compare the genetic diversity of native and invasive bullfrog populations, (3) infer the number of introductions into New Mexico, and (4) identify potential native source populations. While we recognize that the inferential power of a single mitochondrial marker is limited, the use of this widely sequenced locus facilitated the synthesis of range-wide data and serves as an important first step for future invasive population studies of bullfrogs. By integrating genetic data from several previous studies into these analyses, this study provides insight into the mechanisms shaping bullfrog invasion history in the western U.S. and contributes to a broader understanding of how invasive species establish, spread, and persist in novel environments.

Materials & Methods

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.

Table 1: Population genetic diversity metrics of introduced bullfrogs in New Mexico.Diversity metrics reported as haplotype diversity (Hd) and nucleotide diversity (π). Population sample size (n) and number of haplotypes per population (nHd).

Haplotype and haplogroup distribution of bullfrogs sampled across the United States.(A) Cytochrome b haplotype proportions by sampling sites across New Mexico. Different haplotypes are represented by colors. (B) Haplotype network for cytochrome b haplotypes across the native and invasive range. Color indicates region where haplotype is found. Size of pie charts is proportional to haplotype freqency. Dashes on lines between haplotypes represents number of mutations. (C) Distributions of estimated bullfrog haplogroups according to phylogenetic analyses and haplotype network across invasive and native ranges. Silhouette of bullfrog was obtained from phylopic.org. U.S. County-level ranges were obtained from The International Union for Conservation of Nature’s Red List (IUCN, 2022).

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 (Hd) 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 Hd 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 (Hd 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 Hd 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.

Results

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).

Distribution of cytochrome b haplotypes found in New Mexico bullfrogs across native and invasive ranges.Colors represent identical haplotypes to invasive populations found in New Mexico. Other haplotypes not identical to those in New Mexico are shown in grey. Both invasive and native ranges are shown on the map as shaded regions and were originally obtained from The International Union for Conservation of Nature’s Red List (IUCN, 2022) and adjusted according to new and existing museum records. The native range was divided into four meaningful biological areas, visualized by the dashed lines, based on the nested clade analysis published in Austin, Lougheed & Boag (2004). Additional sequences from Funk et al. (2010), Kamath, Sepulveda & Layhee, 2016, Lafond et al. (2022), MacGuigan et al. (2022), Rodgers et al. (2023) were assigned to these regions according to geography (Northeast, East, Overlap, Midwest). The invasive range was split into two regions (Northwest and Southwest) based on geographical boundaries.

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: Hd = 0.830, π = 0.009; CF: Hd = 0.833, π = 0.0034; Table 1; Fig. 1A). We observed the lowest genetic diversity in the Upper Gila River population (UGR: Hd = 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.

Distribution of genetic diversity metrics for native and invasive populations.Distribution of haplotype diversity (A) and nucleotide diversity (B) of native and invasive populations. Mann–Whitney U test statistics (W) and significance values (P) are provided for each comparison.

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.

Table 2: Pairwise ΦST and corresponding P-values between six introduced New Mexico populations.Pairwise ΦST values are found above the diagonal line and corresponding P-values are found below the diagonal line. Bolded text indicates significance after sequential Bonferroni correction for 15 total tests with a new significance level 0.0033.

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).

Table 3: Summary of analysis of molecular variance (AMOVA).Results show the percentage of among-group variation between the introduced New Mexico population and the four native (i.e., East, Northeast, Overlap, and Midwest) and two introduced (i.e., Northwest and Southwest) regions. Lower percentages indicate higher genetic similarity, while higher percentages indicate lower genetic similarity among groups. Bolded text indicates significance (P < 0.05).

Discussion

We investigated genetic differences of bullfrogs in New Mexico compared to other introduced and native populations to elucidate aspects of their invasion history. We found patterns consistent with a single introduction and recent expansion throughout the state of New Mexico, with little cytb differentiation among localities and most haplotypes matching those found in other native and invasive regions. Across New Mexico, genetic diversity was highest in localities with more human activity. Similar to previous studies of bullfrog invasion, we found low genetic diversity across introduced regions when compared to native populations, likely due to genetic bottlenecks associated with founder effects. While we were unable to identify the source population of New Mexico bullfrogs, in part because of limitations of the cytb marker, we found no significant among-group differences between New Mexico and other introduced western populations, which indicates a possible common origin or recent connectivity. Though records of bullfrog introduction into the western U.S. are disparate, our study provides new data and comparisons that suggest recent spread throughout the western U.S. and highlights the need for further research with higher-resolution genomic data.

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.

Conclusions

The invasion history of bullfrogs across the globe has been widely studied, but little work has focused on the Southwestern U.S. Here, we provided the first genetic description of the bullfrog invasion history in this region and placed it in context with other invasive populations. Our results suggest successful bullfrog establishment may not be constrained by low genetic diversity. Specifically, we found that invasive populations consistently had lower genetic diversity than native populations. We also found limited genetic differentiation among sites in New Mexico, which is consistent with a single introduction event. This suggests that a single invasive founding population could have originated from multiple sources across the United States which contributed sufficient genetic diversity for widespread establishment. Given the widespread establishment of bullfrogs across New Mexico, future work should investigate the role of secondary spread. Additional genetic markers are needed to clarify a potential source population and help predict movement across the Southwest. Overall, our results provide insights into the connectivity of invasive bullfrog populations in New Mexico and highlight the importance of understanding the history of invasive species to managing existing populations.

Supplemental Information

10.7717/peerj.20491/supp-1Supplemental Information 1Cytb Phylogenetic treeMaximum likelihood phylogeny estimated with cytochrome b sequence datasets shown as cladograms. Bootstrap values >70 are shown. All sequences are named after Genbank accession numbers and corresponding haplotype. The tree is rooted using Rana heckscheri and Rana clamitans as outgroups.

10.7717/peerj.20491/supp-2Supplemental Information 2Metadata of all New Mexico specimensThis file contains Metadata of all New Mexico specimens used in this study from the Museum of Southwestern Biology. The file also contains links to the Arctos museum database where all information is readily available. The first sheet within the file contains important data for each specimen. Definitions of column headers are included on the second sheet.

10.7717/peerj.20491/supp-3Supplemental Information 3Cytb sequence dataThe data for sequences per site used in the study including diversity indices, sample sizes and geographical data. This also exhibits the number of populations compared, the AMOVA region, and which study each population was from.