PLANS TO LIST populations and species under the U.S. Endangered Species Act (ESA) or under similar legislation in other countries can be informed by several types of information. Morphological studies, especially those focused on subspecies taxonomy, suggest hypotheses for evaluating whether different parts of a species' range are evolutionarily independent. Modern genetic studies can reveal whether populations or subspecies are evolutionarily discrete or significant in relation to the species as a whole and document levels and pathways of historical or ongoing gene flow. Furthermore, genetic data can show whether particular populations are lacking in genetic variability owing to inbreeding or bottlenecks and whether populations have undergone recent range and population increases. Various genetic markers exist, each suited to answering different questions about threatened species (Brito and Edwards 2009). For example, recent restrictions to gene flow might best be revealed by analyzing microsatellite allele frequencies (Barr et al. 2008). Major evolutionary divisions are perhaps best identified by analysis of mitochondrial DNA (mtDNA) or nuclear DNA sequences (Zink and Barrowclough 2008).

Although it has received less attention than morphological and genetic data, quantitative tests of niche divergence can show whether a population is ecologically distinct. Threatened populations that occur in different habitats present a concern for conservation biologists because it must be determined whether such populations harbor evolved ecological adaptations or simply reflect a generalist ecological strategy. Short of long-term, costly, and logistically difficult cross-fostering or common garden experiments, ecological niche modeling provides a basis for making quantitative assessments of ecological differentiation in a hypothesis-testing framework (Warren et al. 2010). If populations from different areas have significantly different niche characteristics, one might conclude that the populations qualify as a distinct population segment under the ESA (U.S. Fish and Wildlife Service [USFWS] and National Marine Fisheries Service 1996), irrespective of whether they were also genetically or morphologically discrete. Information from multiple sources provides a more complete perspective on threatened and endangered species in the context of their preservation (May et al. 2011).

We assessed the genetic and ecological distinctiveness of the subspecies of the California Gnatcatcher (Polioptila californica). The California Gnatcatcher is a small nonmigratory songbird that ranges from Los Angeles, California, to the southern tip of the Baja California peninsula (Atwood and Bontrager 2001). The species occupies coastal sage scrub (CSS) in the northern portion of its range; to the south of the CSS, the habitat occupied is typical of the drier Sonoran Desert and population density is much higher. Three independent taxonomic studies based on morphology supported a subspecies in the CSS, although the proposals differed in the geographic locations of taxonomic boundaries (Miller et al. 1957, Atwood 1991, Mellink and Rea 1994). The northern, coastal subspecies (P. c. californica) is listed as threatened under the ESA (USFWS 1993) because of habitat loss and fragmentation and concomitant population declines. The California Gnatcatcher has served as a flagship species for preservation of the highly fragmented CSS ecosystem (Atwood 1993). The enormous economic value of real estate in this area of already high and increasing human population density has led to scrutiny of the subspecies status of the coastal California Gnatcatcher and whether it should be listed under the ESA (Cronin 1997). Cronin (1997:663) stated that “I believe that subspecies designations of P. californica should be ignored until thorough phylogenetic analyses of geneticallybased characters are done.” The USFWS (2011:66255) concluded “that the coastal California Gnatcatcher constitutes a valid subspecies….”

Comparison of mtDNA sequences obtained throughout the range of the California Gnatcatcher revealed lower genetic diversity in the CSS but did “not support any subspecies scheme, either previously described [Fig. 1] or unforeseen” (Zinket al. 2000:1398). In the present study, we integrate new nuclear gene-sequence data with past studies of mtDNA (Zink et al. 2000) and morphology (Atwood 1991). In addition, we test for significant niche divergence (McCormack et al. 2010) as a proxy for ecological distinctiveness. Our goal is to obtain a more complete understanding of the species' recent history and to provide a perspective on its conservation from multiple sources of data.

METHODS

