Genetic estimates of contemporary effective population size in an endangered butterfly indicate a possible role for genetic compensation.

The effective population size (N e ) is a critical evolutionary and conservation parameter that can indicate the adaptive potential of populations. Robust estimates of N e of endangered taxa have been previously hampered by estimators that are sensitive to sample size. We estimated N e on two remaining populations of the endangered Miami blue butterfly, a formerly widespread taxon in Florida. Our goal was to determine the consistency of various temporal and point estimators on inferring N e and to determine the utility of this information for understanding the role of genetic stochasticity. We found that recently developed 'unbiased estimators' generally performed better than some older methods in that the former had more realistic N e estimates and were more consistent with what is known about adult population size. Overall, N e /N ratios based on census point counts were high. We suggest that this pattern may reflect genetic compensation caused by reduced reproductive variance due to breeding population size not being limited by resources. Assuming N e and N are not heavily biased, it appears that the lack of gene flow between distant populations may be a greater genetic threat in the short term than the loss of heterozygosity due to inbreeding.

Introduction

The concept of effective population size (N) is a critical parameter in evolution and conservation because it, not census size, indicates the rate of loss of heterozygosity in finite populations (Wright 1931; Caballero 1994; Allendorf and Luikart 2007). The N is the size of an idealized population that has the same rate of change in heterozygosity (or inbreeding) of the real population (Crow and Kimura 1970). The idealized population is typically much smaller than the real population because the latter rarely behaves in an ideal fashion (i.e., having equal sex ratios, constant population size, discrete generations, and an equal contribution of individuals to the next generation) (Fisher 1930; Wright 1931). Ratios of N/N are typically expected to be low in wild populations. Meta-analyses have found N/N ratios ranging from ∼0.11 to 0.50 depending on the life history of the organism and method used to generate N (Nunney and Elam 1994; Frankham 1995; Nunney 1996). Low N/N ratios may indicate that population health is at risk of demographic contraction due to lack of genetic variation even if the number of individuals in the population is still large.

One method for calculating N has been through calculation of demographic parameters that requires knowledge of sex ratios, family size variance, and population size fluctuations (Frankham et al. 2003; Wang and Whitlock 2003) and, as such, collecting sufficient data for threatened taxa has proven especially difficult (Harris and Allendorf 1989; Frankham 1996). In contrast, samples for molecular estimates can be relatively easy to collect. However, molecular estimators have traditionally been characterized as having low precision due to the stochastic nature of neutral markers in small populations (Wang 2005) and the difficulty in sampling temporally in many taxa. The latter point is related to the two primary approaches for calculating N values from genetic data. Temporal estimators use allele frequency data from multiple samples, across generations, to calculate N (Nei and Tajima 1981; Pollack 1983; Waples 1989; Wang 2001). In contrast, single point estimators use linkage disequilibrium or heterozygote excess to estimate N but have not been used as much as temporal methods due to relatively imprecise and biased estimates, especially when sample sizes were small (Waples 1991; England et al. 2006; Tallmon et al. 2008). With respect to the latter, one can estimate using a linkage disequilibrium method that measures departures from expected proportions (Weir and Hill 1980; Waples 1991; Bartely et al. 1992). The heterozygote excess method compares expected Hardy-Weinberg values to the increases in the observed number of heterozygotes (Pudovkin et al. 1996). The main advantage of single-point estimators is that one is able to calculate N from a single generation of data, a strategy that may be especially useful for longer-lived organisms for which multi-generational sampling would be prohibitively difficult.

With respect to the utility of estimating N from molecular data for species of conservation concern, numerous, highly polymorphic markers are also becoming more prevalent and new programs are able to capitalize on these advancements and make single point estimators more useful. However the irony remains that understanding genetic stochasticity in small populations is hampered by the ability of genetic estimators to accurately estimate N with small sample sizes. Hence, estimating N has been particularly problematic for most conservation applications due to limited sampling and the vagarities of calculating N from molecular data. Recent computational advancements may permit more reliable and biologically meaningful estimates of N (Peel et al. 2004; Palstra and Ruzzante 2008; Tallmon et al. 2008; Waples and Do 2008) that in turn will provide important information for wild and ex situ conservation.

