Evaluating the effects of red imported fire ants (Solenopsis invicta) on juvenile Houston Toads (Bufo [=Anaxyrus] houstonensis) in Colorado County, TX

Study site

We conducted this study on the 4,265 ha APCNWR (Colorado County, Texas) which is about 120 km inland from the Gulf of Mexico (Morrow et al., 1996). The refuge is predominantly prairie grasslands with loamy prairie, sandy prairie, and coarse sand range sites (Morrow et al., 1996). This site is a unit of the USFWS National Wildlife Refuge System and represents one of the largest remaining patches of coastal prairie habitat. It is managed for the federally endangered Attwater’s prairie-chicken (Tympanuchus cupido attwateri) through prescribed fires, cattle grazing, native seeding, and RIFA control. The climate in these coastal prairies is subtropical with average annual precipitation ranging from approximately 750 mm to 1,100 mm (Kwiatkowski, Saenz & Hibbitts, 2017). The San Bernard River borders APCNWR to the east. For our field experiment, we selected four prairies (Fig. 1) within APCNWR, each of which were coarse sand range sites to enable burrowing by Houston Toads (Caldwell, 2015). These open grasslands were dominated by mid-tall native grass species (Morrow et al., 2015).

Experimental design

We conducted this study over a total area of 739.5 ha. Of that total, 261.5 ha were left untreated and 478 ha were treated for RIFA using 1.7 kg/ha of Extinguish^® Plus (0.25% s-methoprene, 0.36% hydramethylnon; Central Life Sciences, Schaumburg, IL) on 22 October 2014 (Fig. 1). An insecticide and insect growth regulator in this bait-driven suppressant is reported to disproportionately affect RIFA in comparison to native ant species (Calixto et al., 2007). A helicopter was used to place aerial applications of Extinguish^® Plus on two (i.e., treated zones) of the four prairies selected for our study (Fig. 1). We adhered to manufacturer’s protocol with regard to the application rate of the insecticidal bait (Wellmark^®/Central Life Sciences of Central Garden & Pet). RIFA from small colonies have been observed to forage no further than 3 m for food while those from large colonies showed foraging distances of up to 17 m (Stringer et al., 2011). To account for the upper-bound of foraging distances and ensure no overlap between treatments, we created a 50-m buffer from the edge of each treatment. Further, we also restricted the random assignment of experimental unit locations within 250 × 500 m plots that were approximately within the center of each of our four prairies. This sought to control for any overspray of insecticidal bait from treated into untreated prairies, or any adjacency effects of the treatments themselves.

Assessing RIFA and native invertebrate density

We assessed RIFA densities among treated and untreated prairies starting on 8 April 2015, every two weeks through July 2015. We evaluated RIFA activity using 6 bait traps per prairie. Each bait trap consisted of pureed, canned sausages in a petri-dish and we placed these bait traps at randomly assigned locations between access roads and the edge of our treatment sites. Bait traps in each prairie were placed at least ∼10 m apart from each other. Given that these were baits of meat and likely to be monopolized by a dominant, predacious species (Greenslade & Greenslade, 1971), we restricted trap locations, as above, so as to not directly attract RIFA to our outdoor experimental units (i.e., toad exclosures). Bait traps were set for a period of ∼3 h (∼10:00AM–1:00PM). We concurrently sampled for surface-active native invertebrates using wet pitfall traps. These were High Density Polyethylene cups of ≈120 ml volume that contained propylene glycol (commercial non-toxic antifreeze). Pitfall traps were placed 15 m north of toad exclosures in each prairie, to prevent attracting ants to those exclosures. Pitfall traps were open for a duration of 24 h, though not on the same day as bait traps. We expected pitfall traps were less likely to be dominated by predacious species because this is a more passive sampling technique (Greenslade & Greenslade, 1971). We identified native ants to the genus level and classified non-ant native invertebrates as small (1 to 12 mm body length) or large (>12 mm body length). We counted the total number of RIFA, native ants, and non-ant native invertebrates in each pitfall trap. These pitfall traps were primarily used to quantify the number of available prey items (i.e., ground-dwelling/surface active native invertebrates) for Houston Toads with respect to RIFA treatment. Locations for bait and pitfall traps were marked by surveyor flags and recorded in a handheld GPS unit to enable repeated visits.