Samples.—Samples (Table 1) were those used in Zink et al. (2000), although degradation of the DNA reduced the sample size. Some analyses were conducted on divisions of the sample localities into northern and southern, based on the pattern of nucleotide diversity found in Zink et al. (2000), with the dividing point at 28°N, between El Rosarito and San Ignacio (see Fig. 2). Exact sample sizes per locus per locality can be obtained from the GenBank accessions (KC863990-864745).

Genetic analyses.—We sequenced seven introns and one exon (Table 2) using primers from Backström et al. (2008) and Kimball et al. (2009). In addition, as a check on the mtDNA control-region sequences, and for calibrating population size changes, we sequenced 43 individuals for the mitochondrial gene ND2 (23 individuals from the north, 20 from the south, with one sequence of the sister taxon, P. melanura, as an outgroup); these data were analyzed with ARLEQUIN (Excoffier et al. 2005). From initial sequencing runs, we determined whether each nuclear sequence included two or more heterozygous sites; all such sequences were phased into component alleles using specially constructed polymerasechain-reaction primers that would anneal with only one of the two alternate bases at the most 5′ heterozygous site. The resulting sequence allowed us to determine the sequence of one of the two alleles, and the other allele was inferred by subtracting this allele from the ambiguous sites. For each of the seven introns and one exon, we computed the number of gene copies sequenced, the number of sequenced base pairs, the number of alleles observed, nucleotide diversity for each population sample, and the G_ST statistic of Holsinger and Mason-Gamer (1996). We also computed the mean sample size and nucleotide diversity, with bootstrap confidence intervals, over loci for each population. We computed G_ST across loci among all populations and among populations with sample sizes >4, along with bootstrap confidence intervals. These statistics were computed using software written by G.F.B.

For each locus, we used a four-gamete test to determine the longest fragment of DNA consistent with no recombination (Hudson and Kaplan 1985), using PAUP*, version 4.0b10 (Swofford 1998). We constructed minimum spanning networks of alleles for those fragments using PAUP* and computed the frequency of each of the alleles at each sampling locality to show the geographic deployment of genetic variation. We used the Hudson-KreitmanAguade (HKA) test (Hudson et al. 1987) implemented in the HKA program (see Acknowledgments; Hey 2004) to test for departures from neutrality in the ND2 and all nuclear loci and alleles simultaneously (using P. melanura as an outgroup) by employing 10,000 coalescent simulations to assess statistical significance.

FIG. 1.

Examples of patterns of geographic variation at three nuclear loci and the mitochondrial DNA control region of California Gnatcatchers (from Zink et al. 2000). Arrows indicate the root of the network derived from phylogenetic analysis using P. melanura as theoutgroup. Darker shading indicates the range of the California Gnatcatcher. Names of sample locations are given in Figure 2 and Zink et al. (2000). Numbers of individuals per sample are given in Table 1 and in CenBank entries. The dotted line (far right) indicates the division between northern and southern localities based on the pattern of nucleotide diversity discovered in mtDNA.

We used the program STRUCTURE (Pritchard et al. 2000) to determine the number of groups based on the nuclear DNA alleles, as determined above. In brief, this program considers variation among loci simultaneously and attempts to determine whether there is more than one genetically distinct group represented in the sample. The seven variable nuclear loci were formatted as single-nucleotide polymorphisms (SNP; see Manthey et al. 2011). The STRUCTURE analyses assumed an admixture model, correlated allele frequencies, and a fixed lambda value (which was inferred by setting K = 1 and allowing lambda to be estimated in an initial analysis). We analyzed the data for K = 1 to 6 with five replicates for each value of K. Each run contained 100,000 steps as burn-in, followed by 500,000 steps. The δK was calculated ad hoc (Evanno et al. 2005) and used to identify the best estimate of K.

TABLE 1.

Nucleotide diversity (π) in California Gnatcatchers over seven nuclear loci. Localities are shown in Figures 1 and 2.