Understanding the time period to which an estimate of N applies is critically important but is often erroneously interpreted (Waples 2005). Proper interpretation of results of molecular estimates of N depends on the method used to sample populations (e.g. type I or II, see Table 1) and calculate N, as well as the life history of the organism. This last factor is of importance as reproductive variance, sex ratios, and fluctuations in population size all affect N and, as such, must be considered in the interpretation of results and comparisons across taxa. Further, changes in biotic interactions at lower population sizes (genetic compensation) and environmental stochasticity can inflate N/N (Palstra and Ruzzante 2008). Whether we can with any confidence infer patterns of N for highly threatened species remains an important area of empirical research. Our approach here is to explore the utility of N estimators on a highly endangered butterfly that fits a number of useful criteria.

Table 1

Summary of N estimators

Program	Description	Comments/limitations	Methodological reference
NeEstimator	2-sample, moment-based	Variance effective size estimator based on change in F* over generational samples.	Peel et al. (2004); Waples (1989)
		Accommodates for type I and II sample schemes†; estimates the harmonic mean if Ne is not constant.
MNe	2-sample, pseudo-likelihood	May be sensitive to skewed allele frequencies, overestimating Ne.	Wang (2001); Wang and Whitlock (2003)
		Also allows for joint estimate of Ne and m. User defined maximum Ne as input.
		Accommodates type II sample schemes.
TempoFs	2-sample, moment-based	Alleles are weighted to reduce bias in Ne associated with high polymorphism.	Jorde and Ryman (2007)
		Large standard deviation of F.
		User defined estimate of census size.
		Accommodates for type I and II sample schemes.
ONeSAMP	1-sample; uses approximate Bayesian computation; Web-based	User defined priors of Ne.	Tallmon et al. (2008)
		Generates 50,000 simulated populations based on user data, summary statistics close to observed data delineates accepted range of Ne.