We pooled observations for RIFA, native ants, and non-ant native invertebrates by treated and untreated sites, from bait and pitfall traps. The number of observations from pitfall and bait traps (n = 194) were an order of magnitude lower than repeated measures of individual toads (n = 1,339). This difference was because our sampling intensity for toads was higher than invertebrates and did not allow us to use counts of RIFA directly as a predictor of growth and survival in juvenile toads. Such a direct comparison would lose the variation inherent in our toad dataset by pooling repeat measurements or erroneously pseudoreplicate our RIFA counts. For this reason, we employed a syllogistic approach of first testing for differences in counts of RIFA and native ants and non-ant invertebrates among prairies that were treated and untreated for RIFA and then testing for a treatment effect in toad growth and survival.

All analyses were performed using Program R (R Core Team, 2016). Where relevant, we inferred significance at α = 0.05. Mean values are reported with standard error. We fitted separate GLMMs, using Poisson and Negative Binomial response distributions with and without zero-inflation, to our RIFA, native ant, and non-ant native arthropod counts in the R package glmmADMB (Fournier et al., 2012). We attempted to account for repeated measures of our prairies by including individual prairies as a random factor. Candidate models included random variation in intercepts and random variation in slopes and intercepts. We evaluated model fit using AIC_c values in package MuMIn (Barton, 2017). Models with a ΔAIC_c ≤ 2 were considered similar to models that best supported the data. Once we had chosen optimal random structure, we assessed error distributions and evaluated our fixed factors of treatment and time (i.e., week). A subset of candidate models (Tables S1–S3) specified an interaction among time (i.e., week) and treatment to account for continued ant suppression over time and the expected benefits of such RIFA decreases on native invertebrate assemblages (Calixto et al., 2007). The fixed effect structure that best supported the data was chosen using AIC_c values as above.

Juvenile Houston toad growth and survivorship

We treated 478 ha of our study area for RIFA using 1.7 kg/ha of Extinguish^® Plus (0.25% s-methoprene, 0.36% hydramethylnon; Central Life Sciences, Schaumburg, IL) on 22 October 2014. We began our study in March 2015, which was within the estimated 12-month duration of RIFA control by insecticidal bait (Wellmark^®/Central Life Sciences of Central Garden & Pet). We used outdoor terrestrial experimental units (Hopkins, Mendonça & Congdon, 1997) to track growth and survival in juvenile Houston Toads at our treatment sites. We assigned locations for ten experimental units in a 250 × 500 m area of each prairie via random generation of spatial coordinates, for a total of 40 experimental units (Fig. 1).

Each experimental unit, hereafter exclosure, was a 150 × 150 cm area of the ground enclosed by hardware cloth (3.2 mm mesh size, 0.36 mm wire diameter; Fig. 2). We refer to our experimental units as exclosures for their ability to exclude potential vertebrate predators (e.g., snakes, birds). Exclosure walls extended 20 cm below ground and 50 cm above ground. We folded a 10 cm lip of hardware cloth along the length of the top and bottom of each exclosure to prevent toads from escaping by tunneling under or scaling exclosure walls to climb out (Jones, 2015). Additionally, we placed deer netting on top of each exclosure to prevent entry of predators. However, the hardware cloth we used had a mesh size (3.2 mm) that was large enough to enable RIFA to move in and out of exclosures. We mitigated chances of desiccation of individuals by placing two small plastic water tubs full of organic compost and sphagnum moss within the artificial containment. We placed each tub under or in the shelter of the native bunch grasses within each of the exclosures. We recorded exclosure locations in handheld GPS units, which enabled repeated visits to the exclosures.