We used the extended Bayesian skyline plot (EBSP) method (Heled and Drummond 2008) implemented in BEAST, version 1.7.4 (Drummond et al. 2012), to estimate changes in population size through time with multiple loci in a coalescence-based framework (Heled and Drummond 2008). We applied the EBSP analysis to the nuclear plus mitochondrial ND2 sequences for individuals identified as northern or southern by Zink et al. (2000). The ND2 gene was set as the reference locus to calibrate time and population size with a substitution rate of 0.0125 site^-1 lineage^-1 Ma^-1 (Smith and Klicka 2010) using a strict clock prior. By definition, strict clock rates do not incorporate error rates or confidence intervals (adding deviation estimates would convert such a rate into a relaxed clock, even for small ranges) and are the appropriate rates for the estimation of divergences in data sets with low variation (Brown and Yang 2011). The generation time was set at 1 year for small passerine birds. Although we acknowledge that 1 year is an estimate for this species in the absence of a detailed life table, juvenile males and females could start breeding in their first year (Atwood and Bontrager 2001). The EBSP analyses included a Markov chain Monte Carlo (MCMC) run of 100 million steps, sampled each 1,000 steps with a strict clock prior for the mitochondrial locus and relaxed exponential priors for the nuclear loci to account for rate variation among different loci, and a linear demographic model. Trace plots were checked using TRACER, version 1.5 (Rambaut and Drummond 2007), to assess convergence in MCMC analyses. If the 95% high posterior density (HPD) of the estimate of the number of size-change steps (the parameter “demographic.populationSizeChanges”) excluded zero, we concluded that a significant change in population size occurred (Lim and Sheldon 2011). To evaluate the contribution of the mtDNA data to the overall results, we repeated the EBSP analysis using only nuclear loci, to check for population increases at different relative times (e.g., excluding a temporal calibration).

FIG. 2.

Geographic distribution of nucleotide diversity (π) in mtDNA and seven pooled nuclear loci (nuDNA) in California Gnatcatchers. MtDNA data (top panel) suggest that populations north of 28δN were genetically less variable than populations to the south of the line, where mtDNA nucleotide diversity shifted from high (south) to low. The mtDNA pattern was not observed among the nuclear loci.

TABLE 2.

Gene locus characteristics in a population survey of California Gnatcatchers. ^a Note that MC1R is an exon, whereas all other loci are introns.

Ecological analyses.—The occurrence of coastal California Gnatcatchers in the mesic CSS could indicate ecological differentiation sufficient for recognition as a distinct population segment. We constructed correlative ecological niche models (ENM; Peterson et al. 1999, Elith et al. 2011) using breeding records from the Breeding Bird Survey and the Global Biodiversity Information Facility (see Acknowledgments), which were input into MAXENT, version 3.2.2 (Phillips et al. 2006). We divided locality points into those representing CSS (n = 144) and southern populations (n = 48). Climatic data (19 layers) were obtained from the Worldclim bioclimatic database (Hijmans et al. 2005). Each ENM was based on the average of five MAXENT runs and plotted using DIVA-GIS, version 7.1.7.2 (Hijmans et al. 2004). Distribution maps were coded so that predicted presence or absence was based on the logistic threshold for equal training and testing specificity produced by MAXENT. We used the niche identity test and the background test for niche divergence in ENM Tools (Warren et al. 2010); random localities used in the test (Warren et al. 2008) were obtained from Hawth's Tools (Beyer 2004) in ARCGIS, version 9. We show results based on Schoener's D (Schoener 1968), which ranges from 0 to 1 (identical niches). We also estimated the distribution of suitable habitat for the California Gnatcatcher at the Last Glacial Maximum (LGM; 21,000 years ago) using all 192 distribution points.

RESULTS