LDNE	1-sample; Ne estimated based on linkage disequilibrium (Burrow's Δ)	Random mating and monogamous systems.	Waples and Do (2008)
		Separate estimates accommodating rare alleles.
		Confidence assessed via jackknife and parametric CIs.
		Corrects for bias associated with small samples sizes.

F is the standardized variance in allele frequency change (Waples 1989).

Sample scheme I samples from adults with replacement (or following reproduction); scheme II samples before reproduction without replacement (Jorde and Ryman 2007; Nei and Tajima 1981; Waples 1989).

The Miami blue butterfly (Cyclargus thomasi bethunebakeri) was once a common endemic subspecies distributed across southern Florida (Minno and Emmel 1993; Calhoun et al. 2002; Daniels et al. 2008). Sightings decreased during the 1980s, until it was feared to have been extirpated after decades of habitat loss and fragmentation from both natural and anthropogenic forces (FWC 2003). Two isolated populations were recently discovered in the southernmost part of the state; one, a small population in Bahia Honda State Park (BHSP) in the Lower Florida Keys, the other on the uninhabited Marquesas Islands and Boca Grande Island in Key West National Wildlife Refuge (KWNWR), 50 km west of Key West (Ruffin and Glassberg 2000; Cannon 2007) (Fig. 1). Management concerns for the two remaining populations include habitat alteration from hurricanes, impact of drought on host plants, predation, illegal collecting, and the loss of genetic diversity (FWC 2003).

Figure 1

Historical range (inset) and current locations of Miami blue butterfly populations; BHSP, Bahia Honda State Park; KWNWR, Key West National Wildlife Refuge.

Given the precarious nature of the remaining butterfly populations, proper management should incorporate both demographic and genetic aspects because the loss of genetic variation due to inbreeding can adversely affect population persistence (Saccheri et al. 1998; Spielman et al. 2004). Here we present results from a comparative examination of multiple contemporary N estimates on the two remaining Miami blue populations. Our objectives are to calculate contemporary N from both temporal (eight generations) and point samples and to compare and contrast values (as appropriate) from these different methods. Where appropriate, we then compare these values to census counts (N) from both populations to determine how much genetic variation has been retained over the recent range contraction. We also estimate the N from samples taken from the captive colony of Miami blue butterflies and compare these with N estimates from existing natural populations.

Methods

Study populations and natural history

The habitat of Miami blue butterflies on BHSP and KWNWR (Fig. 1) is defined by ephemeral patches of host plants Caesalpinia bonduc and Pithecellobium keyense (Fabaceae). Because of the 70 km separating the two populations and the limited dispersal ability and longevity of adults (below), gene flow between the two populations is unlikely, particularly over a time scale relevant to our data, and these remnant populations are considered closed to immigration. Since 2002 over 30 000 individuals have been produced in captivity (McGuire Center for Lepidoptera and Biodiversity, University of Florida) and over 7000 have been released into suitable habitat within state and national parks in the South Florida mainland and northern Florida Keys (Daniels et al. 2008). However, these release sites are sufficiently far from the two remaining populations to preclude the possibility of gene flow (Fig. 1). Extensive studies of butterflies in BHSP and in captivity have revealed a generation time of 4–6 weeks. Females lay over 100 eggs in their lifetimes although variance is high in captivity (120.46 ± 85.62 eggs, mean ± SD, N = 24) (E.V. Saarinen, unpublished data). Adult butterflies survive less than a week in the wild, but generations are not completely discrete owing to the variance in development time of all life stages. There are approximately 8–10 generations per year in the wild. Adults mate soon after eclosion and females lay eggs on the developing leaves, shoots, and flower buds of C. bonduc or P. keyense. Larvae from both populations develop over the course of 3–4 weeks and pupate for 7–12 days on their host plants (J.C. Daniels, unpublished manuscript).

Genetic sampling and microsatellite fragment generation

Adult butterflies were caught with hand-held aerial insect nets and were nonlethally sampled by removing an approximate 2 mm2 portion of hind wing with forceps (type I sampling, Nei and Tajima 1981). Wing tissues were placed in 90% ethanol and stored at −80°C. We sampled BHSP at two time periods (September 2005, N = 24 and June 2006, N = 39) and at KWNWR during one period (February 2008, N = 27). We extracted DNA from wing fragments using the Qiagen DNeasy Kit (Qiagen©, Valencia, CA, USA) following the manufacturer's protocols but with a final, one time elution of DNA in 50 μL of 55°C 10 mM TRIS, and stored DNA extractions at −20°C. We genotyped butterflies at 11 trinucleotide and one tetranucleotide polymorphic microsatellite markers. Polymerase chain reaction (PCR) conditions and genotype analyses are described in Saarinen et al. (2009).

We detected multiple private alleles between sampling periods and geographic locations. Approximately 10% of all BHSP samples were re-genotyped at nine microsatellite loci (those containing private alleles) to confirm genotypes. All private alleles were individually assessed for validity that they were not the result of genotyping or scoring error. Micro-Checker version 2.2.3 (van Oosterhout et al. 2004) was used to test for the presence of null alleles as problems of allelic drop-out have been identified for PCR-amplified products in lepidopteran microsatellites (Zhang 2004) as well as for studies utilizing noninvasive genetic samples of low quantity (Taberlet et al. 1999). We used FSTAT version 2.9.3.2 (Goudet 1995) to test for linkage disequilibrium and examined departures from Hardy-Weinberg equilibrium using exact tests implemented in Genepop 34.0 (Rousset 2008).

We tested for the relative contribution of stepwise mutations (SMM) and infinite allele model (IAM) to global and pairwise population differentiation using the randomization test of Hardy et al. (2003). Briefly, this test examines the contribution of stepwise mutation, relative to drift, to population differentiation by randomizing allele sizes within a locus while maintaining genotypic states of individuals. If allele size shifts are resulting predominantly from stepwise mutations, then the observed estimates of R should be greater than those estimated from the permutated data set (pR). We conducted global population tests running 1000 permutations each, using SPAGeDi version 1.2 (Hardy and Vekemans 2002).

Data analysis

As estimates of N can be affected by recent population size changes, we first tested for evidence of recent bottleneck in each sample (temporal and geographic) using the program Bottleneck version 1.2 (Cornuet and Luikart 1996; Piry et al. 1999). Bottleneck tests for excess heterozygosity using several mutation models including the stepwise mutation model (SMM), infinite allele model (IAM), and a two phase model (TPM). For the TPM a proportion of mutations consist of more than one-step. We set the proportion of SMM at 70% for the TPM, and ran 1000 iterations for each of the three models. As a control, we also tested a known population bottleneck using data from the captive breeding colony of Miami blue butterflies housed at the University of Florida. This colony was founded in June 2006 from 80 eggs collected from eight discrete habitat areas on BHSP. Butterflies were allowed to breed as a closed population in captivity for four generations before being used to test for bottleneck. We genotyped 40 individuals from generation one, and 10 individuals (only 10 were sampled before adults were released) at generation four for each of the 12 microsatellite loci.

Temporal estimators

For the BHSP samples we applied three temporal estimators of the harmonic mean of N spanning September 2005 to the parental generation of June 2006 sample. At this temporal separation (∼8 generations) the bias due to none-discrete generations is expected to be insignificant (Waples and Yokota 2007). We used a moments-based temporal estimator (Waples 1989) implemented in NeEstimator version 1.3 (Peel et al. 2004). We also used a pseudo-likelihood method implemented in MNe version 2.0 (Wang 2001; Wang and Whitlock 2003). We used MNe with migration set at zero under which conditions MNe estimates the moment estimator of Nei and Tajima (1981) and the likelihood estimator of Wang (2001). For our MNe run we set the maximum N allowed at 1000. Finally we used a more recent estimator developed by Jorde and Ryman (2007) implemented in TempoFs that (unlike the previous two methods) reduces the bias associated with small samples sizes and skewed allele frequencies. For TempoFs, we estimated the harmonic mean census sizes (N = 21) observed from BHSP between August 2005 and May 2006 as these months reflect the time period to which the estimate of N represents (Waples 2005).

Point estimators

We also used two different single-sample (point) estimation methods to calculate N at T-1 for each temporal BHSP sample and the KWNWR sample. OneSamp 1.1 (Tallmon et al. 2008) uses approximate Bayesian computation to estimate variance N from summary statistics that are related to N. Using user-defined estimates of N priors, OneSamp generates 50 000 simulated populations drawn randomly from the distribution of N priors. Samples are drawn from each simulated population with the identical numbers of individuals and loci as the actual data set. Summary statistics are calculated and compared to the actual data and similar summary statistics are retained for use in generating an estimate of N using weighted local regression (Tallmon et al. 2008). We used upper and lower bounds on the prior for N of 2–100, 4–200, and 6–500 to explore the impact of widening priors on N estimation. Each set of priors was estimated across three replicates. OneSamp requires loci that are variable and does not support the use of loci with large numbers of missing alleles. As a result, we removed locus CthC124 from both BHSP samples due to lack of variation in 2005, and locus CthB117 and CthC116 from KWNWR due to too much missing data. Finally, LDNe 1.31 (Waples and Do 2008) uses linkage disequilibrium (LD) information among alleles at different loci caused by genetic drift in finite populations. The method does not assume random mating and corrects for biases associated with small sample sizes (England et al. 2006; Waples and Do 2008). We estimated N for varying levels of inclusion of rare alleles in order to examine potential bias contributed by low frequency microsatellite alleles. We include a summary of software programs to calculate N used in this study (Table 1).

Survey estimates of census size

Adult butterflies were surveyed on BHSP from October 2002 through August 2008 (Fig. 2). Population abundance estimates and the maximum number of individuals observed was determined from Pollard walk transects, a standard method to inventory butterflies in remnant habitats (Pollard 1977). A series of fixed-route transects were established along existing park roads and trails within the core of the metapopulation on the south end of BHSP to minimize habitat impact. Transect length varied with site area, ranging from 165 to 580 m. BHSP was typically sampled on a monthly basis throughout the duration of the study. On KWNWR, butterflies were surveyed in February 2008 by the checklist method. Checklist surveys are employed primarily to confirm the presence of individual species and sometimes the number of individuals at each survey site (Royer et al. 1998). Although procedurally similar to Pollard surveys where individuals walk designated transect routes, the checklist method enables the recorder to wander freely within the habitat. The increased flexibility provides for rapid survey response to local environmental conditions (e.g., changes in wind direction or speed, nectar availability, host quality, etc.) and enables the recorder to seek out more preferable sites for butterflies. Such survey plasticity was particularly critical due to dense, mangled vegetation present on the remote islands of KWNWR. As the captive colony (CC) of butterflies is intensively monitored, a direct adult count is recorded for every generation.

Figure 2

Census estimates of Miami blue butterfly adults at the south end colonies of Bahia Honda State Park, Florida from July 2002 to August 2008.

We calculated the harmonic mean across generational data available from BHSP and CC counts. The harmonic mean population size from census estimates for BHSP August 2005–May 2006 was compared to the temporal estimates of N calculated from September 2005 to June 2006 genetic data. For the captive colony, the harmonic mean was calculated from direct counts of adults from June 2006 to October 2006 and compared to temporal estimates of N calculated from July 2006 to November 2006 genetic samples (see Waples 2005).

Results

Population assumptions

Tests for linkage disequilibrium (followed by sequential Bonferroni correction) of the 12 microsatellite markers showed that loci were unlinked. Test results show seven of 12 loci in HWE in the two generations sampled from the BHSP population, although not the same seven loci. In the September 2005 population, loci CthB11, CthB101, CthB103, CthB115, CthC116, CthC127 and CthD7 were in HWE and in June 2006 loci CthB11, CthB103, CthB117, CthB119, CthC116, CthC124, CthC127 were in HWE. The KWNWR population had ten of 12 loci in HWE (CthB115 and CthB119 were not in HWE).

Null alleles were detected at two loci in BHSP September 2005 (CthB117 and CthB119) and at three loci (CthB119, CthC12, and CthC124) in BHSP June 2006 due to the presence of homozygote excess. Three loci showed evidence of null alleles in the KWNWR population (CthB103, CthB115, and CthB119). Micro-Checker (van Oosterhout et al. 2004) showed no evidence of large allele drop-out or scoring error due to stuttering. We further verified genotyping results by re-genotyping 5–10% of samples at nine loci. In loci where successful re-genotyping rates were less than 100%, samples were re-genotyped a third time to verify results. All statistical analyses were performed on this corrected dataset and homozygosity persisted after corrections. All private alleles fit the expected mutation pattern.

Overall, the SMM estimator R did not perform better than F in estimating population differentiation. Three of the 12 loci tested had an observed R greater than that estimated from random (pR, supplementary information). Results from the program Bottleneck (Cornuet and Luikart 1996) did not show significant evidence of a genetic bottleneck in the September 2005 or June 2006 BHSP samples analyzed. The Wilcoxon signed-rank tests (Luikart and Cornuet 1997) for these populations failed to show excess heterozygosity under the IAM (P = 0.62; P = 0.85 respectively) and TPM (P = 0.91; P = 0.95, respectively) models. Additionally, these two BHSP samples have an L-shaped allele frequency distribution that is characteristic of a population in mutation-drift equilibrium (Cornuet and Luikart 1996). Populations that have undergone a bottleneck show a mode shift in their allelic distribution. The KWNWR population had heterozygote excess (IAM, P = 0.003; TPM, P = 0.088) but showed no evidence of a mode shift. Only the captive colony population, which had undergone a known bottleneck, showed evidence of a mode shift. This population also showed heterozygosity excess under the IAM model (P = 0.003).

Estimates of N

Temporal estimates for BHSP obtained from NeEstimator and TempoFs provided comparable estimates of the harmonic N spanning the time period of August, 2005 through May, 2006 (N = 20 and 28 respectively) with considerable overlap in confidence intervals (Table 2). In contrast, the moment estimates of N generated in MNe was 136.44 and the likelihood estimate was 322 (95% CI 150–1000), with the upper estimate of 95% CI reaching the upper limit set for N (1000). The single-point estimators returned similar results to each other for both the September 2005 and June 2006 BHSP populations. The Bayesian estimator OneSamp (Table 3) varied modestly depending on the priors used although values consistently fell within the margin of estimates based on temporal samples. Linkage disequilibrium point estimates implemented in LDNe (Table 4) were similar to those from OneSamp with slightly lower estimates for BHSP September 2005, though with overlapping confidence intervals. Estimates decreased noticeably when low frequency alleles were omitted from LDNe analyses.

Table 2

Estimates of effective population size (N) of BHSP (September 2005–June 2006) using temporal methods

Program	N	Ne	95% CI
NeEstimator	–	20.9	9.8–62.8
TempoFs	19.7	28	17–83
MNe (moments-based)	–	136.44	–
MNe (likelihood)	–	322	150–1000

N denotes the population size input required under type I sampling for the program TempoFs, based on harmonic mean census results (see Table 1).

Table 3

Bayesian estimates of effective number of breeders calculated using the program OneSamp

Sample	Priors	Mean (SD)	95% lower	95% upper
BHSP Sept 2005	2-100	34.056 (3.748)	22.606	58.185
	4-200	27.775 (4.152)	17.970	47.028
	6-500	34.677 (12.064)	20.558	79.761
BHSP June 2006	2-100	40.701 (6.996)	24.516	68.864
	4-200	41.577 (16.651)	23.679	95.450
	6-500	29.260 (4.290)	20.730	57.572
KWNWR Feb 2008	2-100	28.721 (2.082)	19.920	46.448
	4-200	26.720 (2.179)	18.586	48.513
	6-500	24.327 (0.714)	17.397	42.247
CC July 2006	2-100	51.525 (4.112)	34.045	80.716
	4-200	32.095 (4.961)	21.984	56.410
	6-500	36.220 (1.092)	26.864	61.361
CC Oct 2006	2-100	16.454 (0.810)	11.898	23.530
	4-200	14.054 (0.599)	8.771	20.187
	6-500	16.376 (3.352)	9.957	28.952

The mean and standard deviation from three replicates is given and the 95% lower and upper values are the greatest and lowest of these replicates, respectively.

Table 4

Linkage disequilibrium estimates of effective number of breeders calculated using the program LDNe

Sample	# Independent comparisons	Lowest allele freq.	Ne (95% CI)
BHSP Sept 2005	399	0.05	12.7 (7.4–23.7)
	715	0.02	23.8 (14.2–49.5)
	715	0.01	23.8 (14.2–49.5)
BHSP June 2006	329	0.05	21.3 (12.7–39.8)
	640	0.02	35.9 (22.4–68.5)
	1153	0.01	46.2 (30.6–81.3)
KWNWR Feb 2008	674	0.05	19.2 (11.9–36.2)
	912	0.02	37.4 (20.8–106.9)
	920	0.01	38.2 (21.1–111.2)
CC July 2006	424	0.05	9.7 (6.4–14.2)
	657	0.02	13.3 (9.7–18.2)
	829	0.01	14.1 (10.6–18.8)
CC Oct 2006	243	0.05	106.8 (6.8–∞)
	243	0.02	106.8 (6.8–∞)
	243	0.01	106.8 (6.8–∞)

We also estimated N from the captive colony at F1 and F4 as a comparison with the natural populations. The temporal estimator of the harmonic mean of N between July and October 2006 of the captive colony was 37.2 (95% CI 14.4–244.8) using NeEstimator and 23 (95% CI 12–234) from TempoFs (Table 2). Single point estimates from the July 2006 captive colony ranged between 32 and 52 when estimated by OneSamp (Table 3). The October 2006 generation of the captive colony gave OneSamp estimates of 14–16. LDNe produced July 2006 estimates ranging from 10 to 14, and the October generation estimate was 107 across all frequencies of rare alleles (Table 4).

N/N ratios

We calculated the harmonic mean from BHSP census data between August 2005 and May 2006 as 19.7 resulting in a N/N ratio of ∼1 using NeEstimator and ∼ 1.4 using the TempoFs estimator. Both MNe estimates were exceedingly large and resulted in N/N ratios greater than 7. Point estimates using either Bayesian (OneSamp) or linkage disequilibrium (LDNe) methods produced N/N greater than 1.0 regardless of priors or allele frequencies (not shown). In contrast, N/N estimated from colony samples were 0.02–0.05 in July 2006 (based on June 2006 adult census) and less than 0.01–0.07 in October 2006 (based on September 2006 adult counts). Using the temporal estimators, N/N ranged from 0.05 to 0.09 for TempoFs and NeEstimator methods respectively. The latter was based on the harmonic mean count between June and October 2006 (433). Due to a number of concerns about the colony data (i.e. sample sizes and a demographic bottleneck) the N estimates and corresponding N/N values are estimated for qualitative purposes only (see below).

Discussion

Interpreting N in the two remaining wild butterfly populations

Our results suggest that N is very low in both KNWNR and BHSP populations. These results are consistent with the low census numbers that have been estimated since 2002 (Fig. 2). One concern with using N estimators on small populations has been the bias and low precision of estimators including those implemented in NeEstimator (Peel et al. 2004) and MNe (Wang 2005). Our results reflect these concerns given that the estimates from these two programs were far higher and had large variances. In contrast, the temporal estimator TempoFs, while still having high variance, produced an N value that was consistent with those estimated from point sample methods (Bayesian and gametic disequilibrium) and that was more realistic given census numbers in these populations. Although the temporal estimates reflect the harmonic mean N over the ∼8 generations rather than the effective number of breeders in the previous generation represented by point estimators, the long-term demographic pattern (Fig. 2) suggests that there has not been a major decline over the period dealt with in this paper, although population fluctuations are evident. Such fluctuations are known to have a significant effect on N, where the harmonic mean represents the demographic N. Hence, the temporal N and the point estimates reflect similar patterns. The consistent population average size is further supported by the nonsignificance of the bottleneck tests, suggesting that there has not likely been a recent (∼100 generations) dramatic decline in population size. As such, we consider the estimates produced by the point estimators used here and the TempoFs temporal estimator to be most accurate and likely provide important information on the genetic status of the remaining butterfly populations. Overall, N is critically low in the Miami blue butterfly and we feel that genetic stochasticity as well as demographic stochasticity are likely to be important factors determining the long-term persistence of these populations.

Properties of estimators

The temporal likelihood method of Wang (2001) as implemented in MNe is downwardly-biased when sample size is less than 50 (Palstra and Ruzzante 2008). Assuming a single isolated population with no immigration, the likelihood estimator as implemented in MNe produced the highest estimate of effective population size for the BHSP population, 322.81, and the moment-based method in MNe yielded N = 136.44, a value that is still several times higher than any of the other estimates. When sample sizes are small, likelihood methods may return biologically unreasonable estimates of N (Jorde and Ryman 2007). The other temporal estimators returned more consistent and reasonable results. A strength of the program TempoFs is that it weights alleles to reduce the sensitivity that affected the likelihood temporal estimator when polymorphism is highly skewed. All three temporal methods employed here returned high 95% upper CI. MNe hit the upper limit (1000) defined by us, and TempoFs remained higher than the point estimators for either BHSP temporal sample. NeEstimator returned an upper 95% CI most similar to the point methods and one most biologically reasonable (mean 21, 95% CI 10–63).

Results from this study show that single point estimations of N are consistent across different estimation methods (Bayesian in OneSamp versus LD in LDNe). The latter method is expected to be robust when sample sizes are small, however we encountered incongruent results when N = 10 (captive colony October 2006), suggesting a lower limit to the strength of this method in dealing with such small sample sizes. We also found that the Bayesian single-point estimators gave increased precision (95% CI ranges) as compared to the temporal likelihood method. However, when census size was large (as in the captive colony) the precision decreased in the LD method.

Interpreting N/N ratios

Care must be taken when comparing census data to N estimates because the information provided by the two numbers may not correspond temporally (Waples 2005). Our estimates of N were compared with the parental generation census for point estimators. For the temporal estimators we calculated the harmonic mean of census figures including the parental generation census numbers. Overall, N/N ratios were very high compared to most ratios for small populations. Data from BHSP reveal ratios at or above 1.0 (depending on N estimate used), unexpected values for threatened taxa. This may be a result of (i) imprecision in N, (ii) imprecision of census estimates, or (iii) potential impacts of low population size on reproductive variance. With respect to the first and second scenarios correcting these factors would not likely reduce the ratio to below 0.5 under realistic scenarios. First, some of the N estimators were selected because they are expected to perform better under smaller sample sizes by reducing the upward bias when the true N is greater than the sample size (Jorde and Ryman 2007; Waples and Do 2008). These estimators consistently gave smaller estimates than ‘biased’ estimators and therefore we assume that they have succeeded in reducing the bias associated with small sample sizes (a common consequence of working on endangered taxa). Survey estimates may be biased downward given the difficulty of observing a small insect (wingspan of 25 mm) under nonideal conditions in heterogeneous habitats. Given the potential difficulties in precise estimates an underestimate of N is likely although this bias would need to be very large to explain the N/N rations observed (e.g. counting 50% of existing butterflies). ‘Typical’ average N/N estimates have ranged from 0.1 to 0.5 (e.g., Nunney and Elam 1994; Frankham 1995; Palstra and Ruzzante 2008), although this varies considerably based on the ecology of the organism. Our wild populations have high N/N that, barring an upward bias in our N estimates or downward imprecision in census counts, suggests that the latter mechanism may be playing a role in producing high ratios. Miami blue butterfly habitat may be under stochastic influences which in turn should have an equal impact on all genotypes. This, at least in the short term, could upwardly bias N/N ratios. Similarly, highly reduced adult population sizes correspond to reduced variance in reproductive output and would result in more offspring from more reproductive pairings surviving to maturity. This form of ‘genetic compensation’ has been suggested for populations of salmonids (Ardren and Kapuscinski 2003; Araki et al. 2007) and damselflies (Watts et al. 2007), and may be an important genetic outcome of reduced population size. A further possibility is that all four of these factors lead to our results.

Census estimates of the BHSP population show that there are considerable population fluctuations over time. Environmental conditions like drought and rainfall affect host plant quality and flower nectar availability, both of which are required for large and healthy butterfly populations. Environmental stochasticity can also lead to the high N/N values because stochastic environmental effects should affect all individuals equally at neutral loci, reducing inter-individual reproductive variance (Palstra and Ruzzante 2008). This highlights the importance of gathering biological information and understanding the environmental and demographic history of a population when calculating N and N/N. Regardless of the specific demographic history of either population, severe fluctuations in population size can have major effects on N because of the relationship between the harmonic mean and variance in population size.

Data from the captive colony reveal lower ratios of 0.05–0.09; results that mirror other published estimates of N/N (Frankham 1995). However, our estimate is likely downward biased because of the extremely small sample size obtained for the October 2006 period. The knowledge of, and genetic evidence supporting, the presence of a demographic bottleneck further exacerbates our interpretation of the captive colony data. However, the rapid increase in colony size and relatively low N does suggest that although demographic numbers are large, which is the explicit goal of captive breeding for reintroduction, there has not been a concomitant increase in N. We hypothesize that the high N/N ratios in the two wild populations are the result of genetic compensation (above) and our understanding of inter-individual fecundity comes from the captive colony. Therefore the large variance among individuals may not be consistent in the wild at low population densities.

Conservation implications and applications

Two main concerns for estimating effective population sizes in threatened taxa are (i) accuracy and precision from small sample sizes and (ii) practical ability to sample threatened populations. In some species, repeated genetic sampling is not an option due to negative impacts of sampling, longevity of the organism, or research limitations to sampling. N estimators using temporal variance have been criticized because they assume all changes in allelic frequency are due to genetic drift and do not consider immigration from neighboring populations (Wang and Whitlock 2003). In the case of the Miami blue butterfly, as well as many other endangered taxa, there is no gene flow between populations and thus temporal estimators may be more useful (if multi-generational data are available). In their recent review, Palstra and Ruzzante (2008) reported a median unbiased estimate of N = 260 among 83 studies utilizing temporal methods, with estimates smaller for threatened taxa. The majority of our N estimates yield an effective population size of less than 50, regardless of sampling time or location. We recommend single-point estimations, especially implemented in LDNe for small effective population size estimation and Bayesian estimation as implemented in ONeSAMP when population samples are larger. Other studies (Watts et al. 2007) have found that precision is reduced in genetic estimates of N because of the short interval (one generation) between sampling periods. We have overcome this issue with eight generations between sampling periods but still encountered reduced precision using the temporal likelihood estimator. If temporal estimates are employed by a conservation program, they should be accompanied by point estimates of each sampled generation as well. The inclusion of multiple estimates will likely reduce the upper 95% CI and provide a more robust and useful estimate of the effective population size. Additionally, differences between point and temporal estimates may highlight issues such as bottleneck events or population crashes that may have occurred during the sampling period.

As N is often a criterion for listing endangered species (Mace et al. 2008), a precise estimate of N is often important for determining the legal status of imperiled taxa. Therefore, knowledge of the precision and limitations of the various N estimators will help inform decision-making and aid in the estimation of a biologically-accurate effective population size. Imprecise measurements may have conservation consequences such as premature delisting, down listing, or removal from protected status. Caution should therefore be exercised when interpreting N estimates from genetic data and both the estimation method itself, sample size, and polymorphisms of marker loci need to be considered. The time-scale over which data were collected and the population structure (if present) should also be considered when interpreting N estimates. Based on the N criterion alone, this butterfly meets IUCN standards for a critically endangered taxon. Moreover, there is a growing number of captive-breeding programs (IUCN 2008; Mace et al. 2008) and these provide the additional opportunity of comparing captive and wild N estimates. We recommend that other programs continue this effort as an additional way to assess the genetic health of their captive-bred populations.

As an extinction risk, a greater short-term genetic problem to the remaining two natural populations of Miami blue butterflies may be the lack of gene flow rather than the current N. Even limited gene flow can offset inbreeding depression (Vila et al. 2003; Hogg et al. 2006). Given the highly fragmented and reduced range of the Miami blue, efforts should concentrate on restoring gene flow through the reestablishment of populations within the current dispersal range of existing populations with ex situ propagated individuals. The conservation decision to augment populations should not be made without careful consideration of habitat availability and the potential for the introduction of maladapted genotypes. In the meantime, the severely reduced size of the existing populations suggests that genetic factors along with environmental stochasticity may already be affecting the persistence of the Miami blue butterflies.