We randomly assigned and released 200 full sibling, juvenile toads in exclosures on 25 March 2015. These juveniles were sourced from the captive breeding program at Houston Zoo. We aimed to minimize heterogeneity by using full sibling toads (i.e., from the same clutch). Houston Toads were used immediately after metamorphosis was complete (defined as total absorption of tail) and did not appear to show considerable difference in size among exclosures in treated (Body Mass (W in g) = 0.79 ±  0.03 g; Snout Urostyle Length (SUL in mm) = 17.89 ± 0.26 mm; Mean ± SE) and untreated sites (W = 0.75 ± 0.04 g; SUL = 18.2 ± 0.35 mm; Mean ± SE). Juvenile release densities in each exclosure were based off a previous study conducted in canopied habitat in Bastrop County, Texas (Jones, 2015). Jones (2015) observed that survivorship in Houston Toads was best explained by release density, with a release density of five (ϕ = 0.93) toads per exclosure having the highest survivorship followed closely by exclosures that had six (ϕ = 0.92) and four (ϕ = 0.91) toads. Recapture rates in exclosures that had four (p = 0.91) and six (p = 0.92) toads were higher than those exclosures which had five (p = 0.87) toads (Jones, 2015). We released four juveniles per exclosure in 5 exclosures per prairie (n = 20 exclosures) and six juveniles per exclosure in the remaining 5 exclosures per prairie (n = 20 exclosures). We varied juvenile densities in exclosures as above because these densities potentially ensured study longevity via a trade-off between survival and recapture rates (Jones, 2015) and the number of available captive toads. Secondarily, variation in release density also allowed us to extend a test of these previously determined optimal release densities (Jones, 2015) in coastal prairie habitat. We marked each individual by clipping toes using sterile dissection scissors in correspondence to a specific numbering system (i.e., number 0, 1, 2, 3 for 4 juveniles per exclosure and 0, 1, 2, 3, 4, 5 for 6 juveniles per exclosure; (Ferner, 1979). We weighed and measured juveniles prior to placing them in exclosures. Marked individuals allowed for tracking of individual growth and survival. Our study was authorized by US Fish and Wildlife Service (Permit No. TE039544-0) and Texas Parks and Wildlife Department (Permit No. SPR-0102-191). Our study protocol was approved by the Texas State University Institutional Animal Care and Use Committee (Protocol No. IACUC201648186).

Previous studies report reduced survival in translocated amphibians immediately following release (Cooke & Oldham, 1995; Tocher & Pledger, 2005). For this reason, we conducted daily searches of exclosures for toads over five days immediately following release and sprayed water tubs to keep them moist. Following the initial daily searches, data were collected from the individuals for approximately six months; we recorded W in g, SUL in mm, and head width (HW in mm) for each individual recaptured after searching the exclosure for five minutes. We searched exclosures and measured toads weekly from Week 2 to 6. Subsequently we searched exclosures every two weeks from Week 8 to 16. We conducted the last two exclosure searches on Week 20 and 24. Toad desiccation would likely have been exacerbated by artificial containment and so we mitigated this by moistening water tubs in the exclosures after five days of no precipitation. We sampled more frequently at the start of the study in order to monitor for low survival following release and to account for an inherently low survival rate in Houston Toads of a smaller size class (Cooke & Oldham, 1995; Tocher & Pledger, 2005; Greuter, 2004; Swannack, Grant & Forstner, 2009). As Houston Toads increased in size and thereby survivorship, we decreased sampling frequency to reduce the level of disturbance (Greuter, 2004; Swannack, Grant & Forstner, 2009). Exclosures were searched by multiple observers, including Texas State University students, Student Conservation Association (SCA) interns, and APCNWR staff members. In order to specifically control for observer bias in measuring toads, all observers were trained prior to their participation in the study. In addition, the role of taking measurements was restricted to the same subset of observers over repeat occasions.

Juvenile body mass is likely to vary as a function of both growth and hydration status (Todd & Rothermel, 2006). We therefore used SUL in mm as a more stable proxy for growth. We used linear mixed effects models in package lme4 (Bates et al., 2015) to address the primary question on whether juvenile growth varied among treated and untreated sites, while attempting to account for other sources of variation. We used the beyond optimal approach as described by Zuur et al. (2009), i.e., we included as many explanatory variables in the fixed component of the model- so that the random component did not have any information that we would instead need in the fixed component of the model. We fitted models with all explanatory variables and as many interactions as possible, i.e., a three-way interaction of fixed factors of treatment, time (i.e., week), and release density while varying random factors to determine optimal random effect structure. Models were fit using Restricted Maximum Likelihood (REML) estimation and optimal structure of random factors was assessed using Akaike Information Criterion scores corrected for a small sample size (AIC_c) in package MuMIn (Barton, 2017). Models with a ΔAIC_c ≤ 2 were considered similar to models that best supported the data. Once the structure of our random factor was chosen, we fitted models with varying fixed factor structure using Maximum Likelihood (ML) estimation and assessed fixed factor structure using AIC_c values as above. We then presented our best-fit model using REML estimation, where t tests on such estimates used Satterthwaite approximations to degrees of freedom (Zuur et al., 2009).

We accounted for non-independence of repeated measures of toads in exclosures by testing each individual toad, each exclosure, and each individual toad within an exclosure as a random factor. Candidate models included random variation in intercepts and random variation in intercepts and slopes (Table S4). We tested for potential differences in SUL among fixed factors of treatment (i.e., treated or untreated), initial/release density (i.e., four or six juveniles), and time (i.e., week), as well as higher order interactions of these fixed factors (Table S4).

We estimated juvenile survivorship by applying a Cormack–Jolly–Seber model to live recapture data in Program MARK (Cormack, 1964; Jolly, 1965; Seber, 1965; White & Burnham, 1999). We used the RMark package to interface with Program MARK from the R environment (Laake, 2013; Mazerolle, 2015). Program MARK employs information-theoretic methods in selecting the best supported model estimates of apparent survival (ϕ) and recapture probability (p). We explored variation in ϕ and p using time and treatment (i.e., treated and untreated sites) as well as higher-order interactions of time and treatment as explanatory variables. Additionally, we also included models that considered either or both parameters as constant among sampling periods. We developed 14 candidate models (Table S5), with each model representing a distinct biologically-relevant hypothesis. We estimated the overdispersion parameter ($\hat{c}$) by dividing the global χ² by its corresponding degrees of freedom (Mazerolle, 2015). The overdispersion parameter was incorporated to provide Quasi-likelihood Akaike Information Criterion scores corrected for a small sample size (QAIC_c), based off which model fit was evaluated. Models with a ΔQAIC_c ≤ 2 were considered similar to models that best supported the data; however, we report models over a cumulative weight of ≈0.95.

Assessing RIFA and native invertebrate density

The mean number of RIFA per sample in untreated sites (5.89 ± 1.2) were nearly 9 times greater than that in treated sites (0.695 ± 0.21 ants; Fig. 3). Our best-fit model (Table S1) was one that fit RIFA counts to a type-two negative binomial distribution (i.e., quadratic mean–variance relationship) with intercepts randomly varying among prairies. This model showed a greater rate of decrease in the number of RIFA in treated areas relative to untreated areas (β = −0.257, SE = 0.09, Z =  − 2.93, P = 0.003; Table 1). Standard deviation among prairies was 0.0003.

Counts of native ants across 11 genera as well as counts of native ants whose taxonomic identities could not be ascertained were used. Mean number of native ants per sample were nearly 5 times greater in untreated sites (30.80 ± 6.10) compared to treated areas (6.28 ± 0.96; Fig. 3). A majority of observed native ants were cone ants (Dorymyrmex sp.), crazy ants (Nylanderia sp.), and Forelius sp. ants at treated (93.8%) and untreated (92.2%) sites. We failed to achieve convergence with the model that specified random variation in intercepts and slopes among prairies. The top-ranked model fitted native ant counts to a type-two negative binomial distribution, with intercepts varying randomly among prairies (Table S2). The top-ranked model solely specified a treatment effect with counts of native ants significantly lower in treated sites (β = −1.552, SE = 0.25, Z =  − 6.11, P <0.001; Table 2). Standard deviation among prairies was 0.158.

We observed nearly equal numbers of native non-ant invertebrates in treated (11.51 ± 1.46) and untreated sites (11.52 ± 1.23; Fig. 3). We observed beetles, crickets, juvenile grasshoppers, snails, spiders, and weevils among other native non-ant invertebrates. A majority of captures at both treated (97%) and untreated (95%) sites were small non-ant native invertebrates. Our top-ranked model fit these counts to a type two negative binomial distribution with intercepts randomly varying among prairies (Table S3). The best-fit model did not specify a treatment effect, instead specifying that counts of non-ant invertebrates decreased over time (β = −0.034, SE = 0.013, Z =  − 2.57, P = 0.01; Table 3).

Juvenile Houston toad growth and survivorship

We released a total of 200 juvenile Houston Toads, with 100 juveniles each released in exclosures in sites treated and untreated for RIFA. The optimal random structure for our models allowed random variation in intercepts and slopes of exclosures and individual toads within exclosures, proving a better fit (ΔAIC_c = 68.36; Table S4) than the next best model which allowed random variation in slopes and intercepts of individual toads alone. In terms of fixed effects structure, we observed three competing models (ΔAIC_c ≤ 2; Table 4). The model with the lowest AIC_c score included a higher order interaction of release density with time, while the next model solely included time as a fixed factor (Table 4). We report from the top-ranked model that juvenile toads in exclosures with a lower release density showed a significantly greater rate of increase in SUL (β = 0.245, SE = 0.110, t_34.5 = 2.226, P = 0.033; Table 5; Fig. 4). We saw more variation among individual toads within exclosures (SD = 2.09) and among exclosures (SD = 2.15) relative to the residual standard deviation (SD = 1.6).

The model that best fit our live-recapture data and received majority of support stipulated constant apparent survival and recapture probability that varied by time (ϕ_(.) p_(t); Table 6). The next best and competing model specified time-varying apparent survival and treatment-varying recapture probability (ϕ_(t) p_(g);ΔQAIC _c = 0.33; Table 6). We report estimates from our top-ranked model. Apparent survival was estimated at 0.916 ± 0.008. We approximate overall apparent survival over the 24-week period to 0.12 (i.e., 0.916²⁴). Recapture probabilities ranged from a minimum of 0.164 ± 0.084 to 0.922 ± 0.026 and were lowest over the last two sampling occasions (i.e., Week 20 and Week 24; Fig. 5).