Seven of the eight loci surveyed were variable, although each locus lacked a structured geographic pattern (Fig. 1 and Table 1), as did the ND2 sequences (F_ST = 0.021, P = 0.37). The locus MC1R has been associated with darker-colored phenotypes in some organisms (Baião et al. 2007); this locus also lacked geographic structure in California Gnatcatchers, despite the CSS populations being characterized as having somewhat darker plumage (Atwood 1991). Nucleotide diversity (π) did not differ across localities, unlike the result of the mtDNA study (Fig. 2) and for the new ND2 sequences (π = 0.0013 for north, 0.0030 for south). The pie charts in Figure 1 suggest a greater diversity of rare alleles in the north; however, the average sample size across loci was 67.9 in the north and 36.5 in the south (Table 1), which suggests a sample-size explanation for the apparent higher number of rare alleles in the north. The G_ST values across loci ranged from -0.089 to 0.244 (Table 1), and the overall G_ST was 0.013 (95% bootstrap confidence interval: -0.058 to 0.104), indicatinga negligible level of genetic divergence across localities. A similar result was obtained for mtDNA (control region) data (Zink et al. 2000). Plots of genetic distance versus geographic distance suggested a lack of isolation by distance for mtDNA and autosomal loci (Fig. 3), further supporting an inference of no geographic structure. The HKA test was insignificant (χ² = 8.14, df = 7, P = 0.32), indicating that selection did not influence the pattern of genetic variation.

FIG. 3.

