We interpreted the results of nuclear DNA sequencing to be inconsistent with the recognition of California Gnatcatcher (Polioptila californica) subspecies. McCormack and Maley (2015) suggested that our data did support 2 taxa, one of which was P. c. californica, listed as Threatened under the Endangered Species Act (ESA). We summarize here how 2 sets of researchers with access to the same data reached different conclusions by including different analyses. We included the southern subspecies' boundary from the taxonomy of Atwood (1991), the taxonomic basis for the ESA listing, which resulted in an Analysis of Molecular Variance that provided no support for subspecies. In contrast, using a novel taxonomic hypothesis without precedent in the literature, McCormack and Maley (2015) found statistically significant FST values for 2 loci, which they suggested supports P. c. californica. We propose that our mitochondrial and nuclear data had sufficient power to capture geographical structure at either the phylogenetic (monophyly) or traditional “75% rule” level. McCormack and Maley (2015) suggested that finding an absence of population structure was a “negative result,” whereas we consider it to be the null hypothesis for a species with gene flow and no geographical barriers. We interpret the unstructured mtDNA and nuclear DNA trees, the STRUCTURE analysis supporting one group, the identification of just 26% (and not 75%) of individuals of P. c. californica with the most diagnostic nuclear locus, the overall GST that suggests that over 98% of the variation is explained by nontaxonomic sources, and the lack of evidence of ecological differentiation to indicate that P. c. californica is not a valid subspecies. McCormack and Maley (2015) suggest that statistically significant differences at 2 loci that explained <6% of the genetic variation, and previous morphological data, support recognition of P. c. californica. If ornithology continues to recognize subspecies, these different standards should be reconciled.
Subspecies have long been a controversial taxonomic rank in general (Wilson and Brown 1953), as well as in ornithology (Remsen 2005, Zink 2015). Part of the issue stems from a discrepancy between how subspecies were described historically and some of their modern uses. Many avian subspecies were described from few specimens and few localities, and represented geographic variation in only one or at most a few characters (Zink and Remsen 1986). Consequently, subspecific nomenclature sometimes reflects arbitrary or subjective divisions of single character clines (Barrowclough 1982). If the goal of an investigator is simply to find examples of geographic variation, most subspecies could be useful indicators of potential study systems. Alternatively, if the use of a subspecies requires it to be an evolutionarily significant unit (Moritz 1994), such as in a phylogenetic, comparative, or biogeographical study, many described avian subspecies are inappropriate (Zink 2004, Phillimore and Owens 2006).
Thus, the definition of subspecies is an important issue. Barrowclough (1982) suggested that subspecies should be held to the same standard as other taxa in the Linnaean hierarchy, i.e. that their taxonomic names should be predictive of concordant character variation. Cracraft et al. (1998) suggested that subspecific taxa should be diagnosable. Ornithology has a history (Amadon 1949) of attempting to apply a “75% rule,” in which “75% of a population effectively must lie outside 99% of the range of other populations for a given defining character or set of characters” (Patten and Unitt 2002:27). The lack of a consistent standard for recognition has led to a state in which currently described subspecies range from ahistorical entities to discrete evolutionarily significant units.
Efforts to conserve and preserve biodiversity often take formal subspecies taxonomies at face value, with the assumption that described subspecies represent discrete and evolutionarily independent taxa. It is the case, however, that the assumption cannot be made that all described avian subspecies qualify as discrete evolutionary units (Remsen 2005). We suggest that if a subspecies is selected for listing under the U.S. Endangered Species Act (ESA), then at least one major character system, such as genetics, morphology, ecology, or behavior, should provide diagnostic support (Moritz 1994), especially given the economic costs of species preservation (McCarthy et al. 2012). We further suggest that statistical differences in morphological characters or gene frequencies alone might not support the existence of a subspecies.
The California Gnatcatcher (Polioptila californica) is distributed linearly from southern California, USA, to the tip of Baja California Sur, Mexico. It has been divided into subspecies, although the number of subspecies and their geographical distributions differ among authors (Miller et al. 1957, Atwood 1991, Mellink and Rea 1994). The subspecies scheme of Atwood (1991) was used to frame the listing of the northern subspecies, P. c. californica, as Threatened under the ESA (USFWS 1993). At the species level, Birdlife International (2015) categorizes the California Gnatcatcher as a species of Least Concern. Zink et al. (2000) found that there was no significant geographical differentiation in mitochondrial DNA (mtDNA) within the entire species, and therefore that there was no support for any previously described subspecies, nor any other geographical divisions previously unrecognized taxonomically.
The U.S. Fish and Wildlife Service (USFWS 2011) concluded that mtDNA data were insufficient for testing taxonomic boundaries in the California Gnatcatcher, although they have relied on this data for many prior listing decisions in the past 20 yr. Furthermore, several bird (Zink et al. 2001) and other vertebrate species (Riddle et al. 2000) show distinct mtDNA breaks over the same range as that of the California Gnatcatcher, revealing the potential for mtDNA to detect differentiation if it exists. The rationale for not relying on mtDNA alone is that, in spite of the mtDNA genome including many genes, mtDNA genes are linked without recombination, and whether one analyzes one or many genes, the result is a single gene tree. A single gene tree can be a weak test, irrespective of whether it is from the nuclear or mitochondrial genome. Thus, many researchers have turned to analyzing genetic variation across multiple nuclear loci. Zink et al. (2013) followed USFWS (2011) guidance and evaluated nuclear loci to test the mtDNA results from the California Gnatcatcher. Zink et al. (2013) concluded that the nuclear DNA results were consistent with the mtDNA results, and showed that there was no significant geographical differentiation in the California Gnatcatcher that would support taxonomic recognition.
In contrast, McCormack and Maley (2015) suggested that some of the genetic data published by Zink et al. (2013) supported subspecies designation for the threatened subspecies of California Gnatcatcher (P. c. californica). Below we discuss how 2 different sets of researchers reached different conclusions from the same data.
Choice of loci
McCormack and Maley (2015) expressed concern over the way in which nuclear loci were surveyed by Zink et al. (2013). An important issue in assessing patterns of genetic variation is whether the events of interest, such as isolation of populations, occurred within a time frame relevant to the resolving capability of the marker of choice. It is well established that the time to coalescence is proportional to the effective size of the locus under study. For mtDNA, the effective size is ¼ that of an autosomal locus, which means that, on average, any mtDNA gene tree will coalesce 4 times more quickly than any random nuclear gene tree. Therefore, mtDNA has been very important in discovering recently isolated groups of populations (Avise 2009), many of which are of interest to conservation, whereas most nuclear loci will not distinguish these populations. Given that the mtDNA tree for California Gnatcatchers revealed no geographic structure (Zink et al. 2000), it follows that there would be a low probability of nuclear loci recovering structure missed by the mtDNA analysis (Zink and Barrowclough 2008, Barrowclough and Zink 2009). However, a few instances exist for birds in which mtDNA failed to capture significant evolutionary divergence (McKay and Zink 2010, Toews and Brelsford 2012) that was revealed by nuclear loci. The question then becomes how to analyze genetic variation across nuclear loci. Some investigators have suggested that microsatellite loci “evolve” more rapidly than other single-locus nuclear genes, but this is incorrect. Although microsatellite loci selected for use in population genetic studies have high mutation rates, and often high numbers of alleles, the coalescence time of a locus is independent of mutation rate, depending instead on the effective size. Among nuclear loci, only sex-linked (Z-linked in birds) loci evolve (coalesce) somewhat faster because they have an effective size ¾ that of an autosomal locus (including microsatellite loci, unless they are sex-linked). Furthermore, it is not straightforward to jointly analyze microsatellite allele frequencies that represent populations and mtDNA sequences from individuals. Hence, many researchers (e.g., Brito and Edwards 2009) have advocated sequencing nuclear genes to complement mtDNA gene trees, which, following direction from the USFWS (2011) to assess variation in nuclear genes, was the approach followed by Zink et al. (2013). There are newer genomic techniques that allow for detection of single-nucleotide polymorphisms from thousands of loci, and we have used these to assess geographic variation in the California Gnatcatcher, but we do not discuss these results here.
Across the 7 variable loci, 52 alleles were observed, with an average of 7.4 alleles per variable locus. This level of variation is more than sufficient to detect significant geographic structure if it exists. Although Zink et al. (2013) reported the overall amount of genetic variation distributed among populations across all loci, i.e. GST = 0.013 (i.e., 1.3% of the total), as well as the locus-by-locus values, they did not perform an Analysis of Molecular Variance (AMOVA) for each locus. In this case, an AMOVA evaluates the distribution of genetic variation, with hierarchical levels including individuals within populations, populations within subspecies, and subspecies relative to the total genetic variance. The AMOVA results of McCormack and Maley (2015) comprise the core of their reanalysis of the data of Zink et al. (2013); we assume that they did not simply pool all individuals into 2 groups, which would not constitute an AMOVA. McCormack and Maley (2015) expressed concern that Zink et al. (2013) did not employ an AMOVA approach; however, Zink et al. (2013) cited a previous paper (Zink 2010) that pointed out the potential arbitrariness of such results (see below).
McCormack and Maley (2015) concluded that the FST values for 2 loci (ACON1-I15 and TGFB2-I5) were statistically significant when comparing P. c. californica with the rest of the range. However, this statistic was driven by an excess of rare alleles as a result of larger sample sizes in the north (see below), as well as by population expansion (Zink et al. 2013). We are mostly concerned about the way in which McCormack and Maley (2015) constructed their tests. McCormack and Maley (2015:382) stated that “For test 1, we separated the recognized subspecies californica (Los Angeles south to San Telmo) from more southern populations (Misión San Fernando south to Cabo San Lucas), based on Atwood's (1991) quantitative morphological subspecies boundaries and the USFWS initial listing boundary of 30°N latitude.” This misrepresents the subspecies scheme of Atwood (1991; see McCormack and Maley  figure 1), which includes 3 subspecies, not 2. We find it inappropriate to accept Atwood's (1991) northern subspecies boundary and then to exclude the southern boundary and pool the 2 southernmost subspecies, because in effect this manipulates sample sizes and could lead to a spurious result. McCormack and Maley (2015) therefore tested a taxonomic hypothesis, without justification, that has not been previously proposed.
We divided population samples from Zink et al. (2013) into 3 groups to match the subspecific taxonomy of Atwood (1991): (1) Los Angeles, Riverside, Orange County, San Diego, Ensenada, and San Telmo; (2) Misión San Fernando, El Rosarito, San Ignacio, Mulegé, and Villa Insurgentes; and (3) La Paz and Cabo San Lucas; and we used the program Arlequin 188.8.131.52 (Excoffier and Lischer 2010) to compute AMOVAs. A significant among-group variance component would be considered to indicate population or subspecies structure. If a significant value were to be found, it would not, however, identify which particular group or groups were different. In addition, the pattern across all loci should be examined, and whether the data support other groupings not encompassed by prior subspecies boundaries should be determined.
In our hierarchical AMOVA, the FST of 0.034 for TGFB2-I5 was not significant (P = 0.13), which suggests that McCormack and Maley's (2015) result for that locus was an artifact of using a novel taxonomic hypothesis. For ACON1-I15, the overall FST of 0.053 was statistically significant (P = 0.034). However, this masked the actual pattern of divergence. Computing pairwise AMOVAs for the ACON1-I15 locus, we found that P. c. californica was significantly different from P. c. margaritae (FST = 0.056, P = 0.048), the subspecies immediately to the south, but P. c. californica was not different from the southernmost subspecies, P. c. abbreviata (FST = 0.034, P = 0.20), nor were P. c. abbreviata and P. c. margaritae different (FST = 0.077, P = 0.07). Thus, following the logic of McCormack and Maley (2015), the samples from La Paz and Cabo San Lucas at the southernmost extent of the distribution, representing P. c. abbreviata, would have to be pooled with those from northern P. c. californica, creating a leapfrog taxon, which would be inconsistent with traditional taxonomic schemes (including that of Atwood ). To explore heuristically the effect of relatively increased sample sizes in the north, we took a random sample of 8 alleles from each of the larger Los Angeles and Orange County samples to standardize them to the average sample size for the other localities. The resulting AMOVA returned an insignificant FST (0.061, P = 0.08). Although the magnitude of the FST was the same as that derived from the analysis of all individuals, it showed that assessing statistical significance was dependent on uneven sample sizes in the northern samples representing P. c. californica. Thus, the ACON1-I15 locus also does not support the evolutionary distinctiveness of P. c. californica. The statistical significance at 2 loci found by McCormack and Maley (2015) was an artifact resulting from uneven sample sizes and the use of a taxonomic hypothesis that was not used in framing the ESA listing decision. Furthermore, AMOVAs and other allele frequency–based tests assume that populations are at equilibrium. The nuclear DNA data suggested that California Gnatcatcher populations were expanding (Zink et al. 2013), hence not at equilibrium (Slatkin 1993), which can affect interpretation of statistical significance (Boileau et al. 1992).
There is an additional concern around the strategy employed by McCormack and Maley (2015), identified by Zink (2010). If we followed McCormack and Maley's (2015) 2nd subspecies scheme (from Mellink and Rea ), which we argue is irrelevant for conservation because the ESA listing is based on Atwood (1991), and divided our northern samples into 2 new sets (Los Angeles to San Diego and Ensenada to San Telmo), FST increases to 8.7% of the variance (P = 0.003). If we create a new subspecies for Los Angeles and another one for Riverside County, FST becomes 11% with a P-value of ~0. What we have done, in our opinion, is use the small samples for a single locus to promote sampling error into biologically meaningless statistical significance. In fact, with just a small number of loci, one could adjust hierarchical schemes in such a way that most populations could merit subspecies recognition. In our view, this is not the best way to test the genetic distinctiveness of subspecies, especially those listed under the ESA.
A major advantage of using molecular characters to study geographic variation is the conceptual link between allele frequencies and population genetics theory. Our results (figure 3 in Zink et al. 2013) indicated that genetic differentiation did not increase significantly with geographic distance; Rousset (1997) described the causal relationship between this lack of differentiation and high levels of gene flow. Contra McCormack and Maley (2015), we tested the relationship with the appropriate Mantel permutation test (Dietz 1983) using a Spearman rho-statistic, for the 9 populations with average samples sizes ≥5. Our analysis suggested that gene flow was the process resulting in the lack of monophyletic groups, or any significant geographic differentiation, within the California Gnatcatcher.
Given historical precedence, we evaluated the “75% rule” with our genetic data. We used the data for ACON1-I15 and TGFB2-I5I5 and computed how many individuals in P. c. californica could be uniquely diagnosed using unique alleles for either locus. For ACON1-I15, 26% of the individuals in P. c. californica could be identified by the alleles they possessed, and for TGFB2-I5 only 12%. Therefore, the genetic data that McCormack and Maley (2015) suggest support the taxonomic and evolutionary distinctiveness of P. c. californica identify at best 26% of individuals, ⅓ of the traditional “75%” rule. This reveals why, in our opinion, statistically significant but small values of FST are not biologically useful for taxonomic and conservation decisions (see Björklund and Bergek 2009).
A disagreement between Zink et al. (2013) and McCormack and Maley (2015) concerns the way in which genetic data are used to test subspecies limits. We believe that to test subspecies boundaries, locus-by-locus analyses are suboptimal because the actual question is whether there are 3 geographically and evolutionarily independent lineages. Hence, we used BEAST (Drummond and Rambaut 2007) to combine all loci (excluding bFib-I7) into a lineage tree (see Appendix for details). The multilocus tree (Figure 1) does not support Atwood's (1991) scheme, nor any other geographically structured taxonomic hypothesis. That is, not only are none of the subspecies monophyletic, but only 2 pairs of geographically adjacent populations (out of 12 possible) are each other's closest genetic relative. The posterior density of allele relationships plotted as a cloudogram (Appendix Figure 2) further exhibits the lack of population structure expected at the earliest stages of differentiation. Cloudograms, unlike traditional evolutionary trees, allow for the visualization of conflicts between tree topologies, instead of showing a single consensus tree. Briefly, with little to no conflict, cloudograms show topological congruence in the shape of well-delineated clusters of branching clouds of the different trees due to similar patterns of allele sorting between loci. On the other hand, a lack of branch congruence is visualized as a web reflecting incomplete lineage sorting, as is the case in the California Gnatcatcher.
Lastly, McCormack and Maley (2015) excluded mention of our STRUCTURE analysis based on all loci jointly (figure 4 in Zink et al. 2013), which also shows but a single genetic group across the range. That is, a standard multilocus analysis failed to return evidence of multiple groups that could correspond to subspecies, and, for whatever reason, McCormack and Maley (2015) decided not to mention it in their critique. In our opinion, corroboration of subspecies limits comes from analyses of all loci simultaneously, not from locus-by-locus searches designed to find any statistical support for a particular (and not necessarily optimal) taxonomic scheme.
We find it surprising that McCormack and Maley (2015) disagree with our interpretation of our test for niche divergence, given that we followed the method used by McCormack et al. (2010), and more recently by Jezkova et al. (2015), who used the same climatic variables for the same biogeographical region. In short, we did not find significant niche divergence in the California Gnatcatcher, subject to the caveats relevant to all such niche studies that the results depend on the environmental (climatic) layers and sampling sites used to build the model. Most widespread species span variation in climatic space. Simply being widespread and occupying differing habitats does not mean that there are evolved niche differences. In the case of the California Gnatcatcher, observing that northern populations exist in a more mesic environment than populations to the south does not necessarily mean that they are ecologically differentiated, and the niche background tests suggest that the species has a wide ecological tolerance. Indeed, separating these 2 alternatives was the point of the original development of the background test (Warren et al. 2010). However, niche tests are in an early stage of development, and it is possible that considering other dimensions of the niche would produce a different result (Soberón 2007).
Given the lack of genetic support, in our opinion, for P. c. californica, the question arises whether morphological data support subspecies limits. This is not the place for a thorough review of morphology, and neither Zink et al. (2013) nor McCormack and Maley (2015) performed any morphological analyses. However, McCormack and Maley (2015) suggested that prior morphological work supported the validity of P. c. californica. We present a brief overview to illustrate past work that we suggest provides tenuous support for subspecies.
McCormack and Maley (2015:380–381) stated that there exists “…a century of prior work documenting the occurrence of a distinct population of gnatcatchers in southern California based on evidence of physical differences (summarized in Mellink and Rea 1994).” Furthermore, McCormack and Maley (2015:384) stated that the standard for recognition of subspecies is that character “differences must not vary smoothly from one population to another, but instead must show a discontinuity (i.e. step cline) in their character values. Over the years, many studies on the California Gnatcatcher have delimited changes in phenotypic characters consistent with discrete variation, affirming the distinctness of P. c. californica from subspecies to the south (Grinnell 1926, Van Rossem 1931, Phillips 1980, Atwood 1991).” We find these assertions to be misleading.
The original subspecies descriptions were based on 3 specimens (P. c. californica; Brewster 1881), 2 specimens (P. c. margaritae; Ridgway 1904), and 9 specimens (P. c. abbreviata; Grinnell 1926), with no quantitative assessment of geographical patterns of variation, which was standard practice in the early years of avian subspecies description. Subsequent studies (Grinnell 1926, Van Rossem 1931, Phillips 1980) also provided no evidence for step clines. More recently, Atwood (1991) studied geographical variation in several characters based on study skins and found a color character that appeared to support P. c. californica, although it excluded a sampling location at the southern extent of the range of the subspecies. However, J. L. Atwood later stated that he would not consider that character valid because of the inadequacy of the spectrophotometer that he used at the time ( http://www.pacificlegal.org/document.doc?id=1477). Thus, although Atwood's (1991) unrooted UPGMA (Unweighted Pair Group Method with Arithmetic Mean) phenogram appears to support geographic differentiation, it includes the same discredited color character(s).
The suggestion of McCormack and Maley (2015) that Mellink and Rea (1994) provided diagnostic support for a new subspecies is difficult to understand, given that Mellink and Rea (1994) did not report the results of a single statistical or quantitative test. They compared plumage colors using visual inspection and comparison with color charts, a method criticized by other traditional taxonomists (e.g., Browning 1993). Hence, the results of Mellink and Rea (1994) are not consistent with a standard of concordant step clines. Although McCormack and Maley (2015:384) note that Mellink and Rea's (1994) subspecies “atwoodi” is accepted by some “authoritative taxonomic references,” Remsen (2005:409) noted that these references are “only cursory summaries of those largely qualitative taxonomic judgments.” We stand by our interpretation of Skalski et al. (2008) that the morphological data support gradual or conflicting clines instead of concordant step clines, a conclusion that echoes that of Cronin (1997). We do agree, however, with the final point of McCormack and Maley (2015) that a comprehensive phenotypic analysis, using modern methods and all available specimens and controlling for potential sources of error such as specimen age, would be useful.
McCormack and Maley (2015) concluded that Zink et al. (2013) had a conflict of interest (“sponsorship bias” in their wording; p. 386) based on a quote in an article in the L.A. Times ( http://www.latimes.com/science/la-me-gnatcatcher-20140630-story.html), wherein one of us (R. M. Zink) appeared to suggest that our 2013 study was funded in part by “land developers.” In fact, that comment was meant to refer to an earlier published study, Zink et al. (2000), that was also referenced in the L.A. Times article. Zink et al. (2000:1403) acknowledged financial support from “the U.S. Navy, the National Fish and Wildlife Foundation, Southern California Edison, trustees of the Manomet Center for Conservation Sciences, the Building Industry Association of Southern California, the Transportation Corridor Agency, Chevron Land and Development, the University of Minnesota, and the National Science Foundation (DEB-9317945).” Thus, the quote from the L.A. Times article was in no way intended to reflect that the research reported in our 2013 paper was funded or influenced by “land developers.”
To address McCormack and Maley's (2015) concern about conflict of interest for the 2013 work, we note that a lawyer, Mr. Robert Thornton, approached the authors in 2006 and stated that he represented clients who wished to sponsor the research on geographic variation in nuclear genes that the USFWS (2011) stated was required to test subspecies limits in its published denial of a previous listing petition. The resulting contract, which supported a laboratory technician and supplies for 9 mo, listed no funding source other than Mr. Thornton's firm, and, in fact, the authors were unaware of the identity of the funders, as reflected in the acknowledgments of the 2013 publication. However, given the concerns of McCormack and Maley (2015), we queried Mr. Thornton as to the funding source behind the 2006 contract, and he revealed (R. Thornton personal communication) that the funds were provided by the Transportation Corridor Agency, a public agency formed by the legislature of the State of California. In retrospect, we could have discovered and disclosed the funding source for the 2013 paper, although it would not have influenced our analyses and writing of the paper in any way.
We understand that our failure to discover and disclose the fact that funding came from the Transportation Corridor Agency created a conflict of interest because Mr. Thornton has provided legal counsel in opposition to listing the California Gnatcatcher in 1994, has represented various land developers in southern California that have contested the listing of the California Gnatcatcher, and was one of 2 lawyers representing 7 plaintiffs in a petition to remove the California Gnatcatcher from the list of threatened species under the Endangered Species Act (filed June 11, 2014; http://www.pacificlegal.org/document.doc?id=1477). However, McCormack and Maley (2015) imply that this conflict of interest might have influenced our interpretation of our data. We note that Mr. Thornton was not an author on the 2013 paper, and that neither Mr. Thornton nor the California state agency were given any opportunity to comment on our study design (e.g., choice of loci), analyses, interpretation of results, and design and preparation of the resulting manuscript. Furthermore, the University of Minnesota, through which the funds were disbursed, has strict rules that prohibit funding agencies from having any influence on the interpretation of results and subsequent publication of research. Thus, despite any apparent connections between 3rd parties with vested interests and our research, the conclusions in Zink et al. (2013) flowed exclusively from the data that we collected and analyzed to address the subspecific status of P. c. californica.
McCormack and Maley (2015) stated that if the California Gnatcatcher were to be delisted, 197,000 acres (797.2 km2) would potentially be subject to development, rendering it unsuitable for the gnatcatcher. Of course it is possible that mitigation would be able to secure additional suitable land elsewhere. However, the words “acre” or “acreage” do not occur in Zink et al. (2013), and our paper had nothing to do with this issue, only with the scientific issue of whether the subspecies was distinct. In addition, it is unclear to us that a particular acreage makes an ESA listing decision more or less important, although its geographic location is clearly a factor.
Investigator philosophy and subspecies
McCormack and Maley (2015) implied that because one of the authors (R. M. Zink) of Zink et al. (2013) has repeatedly criticized the subspecies category, it was inappropriate for the paper to address the status of P. c. californica. First, the paper was also authored by G. F. Barrowclough, who long ago (Barrowclough 1982:602) wrote that “most subspecies are not to be taken too seriously.” Secondly, Zink (2004:563) wrote: “Only taxa defined by the congruence of multiple morphological or molecular characters should be recognized at some rank. Over 90% of continental avian subspecies fail this test.” That is, Zink (2004) suggested that some subspecies are valid taxa, although whether they are ranked as species or subspecies is debatable. In addition, the studies of Zink et al. (1995, 1997, 2001) have supported the validity of many taxa currently ranked as subspecies. Just because a scientist has a history of finding little support for a concept, be it subspecies, molecular neutrality, or competitive exclusion, it does not disqualify him or her from conducting legitimate scientific research in these areas of inquiry.
The subspecies concept has long been debated among ornithologists with regard to its value in classifying taxonomic variation within species, and divergent opinions are commonplace. In the original Linnaean system of classification, all units of one taxon were more similar to each other than they were to members of other taxa. With the discovery of phylogenetic methods in the mid-20th century, clustering into taxa based on overall phenotypic similarity was replaced by the concept of genetic relatedness and hierarchical monophyly. That is, all members of a species share a common ancestor and a more recent evolutionary history with each other than with members of any other species. This standard now holds at all levels of the taxonomic hierarchy except for subspecies. If subspecies were diagnosable, a trait shared by all individuals in that subspecies would reveal their monophyletic origin. However, avian subspecies are currently described on the basis of some less-than-universally shared trait and are defined on the basis of geography, i.e. all individuals breeding in a region are members of the same subspecies, whether or not their genetically closest relatives occur there. This lack of agreement on basic Linnaean standards for subspecies leads to disagreements over the validity of subspecies as taxa.
At any level of geographic or taxonomic organization, gnatcatcher populations and subspecies are not monophyletic, and, for example, some individuals from California have a more recent history of ancestry with gnatcatchers in Cabo San Lucas (Baja California Sur) than they do with other individuals from California. Consequently, if subspecies need not be monophyletic, then depending on the trait used to define them, membership is possible in multiple, overlapping, and nonhierarchical groups. Thus, as Wilson and Brown (1953) pointed out, a subspecies based on color could easily conflict with an alternative subspecies based on wing length. This is the Achilles Heel of many subspecies, including the California Gnatcatcher.
Negative results or different null hypotheses?
McCormack and Maley (2015) built their argument around the notion that the results of Zink et al. (2013) were “negative.” Most philosophers of science would take issue with this; one does not prove something is true, one falsifies hypotheses. If a series of populations, such as those analyzed for the California Gnatcatcher, is connected by gene flow, as our earlier mtDNA data showed (Zink et al. 2000), the null hypothesis from population genetics is that the populations will show no geographic structure. Therefore, finding a result of no differentiation is not “negative”; rather, our data reflected the expectation of the null hypothesis and thereby falsified a hypothesis of the existence of any monophyletic units within the gnatcatcher (i.e. of any hierarchical Linnaean taxon), or of any unit consistent with a traditional 75% rule for subspecies. McCormack and Maley (2015) adopted the null hypothesis that the subspecies are morphologically distinct, which in our opinion (see above) is based on tenuous grounds, but nonetheless which our data also falsified. In fact, our mtDNA and nuclear DNA data did have sufficient power to find evolutionarily distinct lineages, such as those that exist in other species whose ranges span the distribution of the California Gnatcatcher (Zink et al. 1997, 2001, Vázquez Miranda 2014); however, neither mtDNA nor nuclear DNA analyses supported geographical divisions.
Biology, politics, and conservation
We believe that conservation scientists should interpret objectively and impartially all of the best available scientific and commercial data (the ESA standard) to retain credibility in the public realm. We have based our conclusions on several types of data and analyses. We disagree that the morphological basis for the subspecies is well established. We have shown that the mtDNA gene tree, the nuclear gene lineage tree, and our STRUCTURE analysis do not support any subspecies or geographical groupings of the California Gnatcatcher. Furthermore, only 26% of individuals of P. c. californica were correctly identified by the most diagnostic locus and, across all nuclear loci, only 1.3% of the total genetic variation was “explained” by geography. The niche models that we used did not show evidence of niche divergence. These observations reveal to us that P. c. californica is not a valid taxon. In contrast, McCormack and Maley (2015) believe that the statistical significance of FST values explaining 6% of the genetic variance at 2 loci between subspecies (albeit from a nonexistent taxonomic hypothesis) and prior morphological work is sufficient to retain P. c. californica as a subspecies. In our opinion, McCormack and Maley's (2015) stance seems more aligned with a societal or political goal rather than a biological one, namely their laudable concern that the delisting of P. c. californica could lead to further loss of native habitat. The scientific question, however, is not about potential acreage lost, or a general debate about subspecies, it concerns only whether P. c. californica is a valid taxon. In our opinion, all of the best available data do not support the validity of P. c. californica or any other California Gnatcatcher subspecies.
Our interpretation of the data does not mean that we see no value in preserving remaining tracts of the gnatcatcher's habitat, coastal sage scrub. We suggest that to maintain both scientific credibility and to aid the preservation of coastal sage scrub it is important to interpret all of the analyses and data, rather than to exclude relevant analyses. Thus, environmentalists need to concentrate on other species or other ways to preserve this habitat, rather than risking the erosion of scientific credibility by attempting to defend the validity of P. c. californica. In many ways, the fact that there is no accepted standard for the recognition of subspecies is the root of the differences between Zink et al. (2013) and McCormack and Maley (2015). We suggest that such debates will continue unless a consensus is reached on minimal standards for subspecies recognition.
We thank S. M. Lanyon, J. Cracraft, J. Griffin, R. Gutiérrez, two anonymous reviewers, and K. Bostwick for useful comments on the manuscript.
Funding statement: See ‘Funding' section above. The authors received no further funding for this work. None of our funders had any influence on the content of the submitted or published manuscript, and none of our funders required approval of the final manuscript to be published.
Ethics statement: Our research followed our respective institutional guidelines.
Author contributions: R.M.Z. and G.F.B. conceived the idea and design; R.M.Z., G.F.B., J.F.F., and H.V.M. conducted the research; R.M.Z., G.F.B., and H.V.M. wrote the paper; G.F.B., J.G.G., and H.V.M. designed the methods; and G.F.B., H.V.M., and R.M.Z. analyzed the data.
Methods Used to Create a Multilocus Lineage Tree for California Gnatcatchers
The early stages of speciation are characterized by pervasive incomplete lineage sorting. Classic phylogenetic and probabilistic algorithms for individual clustering for population delimitation often do not account for this phenomenon. We tested for the possibility of discrete clusters of populations—“subspecies” sensu McCormack and Maley (2015)—in the presence of incomplete lineage sorting by employing a full multispecies coalescent approach (Heled and Drummond 2010) in BEAST 1.8 (Drummond and Rambaut 2007). We mapped alleles onto the species tree based on the sampled populations in Zink et al. (2013) and included 1 Black-tailed Gnatcatcher (P. melanura) as an outgroup. In taxa in which discrete units exist despite a shallow evolutionary scale, populations can be used as terminals for species tree estimation (Smith et al. 2014). We included all data used by Zink et al. (2013) with the same settings as for their extended Bayesian skyline plot, but excluded the locus beta-Fibrinogen Intron 7 (see McCormack and Maley ). Briefly, species-tree and gene-tree estimation were achieved with 3 combined runs of 108 generations, sampling every 1,000 genealogies, with a speciation Yule tree prior, and a strict clock model (see Zink et al.  for justification). Our estimation of the species and population genealogy is the maximum clade credibility tree after a 10% burn-in. Single “summary” trees, though easy to visualize, do not capture the complexity of the Bayesian tree probability space exploration and the density distribution of allele sorting. Moreover, to portray the rampant amount of allele sharing at the earliest stages of population divergence and the uncertainty of gene tree estimation, we plotted trees in the 95% posterior clade credibility distribution as a cloudogram in DensiTree 2.0 (Bouckaert 2010) using consensus tree plots that reflected the frequency of a given topology by its intensity. Cloudograms allow the finding of dominant topologies and delimiting of clades—if groupings do exist—despite the state of allele sorting in a given group (Bouckaert 2010).