Relationship in California Gnatcatchers between genetic and (log) geographic distance, showing modest isolation by distance (P = 0.37 for nuclear DNA; P = 0.08 for mtDNA; Spearman's rank correlation). This plot also reveals a lack of geographic breaks, which would appear as displaced groups of points.

Analyses of nuclear loci combined from STRUCTURE (Fig. 4) identified no geographic groupings that corresponded with any previously suggested subspecies, nor any other significant evolutionary divisions. The δK statistic was nonsensical because the most likely value of K was 1. Thus, the nuclear gene-sequence data do not support individual population assignment or clustering to named subspecies or any other geographic groupings.

The EBSP analyses based on all loci (Fig. 5) suggested that the southern population underwent an increase that began ∼3,500 years ago, whereas the northern population increased to a lesser degree, and more recently. Because the timing of the increases depended on the calibration used, we do not interpret the dates precisely but find only that there was at least one expansion in the south and, potentially, a more recent one in the north. Using only nuclear loci (not shown), the results for both northern and southern populations had HPD confidence intervals for the demographicsize-changes parameter that overlapped zero, and we therefore could not infer a population increase, although the most frequent size-change factor found by the Markov chain in the south was 1.

Niche models (Fig. 6) for the two groups of populations did not cross-predict the majority of each other's ranges, although each predicted co-occurrence around 28–30°N, where the more mesic CSS meets the drier southern vegetation. The observed niche identity (0.343) was less than the distribution of randomized values (Fig. 7), which suggests that the niches differ more than expected by chance. Background tests suggest that the niches are not significantly divergent (Fig. 8). The estimated distribution of California Gnatcatchers in the LGM (Fig. 9) corresponded to the southern half of the present-day distribution and involved considerable area to the west of the current coastline that was above water at that time.

FIG. 4.

Results of STRUCTURE analysis for seven variable nuclear loci showing a lack of geographic structure and assuming that K = 2 (i.e., two groups of California Gnatcatchers). The same result was obtained when K = 3. The highest probability of the data occurred when K = 1. Numbers across the bottom correspond to the samples from north (left) to south in Figure 2. The scale on the y-axis is percentage.

DISCUSSION

Genetic results and comments on molecular markers.—Some authors have suggested that natural selection biases the pattern of variation in mtDNA genomes (Ballard and Whitlock 2004, Galtier et al. 2009:4546, Balloux 2010), but the question is whether selection is sufficiently strong to obscure evolutionary patterns and processes, which several studies suggest is not the case for birds (Zink 2005, Zink et al. 2006, Hung et al. 2012). Our HKA test revealed no evidence of strong selection that would bias mitochondrial or nuclear loci. We therefore believe that the loci we analyzed provide a basis for understanding the recent demographic and evolutionary history of the California Gnatcatcher.

FIG. 5.

Bayesian skyline plots for northern (gray lines) and southern (black lines) populations of California Gnatcatchers (dashed lines indicate 95% HPD intervals). The southern population underwent an expansion beginning ∼3,500 years ago (or earlier, depending on the calibration used), whereas the northern population increased to a lesser degree and more recently. The closeness of the 95% HPD to zero for the northern sample suggests caution in inferring a population increase.

Our assessment of genetic population structure of the California Gnatcatcher was the same for mtDNA and nuclear loci, both of which suggest that no evolutionarily significant divisions exist within the species. Thus, the mtDNA data provided a proper assessment of genetic structure. In general, relatively few avian subspecies qualify as valid evolutionary entities (Zink 2004, Phillimore and Owens 2006). We recommend that an mtDNA survey should be part of attempts to determine whether a species includes multiple lineages, owing to the rapid coalescence time for mtDNA gene trees. However, any single gene tree could misrepresent the lineage tree (Toews and Brelsford 2012), and, hence, both phylogeography and conservation genetics have become multilocus efforts. The question is what type of nuclear gene information should be used in concert with mtDNA (recognizing that all nuclear loci, including microsatellites, have longer coalescence times than mtDNA). We advocate the use of sequences from nuclear loci because we think that the potential to compare coalescence analyses based on both nuclear and mtDNA sequences provides a sounder basis for making evolutionary inferences than comparison of mtDNA sequences with microsatellite allele frequencies (Brito and Edwards 2009, Zink 2010).

FIG. 6.

Niche models showing predicted occurrence of California Gnatcatchers in coastal sage scrub (gray) and southern populations (black); hatched area indicates predicted overlap.

FIG. 7.

Distribution of 100 random values of Schoener's D for the samples from the coastal sage scrub and the south. Arrow shows observed Schoener's D value, which indicates that the two groups of California Gnatcatchers do not have identical niches.

Populations north of 28°N harbored less mtDNA variability than those to the south (Zink et al. 2000). Furthermore, several southern mtDNA haplotypes rooted basally on an otherwise unstructured haplotype tree, also suggesting a southern refugium. The geographic break in level of variability coincides with a biogeographic split in many species, potentially owing to a midpeninsular seaway (Lindell et al. 2006, Leaché et al. 2007). Zink et al. (2000) suggested that California Gnatcatchers recently invaded the northern part of their current range, which could explain the lower nucleotide diversity in the north, owing to leading edge expansion (Hewitt 2000). This is consistent with the estimated distribution of California Gnatcatchers at the LGM (Fig. 9), which suggests that they were limited to the southern half of the current range. In contrast to the mtDNA data, we found no difference in nucleotide diversity at nuclear loci (Fig. 2) between southern and northern parts of the range. We suggest that because of the larger effective population size of nuclear loci, the signature of a northward expansion might be “erased” more quickly at nuclear loci than for mtDNA, owing to the fact that dispersing individuals carry two copies of nuclear genes, whereas only females carry a haploid mitochondrial genome.

FIG. 8.

Results of test for niche divergence between coastal sage scrub and southern populations of California Gnatcatchers. Gray bars indicate the distribution of random samples from the coastal sage scrub population, and black bars indicate those for the southern population. Arrow indicates observed Schoener's D value for the two populations, indicating neither niche divergence nor conservatism.

FIG. 9.

Predicted distribution of California Gnatcatchers at Last Glacial Maximum based on locality points from coastal sage scrub and southern parts of the current range. Considerable habitat existed offshore of western Baja California Sur, which is now under water. Although the analysis shows the occurrence of suitable niche space on the west coast of Mexico, the species currently does not occur there.

Mismatch distributions based on mtDNA control-region data (Zink et al. 2000) and ND2 (not shown) suggested that populations increased at different times, with a southern expansion preceding a northern one. Our EBSP plots (Fig. 5) based on all loci are consistent with this finding and suggest that the population expansion began not immediately at the end of the LGM but perhaps as recently as 3,000 years ago. We did not find a similar signal in EBSP analyses based on nuclear data alone, especially for the northern population. This suggests that it is important to include mtDNA in such analyses because it likely provides more mutations in the relevant time frame of population expansion (Keinan and Clark 2012).

Morphological support for subspecies.—It is possible that neither mtDNA nor nuclear loci coalesce rapidly enough to capture recent geographic isolation, owing to the lag time inherent in all molecular markers (Zink and Barrowclough 2008). If polygenic morphological characters were under selection and they possessed more additive genetic variance than a single locus (such as mtDNA), they might provide support for subspecies that originated in less time than required for reciprocal monophyly at mtDNA (i.e., 2Ne_f generations). Evidence for evolutionarily significant morphological units would include concordant character step-clines, consistent with a morphology-wide response to historical isolating events (Barrowclough 1982). The molecular results suggest scrutiny of the morphological basis of the subspecies upon which the USFWS based its opinion.

Reanalysis of the morphological data (Atwood 1991) on California Gnatcatchers led Skalski et al. (2008) to conclude that the coastal California Gnatcatcher was “incorrectly listed under the ESA due to misinterpretation of morphological data.” Instead of a concordant pattern of discrete character gaps across common geographic localities, these authors found a pattern of gradual morphological change across the range, with inconsistent patterns among characters. This geographic pattern of morphological variation is consistent with genetic data, which suggests that the subspecies were arbitrary divisions of idiosyncratic morphological gradients, and not equivalent to discrete evolutionary (listable) entities.

Ecological distinctiveness.—If the sole criterion for ecological distinctiveness is successful persistence in two or more environments, then the gnatcatchers could qualify as two DPSs (Fig. 6). However, we suggest that the intent of this criterion is to protect populations that have differentially adapted to novel environments. Populations that occupy different environments across a species' range could simply represent a species with broad ecological tolerance, and not indicate that each population possesses evolved ecological adaptations to differing habitats. For example, McCormack et al. (2010:1231) discovered that although allopatric populations of jays in the genus Aphelocoma occupied different habitats, niches were not significantly different because “the allopatric environments they occupy are not significantly more divergent than expected under a null model.” We found a similar result in California Gnatcatchers (Fig. 8). Although the species occupies the distinctive CSS habitat in the north (Fig. 7), the two groups do not exhibit significant niche divergence if the backgrounds of each environment are taken into account. In other words, the species appears to be a habitat generalist. To our knowledge, this is the first attempt to quantify the ecological distinctiveness criterion of the ESA using the niche background test.

As a caveat, we add that the methods for testing niche divergence are in a relatively early stage and that the test is only as good as the models and input data. For example, we tested for differences in what has been termed “Grinnellian” niche dimensions (Soberón 2007), and inclusion of other types of variables could yield a different result.

Conservation implications.—The U.S. Congress directed the USFWS to list taxa under the ESA on the basis of the best available scientific and commercial data. In denying a petition to delist the California Gnatcatcher, the USFWS (2011:66258) relied on statements such as these:

[T]he argument that the California Gnatcatcher is not distinct from other populations is based on a single genetic character, mtDNA, and this is a far too narrow and limited technique for making determinations of taxonomic validity…. [The mtDNA data set alone] does not present substantial information that the current subspecies taxonomic classification of the coastal California Gnatcatcher may be in error.

These arguments could be construed to mean that mtDNA does not qualify as the best available evidence if used alone. One of the litmus tests for assessing the scientific value of a particular data set is its degree of congruence with other data sets. We show that results from morphology, mtDNA, nuclear DNA, and ecological niche modeling agree that the California Gnatcatcher is not divisible into discrete, listable units. For listing, we believe that at least one data set should provide clear evidence of distinctiveness and significance. Furthermore, the U.S. Congress advised that the ESA should be applied “sparingly” when listing DPSs (Bernhardt 2008). Given that the coastal California Gnatcatcher lacks morphological, genetic, and ecological significance, it becomes difficult to justify its listing. Although it is difficult to prove a hypothesis of no divergence, the congruence of evidence suggests that this is the most strongly supported hypothesis given the best available data. This outcome is unfortunate, in that preservation of CSS has benefited from listing of this subspecies. Our analysis refocuses attention on the importance of ecosystem-wide preservation.