Evaluating the taxonomic status of the Great White Heron (Ardea herodias occidentalis) using morphological, behavioral and genetic evidence

Abstract

The Great White Heron (GWH) has an all-white plumage and occurs in the Gulf of Mexico and Caribbean. Described originally as Ardea occidentalis, it is now considered a subspecies of Great Blue Heron (GBH; A. herodias). GWH and GBH meet in Florida Bay at the southern tip of Florida, providing the opportunity to evaluate their interaction and species status. To this end, we examined size variation and mate choice across their contact zone and genetic variation range-wide. Measurements of 7 morphological characters indicate trends, but not a significant difference, in size between GBH and GWH in southern Florida. GBH and GWH nest mainly in different places (mainland vs. islands) and at different peak times. In Florida Bay, mixed pairs occur, but white-white and blue-blue pairs are more common than in a randomly mating population. Assessing mating, however, is complicated because most, if not all, nesting blue birds are of mixed parentage. Microsatellite DNA analysis indicates that white and blue herons in Florida Bay and the outer Keys (outside Florida Bay) form a group distinct from blue forms on Florida Peninsula and elsewhere in North America. However, some gene flow occurs from white herons on the outer Keys to white and blue herons in Florida Bay, and from blue herons in Florida Bay to GBH on the Florida Peninsula. GWH alleles occur in all North American populations, but whether this is from gene flow or incomplete lineage sorting is unknown. Deciding GWH's species status is difficult. GWH and GBH meet in an ecotone where some gene flow occurs, but behavior and habitat largely isolate them. We argue in favor of splitting GWH from GBH. Regardless of how it is ultimately classified, the GWH's small population needs to be actively managed as an isolate in an extremely vulnerable environment.

RESUMEN

Ardea herodias occidentalis (garza blanca) tiene un plumaje completamente blanco y se encuentra presente en el Golfo de México y el Caribe. Descripta originalmente como A. occidentalis, actualmente es considerada como una subespecie de A. herodias (garza azulada). A. h. occidentalis y A. herodias se encuentran juntas en la Bahía de Florida en el extremo sur de Florida, brindando la oportunidad de evaluar su interacción y estatus de especie. Con este propósito, evaluamos la variación en tamaño y la elección de pareja a lo largo de la zona de contacto y la variación genética en todo el rango. Las mediciones de siete caracteres morfológicos indican tendencias, pero no una diferencia significativa, en tamaño entre A. h. occidentalis y A. herodias en el sur de Florida. A. h. occidentalis y A. herodias anidan principalmente en diferentes lugares (continente versus islas) y en diferentes momentos pico. En la Bahía de Florida, existe la presencia de parejas mixtas, pero las parejas blanca/blanca y azul/azul son más comunes que en una población reproductiva al azar. La evaluación del apareamiento, sin embargo, es complicada debido a que la mayoría, sino todas, las aves azules que anidan tienen un parentesco mixto. El análisis de ADN microsatelital indica que A. h. occidentalis y A. herodias en la Bahía de Florida y en los Callos externos (afuera de la Bahía de Florida) forman un grupo distintivo de las formas azules de la Península de Florida y de otros lugares de América del Norte. Sin embargo, algo de flujo génico ocurre desde A. h. occidentalis en los Callos externos hacia A. h. occidentalis y A. herodias en la Bahía de Florida, y desde A. herodias en la Bahía de Florida hacia A. h. occidentalis en la Península de Florida. Los alelos de A. h. occidentalis están presentes en todas las poblaciones de América del Norte, pero se desconoce si esto se debe a flujo génico o a una separación incompleta de linaje. La decisión del estatus de especie de A. h. occidentalis es difícil. A. h. occidentalis y A. herodias se juntan en un ecotono donde ocurre algo de flujo génico, pero el comportamiento y el hábitat las aísla en buena medida. Nos manifestamos a favor de separar A. h. occidentalis de A. herodias. Independientemente de cómo se clasifiquen en última instancia, la pequeña población de A. h. occidentalis necesita ser manejada activamente como algo asilado en un ambiente extremadamente vulnerable.

INTRODUCTION

The taxonomic relationship between the Great White Heron (Ardea herodias occidentalis) and the Great Blue Heron (A. herodias) presents interesting intellectual and practical challenges. Currently, the Great White Heron (GWH) is treated as one of 4 recognized subspecies of the Great Blue Heron (GBH) in North America (Eisenmann et al. 1973, Dickerman 2004, Dickinson and Remsen 2013, Clements et al. 2017). Unlike other GBH subspecies, which have a “blue” plumage and are widely distributed throughout North America, GWH has an all-white plumage and is restricted mainly to southernmost Florida and the Florida Keys, where the breeding population reaches a maximum of ~1,300 pairs (Powell and Bjork 1996, Hunter et al. 2006). In the rest of the Caribbean and Gulf of Mexico, the GWH is rare (Stevenson and Anderson 1994, Raffaele et al. 1998, Vennesland and Butler 2011), and its distribution and breeding are poorly documented. Because of its concentration in the Florida Keys, the GWH population is vulnerable to natural catastrophic events (e.g., hurricanes), habitat degradation from development and introduced predators (e.g., Burmese python [Python bivittatus]), and thus is of conservation concern. However, despite its distinctiveness and vulnerability, the GWH receives remarkably little attention from conservationists. As a formally recognized subspecies it merits as much conservation attention as a full species (Haig et al. 2006, Phillimore and Owens 2006, Winker 2010), but GWH seems to be largely neglected, perhaps because of the impression that it is simply a variant of the common GBH and occurs regularly throughout the Caribbean and Gulf of Mexico.

GWH was treated as a separate species from GBH by the American Ornithologist’s Union (Holt 1928, AOU 1957) until 1973, when it was changed to the rank of subspecies by Eisenmann et al. (1973), who cited Mayr (1956), Meyerriecks (1957) and Bond (1961) for support. Mayr (1956) summarized the evidence for including GWH with GBH: (1) the occurrence of mixed blue and white breeding pairs in the Florida Keys; (2) the occurrence of an intermediate plumaged form, Wurdemann’s Heron (previously A. h. würdemanni), in Florida Bay, indicating successful interbreeding between the blue and white birds; (3) the presence of pure blue GBH in Florida Bay during the breeding season (although Mayr admitted that such records were not backed by specimens); and (4) dimorphic populations containing both blue and white individuals (purportedly) in Cuba, The Isle of Pines, Jamaica and the Yucatán. Meyerriecks (1957) interpreted the occurrence of blue and white breeding pairs in the Florida Keys as providing support for random mating, but cautioned against accepting this hypothesis without further study because he had documented the plumage color of both members of mated pairs at only 9 nests. Bond (1956) described the Caribbean subspecies of GBH as dimorphic, and Bond (1961) noted a mixed colony of “gray” and white herons along the Venezuela coast, as well as reports by fishermen “of both gray and white young” occupying the same nest in that colony (but see the next paragraph). Bond also gave the opinion that GWH originated in the Caribbean and constituted an Antillean element in the Florida Keys.

Thus, original support for GWH as a subspecies of GBH stemmed mainly from the impression that interbreeding blue and white birds occur throughout the Caribbean and that mate choice between these birds is random with respect to plumage. However, although GWH may occasionally occur in GBH breeding colonies, e.g., on mainland North America (Bancroft 1969, McHenry and Dyes 1983), there is no evidence that major populations of GBH and GWH meet and interbreed anywhere except in Florida Bay. Bond’s (1961) reference to GWH mixing with “gray” individuals on the coast of Venezuela would, if true, refer to Cocoi Heron (A. cocoi), an allospecies of GBH, but the interbreeding of GWH and Cocoi Heron has never been verified. Moreover, records of GWH-GBH interbreeding elsewhere in the Caribbean are unsupported by publications or specimens. In fact, our understanding of GWH and GBH occurrence outside of the Florida Keys is muddled. In addition to the vague historical assertions of Mayr (1956) and Bond (1961), one reason may be that ornithologists and birdwatchers visiting the Caribbean nowadays do not distinguish between GWH and GBH, calling both GBH on lists, thus diminishing current records of GWH distribution. Another reason is that migrating (non-breeding) GBH which reach the Caribbean may be incorrectly considered to be breeding inhabitants, sympatric with GWH. The result of these problems is that we have a remarkably poor understanding of the dispersion and nesting of GWH outside of the Florida Keys.

The contact between GWH and GBH in Florida Bay provides an opportunity to study the degree to which these 2 taxa are reproductively isolated. Under a Biological Species Concept (BSC) framework, a key question for taxa in direct contact is whether the 2 mate randomly or assortatively (Mayr 1963, Price 2008). For GWH and GBH, the question remains open: neither random mating within the Florida Keys population nor breeding between blue and white forms elsewhere has been established. Despite a lively historical debate regarding the GWH’s taxonomic status (Holt 1928, Mayr 1956, Meyerriecks 1957, Lazell 1989, Stevenson and Anderson 1994), remarkably little attention has been paid to the interactions between sympatric white and blue herons. In addition, variation in their morphology has barely been assessed—only one morphological character (occipital plume) has been measured in blue and white herons in Florida Bay (Holt 1928)—and no genetic comparisons have been published.

To help determine the extent of reproductive isolation between GBH and GWH, we collected and compared relevant morphological, behavioral and genetic data. Our objectives were to clarify the taxonomic status of GWH and to provide information that may guide conservation efforts. In addition to coloration (white vs. blue), we have examined size variation, mating preference, mating localities and gene flow across the GWH/GBH contact zone in southernmost Florida and the Florida Keys.

Background: Morphology

In a comprehensive revision of the subspecies of blue GBH populations in North America, Dickerman (2004) examined plumage and size variation. He found size to be an unreliable taxonomic indicator, although he did use it to some degree in making decisions. Mostly, he based his revision on plumage, reducing the number of GBH subspecies from a high of 9 (Oberholser 1912) to 4 (Figure 1). Of these 4, 3 have a similar “blue” plumage: (1) A. h. fannini in the Pacific northwest; (2) A. h. herodias in the north from eastern Washington to the southeastern coast; and (3) A. h. wardi south of the other 2 subspecies and including the Florida Peninsula. Dickerman (2004) did not examine GWH per se, noting only that color distinguished it as a subspecies (A. h. occidentalis) from other GBH populations.

FIGURE 1.

Distribution of the 4 North American Great Blue Heron subspecies that were sampled to evaluate geographic variation. Red dots mark genetic sampling localities (listed in Table S2). See Figure 3 for details of area D.

Of these subspecies, only A. h. occidentalis contains individuals with 2 plumage types: the all-white form and a rarer blue form, Wurdemann’s Heron, which is a variable intermediate between all-white and typical GBH plumages. Formerly treated as a separate subspecies, A. h. würdemanni, Wurdemann’s Heron occurs entirely within the range of A. h. occidentalis in Florida Bay and is generally presumed to represent a hybrid between A. h. occidentalis and A. h. wardi. Its existence has caused a century-long debate about the degree of gene flow between GWH and GBH. As part of this debate, taxonomists investigated size differences among 3 subspecies of GBH that occur in the east (Ridgway 1878, Oberholser 1912, Holt 1928, Mayr 1956), and Zachow (1983) documented a significant difference: A. h. herodias was smallest, A. h. wardi was intermediate in size and A. h. occidentalis was largest. However, only Holt (1928) compared white and blue herons within the Florida Keys’ breeding population, and he measured only a single character, occipital plume length. Among females, he found plumes of A. h. wardi to be longest, those of Wurdemann’s Herons to be intermediate and those of white individuals to be shortest. Among males, A. h. wardi has longer plumes than Wurdemann’s and white individuals; the latter 2 do not differ in plume length.

White vs. dark morphs are relatively common in herons. Mock (1978) identified 6 species of dichromatic herons (including GBH) in which adults are either white or dark, e.g., Reddish Egret (Egretta rufescens) and Western Reef Heron (E. gularis). A seventh species, Gray Heron (A. cinerea), is dark throughout most of its range (Europe, Asia and Africa) but includes an isolated population of nearly white individuals (A. c. monicae) on small islands off the coast of Mauritania. Ardea herodias and A. cinerea are usually treated as closely related members of a superspecies (Bock 1956, Dickinson and Remsen 2013). They both have continent-wide distributions and, with 2 exceptions (A. h. occidentalis and A. c. monicae), are monochromatic throughout their ranges. These exceptions may be the result of similar evolutionary events, i.e. geographic isolation followed by drift or adaptation to local conditions causing phenotypic divergence. With respect to adaptation, whiteness in dimorphic herons has been shown to confer a feeding advantage in bright light (Tickell 2003, Green and Leberg 2005), an environmental situation expected in coastal areas of tropical islands. However, researchers have documented that GWHs feed primarily at night (Powell and Powell 1986, Powell 1987). So, if white plumage is adaptive, some other selective force (e.g., heat tolerance) is likely to be at play.

GBH differs from other dichromatic herons in at least 2 respects. First, although the ratio of white to dark individuals varies among populations in most dichromatic species, both color phases are generally present within any given population. However, white individuals are rare or unknown in mainland GBH populations, and when they occur they are generally assumed to be vagrant GWH (Mitra and Fritz 2002) or leucistic individuals. Second, other dichromatic species contain few intermediates; individuals are usually either all white or uniformly dark. The Florida Keys population exhibits a wider range of intermediate plumages than in any other dichromatic heron species (Mock 1978); indeed, most breeding blue-colored herons in Florida Bay have intermediate (Wurdemann’s) plumage (Holt 1928, Mayr 1956). The concentration and substantial variation of blue phenotypes in Florida Bay is more consistent with the hypotheses that this area is: (1) a contact zone between 2 previously isolated taxa, and the intermediates are hybrids with a variable degree of backcrossing; or (2) a cline across a freshwater/saltwater ecotone; than (3) a site encompassing a dichromatic subspecies of the GBH.

Background: Breeding

Approximately 800–1,300 pairs of GWH breed in the shallow marine and coastal mangrove environment of Florida Bay and the Florida Keys (Powell and Bjork 1996, Hunter et al. 2006). Florida Bay is a large shallow estuary that receives freshwater from the Florida Everglades. It is open to the Gulf of Mexico on its western boundary, and lies between the Florida Everglades to the north and the Florida Keys to the south and east (Figure 2). “Outer Keys” refers to the portion of the Florida Keys that extends beyond Florida Bay’s western boundary (a line drawn approximately between Cape Sable on the southwestern Florida Peninsula and Long Key in the Florida Keys). Florida Bay is composed of a series of shallow basins separated by a network of shoals, mudflats and hundreds of small mangrove islands. GWHs build nests on these islands within Florida Bay and along the outer Keys; they rarely breed on Florida Peninsula or on the main Keys themselves (Robertson 1978). Although some GWHs move to freshwater wetlands on the southern Florida Peninsula during the non-breeding season, the population is essentially non-migratory, and many birds spend the entire year within the Florida Keys ecosystem (Powell and Bjork 1990).

FIGURE 2.

Map of southern Florida including Florida Bay and the Florida Keys. The area in pink indicates the southern end of the general breeding range of Ardea herodias wardi (Florida Fish and Wildlife Conservation Commission 2003, Sibley 2007). The area in green indicates the general breeding range of A. h. occidentalis in south Florida. Red dots mark breeding sites of blue-colored herons, including small dots from this study and large dots from the Florida Bird Breeding Atlas (Florida Fish and Wildlife Conservation Commission 2003); black dots mark breeding sites of white-colored herons examined in this study. Numbered sites indicate the number of genetic samples from each site (B = blue, W = white). Table S2 provides the names and locations of sites with numbers.

Florida Bay’s population breeds asynchronously; nests of white herons can be found at any time of the year. Most breeding activity, however, coincides with south Florida’s dry season, approximately from October through April, with peak nesting of both blue and white individuals of GWH occurring from December to February (Gawlik 1998). Breeding of A. herodias wardi occurs on the Florida Peninsula, just a few km away from the Florida Bay population (Figure 2; Florida Fish and Wildlife Conservation Commission 2003, Sibley 2007). However, its nesting localities are physically distinct (inner mainland vs. islands) and its peak breeding occurs at a different time (from March to June). Thus, to some degree, there are physical and temporal barriers to most GWH-GBH breeding.

Breeding behavior provides many opportunities to observe both members of a mated pair together at their nest. Shared duties and characteristic behaviors between mates make it possible to assign pair status—white-white, blue-blue and white-blue—with confidence (Meyerriecks 1960, Mock 1976, Butler 1992, Vennesland and Butler 2011) and, thus, to determine approximately the degree of assortative mating. However, most mating blue-colored herons in Florida Bay have intermediate plumage (i.e. are Wurdemann’s Herons) and therefore probably vary in the degree of backcrossing, which is impossible to assess in the field. Backcrossing could influence mate choice, such that “bluer” herons favor blue mates and “whiter” herons prefer white mates, causing an underestimate of assortative mating if birds are scored simply as blue and white as we have done (see Methods). Additionally, it is difficult to distinguish male and female birds in the field, which may also influence patterns of mate choice.

Background: Conservation

The challenges imposed on south Florida’s ecosystems by a growing human population provide compelling reasons to study GWH taxonomy, regardless of how we ultimately classify the taxon (full species, subspecies or distinct population segment). GWHs constitute one of the world’s many small populations that persist in remnant ecosystems increasingly influenced by surrounding urban, agricultural and recreational landscapes. The Florida Keys population is largely dependent upon heavily disturbed south Florida ecosystems, including Florida Bay, the Florida Keys and Everglades. Because of its narrow geographic distribution, this population is vulnerable to natural catastrophic events as well as to habitat loss or deterioration resulting from human activities. The more we learn about the evolutionary history of the GWH then the better informed we will be in making management decisions about its future.

METHODS

Morphology

Morphological characters were measured in the following groups: northern GBH (A. h. herodias, hereafter “B-N” for blue north), Florida Peninsula GBH (A. h. wardi, hereafter “B-FP” for blue Florida Peninsula) and the Florida Bay breeding population (hereafter “W-FB” and “B-FB” for white and blue Florida Bay, respectively). No data were available for white herons from the outer Keys.

To minimize the possibility of including winter migrants, B-FP sampling was restricted to A. h. wardi collected in Florida, excluding the Florida Keys, between April 1 and September 30. B-FB and W-FB included blue herons collected during the summer, and GWHs and Wurdemann’s Herons regardless of collection date. B-N included GBHs collected within the A. h. herodias subspecies breeding range regardless of collection date. The characters measured (from museum skins) were: length of exposed culmen, depth of bill at base, length of tarsus, wing chord, length of tail, length of middle toe and length of longest occipital plume. The first 5 are commonly measured in birds; the length of the middle toe is often reported for herons and other large birds, and existing data suggest that the longest occipital plume varies among GWH and GBH populations (Holt 1928). H.L.M. performed all measurements and followed the methods of Proctor and Lynch (1993), using a ruler for the longest occipital plume to the nearest mm, digital calipers for the depth of bill to the nearest 0.1 mm, a ruler with an upright stop for the wing chord to the nearest 0.5 mm, and dividers and a ruler for the remaining variables to the nearest 0.5 mm. Each character per individual was measured at least twice and averaged (Supplemental Material Table S1). Because males are larger than females, we analyzed them separately and omitted birds whose sex was unknown. Because first year birds have shorter plumes than older individuals (e.g., Martínez-Vilalta et al. 2018), we omitted them from the analyses. Only 7 adult B-FP (3 males, 3 females, 1 unknown sex) met the selection criteria. Raw data for these 7 herons are included in Table S1 but were dropped from subsequent statistical analyses because of the small sample size.

For each character, we calculated mean, variance and SD, and tested ANOVA and t-test assumptions of normality and homogeneity of variance. We used the W-test for normality to determine whether observations for each character within each group were normally distributed (Shapiro and Wilk 1965). For females, we used an F-test to test for homogeneity of variance between 2 groups (Sokal and Rohlf 1995). For males, we used Bartlett’s test for homogeneity of variance to test among 3 groups (Sokal and Rohlf 1995). We used appropriate parametric or non-parametric statistics to test for differences between or among means.

Nesting Behavior

To determine whether white and blue herons pair randomly with respect to plumage color, we monitored GWHs and GBHs in Florida Bay during the peak of the 1998–1999 breeding season (from October through to February). We used Leica 8x44 binoculars or a Swarovski 60x spotting scope to view nests from a distance until we observed adults engaged in activities that identified them as a pair (e.g., switching incubation duties) and were able to confirm the phenotypes of both members of the pair. For each blue adult, we attempted to determine whether its plumage was characteristic of a “typical” GBH or a Wurdemann’s Heron. Because of the difficulty in doing this (discussed below), adult phenotype is reported here as either blue or white. We used the number of blue and white adults from the sample of observed nests to estimate the color proportion in the breeding population. We used this estimate to generate expected values for each of the 3 pairing categories (white-white, white-blue and blue-blue). We used a χ² goodness of fit test for the difference between observed and expected values (Sokal and Rohlf 1995), subtracting one degree of freedom for the total sample size and one degree of freedom because sample frequencies were used to generate expected values.

Whenever possible, we recorded nestling phenotypes in nests where adult phenotypes were known. We could not discern any differences in plumage among blue nestlings even upon close inspection (we handled nestlings to take blood and feather samples for genetic analysis). This made it impossible to infer whether the adult phenotype of a nestling would be characteristic of a “pure” blue GBH or Wurdemann’s Heron. Nestling color, therefore, was recorded as either blue or white.

Genetics

To evaluate geographic patterns of genetic variation among GWH and GBH populations we used microsatellite loci. Microsatellites are non-coding sequences composed of repeat units. These units are short, generally 1–6 base pairs (bp) long. Microsatellites evolve through gain or loss of repeat units caused by slippage and misalignment during DNA replication (Schlötterer and Tautz 1992, Schlötterer and Pemberton 1994). Alleles at a locus, therefore, differ in length and are easily resolved by electrophoresis (Queller et al. 1993). We genotyped 13 microsatellite loci in 6 A. herodias groups. These include, in addition to the 4 morphological geographic groups defined above (B-N, B-FP, W-FB and B-FB), A. h. fannini of the Pacific Northwest (hereafter “B-PNW”) and A. h. occidentalis of the outer Keys, all of which are white (hereafter “W-OK”). To obtain DNA, we collected blood and feather samples in the field and requested preserved muscle tissues from museum collections (Figure 2, Supplemental Materials Table S2). We collected B-FP, B-FB, W-FB and W-OK samples from nestlings only, thereby minimizing the possibility of including winter migrants from northern populations in these groups. Only one nestling per nest was included in the genetic analysis (which minimized the chance of including siblings). B-FP samples were collected from nestlings in the southern portion of the A. h. wardi range (Water Conservation Area 3A, Dade County, Florida; the northernmost point in Figure 2), where breeding adults were unlikely to interbreed directly with A. h. herodias (B-N). Approximately 0.1 mL of blood (collected from the tibio-tarsal vein) was mixed with 1.0 mL 10% EDTA anticoagulant/preservative buffer. We allowed red blood cells (RBCs) to settle overnight, removed and discarded the plasma/EDTA supernatant, re-suspended the RBCs in 1.0 mL 10% EDTA, and refrigerated the samples at 4^oC. Feathers were refrigerated at 4^oC during the field season and then stored at –80^oC.

We isolated genomic DNA from tissue samples (muscle, feather or blood) with a DNeasy Tissue Kit (Qiagen Inc., Valencia, California, USA). For muscle samples, we used ~50 mg of tissue, and for feathers ~5 mm from the root of the feather shaft. We identified microsatellite markers (Supplemental Table S3) for A. herodias using the enrichment methods described in McGuire and Noor (2002). We selected a final panel of 13 microsatellite loci that excluded monomorphic loci and sex-linked loci, which can produce a false rejection of the null hypothesis when testing for genetic differentiation (H₀: no difference in allele frequencies among populations) if a bias in sex ratios exists within any of the heron populations sampled. Finally, we excluded individuals that amplified at fewer than half the loci (n = 3).

We checked microsatellite data for evidence of null alleles, allelic dropout and potential scoring problems with MICROCHECKER 2.2.3 (Van Oosterhout et al. 2004). We estimated population genetic variation (expected and observed heterozygosity, allelic richness) in GENETIX 4.05.2 (Belkhir et al. 1999) and FStat 2.9.3.2 (Goudet 2002). We calculated Hardy-Weinberg equilibrium and linkage disequilibrium in GENEPOP 4.1.0, and Factorial Correspondence Analysis (FCA) in GENETIX 4.05.2 (Belkhir et al. 1999).

To estimate genetic structure, we used STRUCTURE 2.3.4 (Pritchard et al. 2000) with an admixture model, 200,000 burn-in/200,000 additional iterations (Gilbert et al. 2012), correlated allele frequencies, and otherwise program default settings. We performed 20 runs for each hypothesis of K = 1–7 (Evanno et al. 2005). We processed STRUCTURE output with STRUCTURE HARVESTER (Earl 2012) and used the CLUMPAK server (Kopelman et al. 2015) to condense data from multiple runs and create a figure of population clusters.

We implemented analysis of molecular variance (AMOVA), isolation-by-distance (IBD) via Mantel tests, and pairwise population differentiation in ARLEQUIN 3.5.2.2 (Excoffier and Lischer 2010) with the allowed level of missing data set to 0.1. STRUCTURE results were used to inform the number of groups created in the AMOVA. AMOVA was performed with 16,000 permutations (Guo and Thompson 1992) in a locus-by-locus analysis (to account for missing data) with results weighted as an average over loci (Excoffier and Lischer 2010). Mantel tests were performed for: (1) all 6 groups (B-PNW, B-N, B-FP, B-FB, W-FB and W-OK); (2) all 4 Florida groups (B-FP, B-FB, W-FB and W-OK); (3) 5 groups in which the B-FB and W-FB birds were combined; and (4) 3 Florida groups in which the B-FB and W-FB birds were combined. All tests were performed using 1,000 permutations and Euclidean geographic distances between population pairs. Pairwise population differentiation was tested using F_ST with the number of Markov chain steps set at 100,000 and the number of dememorizations set at 10,000.

We employed BAYESASS 3.0 (Wilson and Rannala 2003) to examine gene flow among groups of Florida herons, using results from STRUCTURE analyses to inform groupings. We first compared migration rates (m, which is defined as the proportion of first generation migrants) between B-FP, Florida Bay (all herons combined) and W-OK. We adjusted the mixing parameters to achieve acceptance rates in the optimal range (Rannala 2015). Final mixing parameters were 0.30 for the migration rate, 0.60 for allele frequencies and 0.90 for inbreeding coefficients to produce acceptance rates of 0.32, 0.31 and 0.36 for each parameter, respectively. Second, we compared migration rates between B-FP, B-FB, W-FB and W-OK with mixing parameters adjusted to 0.8, 0.7 and 0.95 for migration rate, allele frequencies and inbreeding coefficient, respectively, and acceptance rates of 0.35, 0.31 and 0.43, respectively. In both scenarios we used a 30,000,000 iteration burn-in followed by 30,000,000 additional iterations with a sampling interval of 1,000. We ensured convergence by examining the trace file for log probability in TRACER 1.6 (Rambaut et al. 2014).

For all analyses above, results were considered significant when α ≤ 0.05.

RESULTS

Morphology

We measured 7 morphological characters in 101 GWH and GBH specimens (Supplemental Materials Tables S1, S4 and S5; Figure 3). After excluding first year individuals, specimens of unknown sex and groups with inadequate sample sizes, the final dataset contained 75 herons: 8 B-N (8 male, 0 female), 26 B-FB (14 male, 12 female) and 41 W-FB (24 male, 17 female). No specimens collected from Florida Bay had the “typical” blue GBH plumage; all were Wurdemann’s Herons.

FIGURE 3.

Mean values (±95% CI) for 7 morphological variables measured from museum skins of northern Great Blue Herons (B-N), blue herons from the Florida Keys (B-FK) and white herons from the Florida Keys (W-FK).

Among males, ${\overline{\mathrm{X}}}_{\mathrm{B}-\mathrm{N}}$ differed significantly from both ${\overline{\mathrm{X}}}_{\mathrm{B}-\mathrm{F}\mathrm{B}}$ and ${\overline{\mathrm{X}}}_{\mathrm{W}-\mathrm{F}\mathrm{B}}$ at 5 of the 7 characters: length of exposed culmen, depth of bill at base, length of tarsus, length of middle toe and length of longest occipital plume. Wing chord and tail length did not differ among groups. ${\overline{\mathrm{X}}}_{\mathrm{B}-\mathrm{F}\mathrm{B}}$ and ${\overline{\mathrm{X}}}_{\mathrm{W}-\mathrm{F}\mathrm{B}}$ did not differ significantly at any of the 7 variables in either males or females. Although not significantly different, mean values for these 2 groups followed a consistent pattern. In males and females, ${\overline{\mathrm{X}}}_{\mathrm{W}-\mathrm{F}\mathrm{B}}$ was larger than ${\overline{\mathrm{X}}}_{\mathrm{B}-\mathrm{F}\mathrm{B}}$ at all variables except longest occipital plume (where ${\overline{\mathrm{X}}}_{\mathrm{B}-\mathrm{F}\mathrm{B}}>{\overline{\mathrm{X}}}_{\mathrm{W}-\mathrm{F}\mathrm{B}}$>). ${\overline{\mathrm{X}}}_{\mathrm{B}-\mathrm{F}\mathrm{B}}$ was intermediate between ${\overline{\mathrm{X}}}_{\mathrm{B}-\mathrm{N}}$ and ${\overline{\mathrm{X}}}_{\mathrm{W}-\mathrm{F}\mathrm{B}}$ at 5 of the 7 variables: length of exposed culmen, depth of bill at base, length of tarsus, length of middle toe and length of longest occipital plume.

Nesting Behavior

We determined adult plumage color at 114 nests on 14 islands within Florida Bay during the 1998–1999 breeding season. White and blue individuals were clearly distinguishable. Among blue adults a continuum of phenotypes ranged from those with plumage indistinguishable (under field conditions) from other blue North American GBHs to obvious intermediates. If we define, for a moment, 3 somewhat arbitrary plumage categories (blue herons at one end of the blue plumage continuum, intermediate herons at the other end of the blue plumage continuum, and white herons), all pair combinations were observed and all combinations produced viable offspring. Although some blue adults had plumage that was clearly intermediate and others had plumage that was indistinguishable from “typical” blue GBHs, the continuum of blue phenotypes made it extremely difficult to devise any meaningful criteria to categorize blue adults as either “blue” or “intermediate.” Adult phenotype, therefore, is reported as either “white” or “blue” (Table 1). We rejected the random mating hypothesis using a χ² goodness of fit test for the difference between observed and expected values (χ² = 31.32, df = 1, P < 0.001).

TABLE 1.

Observed (O) and expected (E) values for each pairing category expressed as number of mated pairs and percentage of total.

Adult phenotypes		No. of mated pairs		Percentage of total
		O	Ea	O	Ea
White	White	83	73	73	64
White	Blue	17	36	15	32
Blue	Blue	14	5	12	4

^a Expected values assume p(white) = 0.8, q(blue) = 0.2 and that adults pair randomly with respect to plumage color.

We found only white nestlings in nests where both adults were white (Table 2). Mixed pairs produced broods with all blue offspring, all white offspring and mixed offspring. Blue-blue pairs produced broods that were either all blue or mixed. We did not find any blue-blue pairs with all white offspring; however, the number of nests in this category was small and did not preclude the possibility that 2 blue parents could produce a brood of all white offspring.

TABLE 2.

Number and color of offspring in nests for which adult phenotypes are known.

Adult phenotypes			No. of offspring
		Number of nests	White	Blue
White	White	51	113	0
White	Blue	14	15	13
Blue	Blue	14	6	25

Genetics

We genotyped 210 individuals at 13 microsatellite loci (Supplemental Material Table S6). We found no evidence that any loci were sex-linked: all loci amplified in males and females and heterozygous males and females were found in all individuals. With MICROCHECKER, we found evidence of null alleles (homozygote excess) in the B-N population at locus Ah536 and in the B-FB population at locus Ah630. Otherwise, all other populations and loci showed no evidence of scoring errors, allelic dropout or null alleles. With GENEPOP, we detected 2 deviations from Hardy-Weinberg Equilibrium (HWE) following Bonferroni correction: locus Ah 526 in W-FB and locus Ah 536 in W-OK (both heterozygote deficits). Given that departures from HWE were only found at a single locus in each of 2 populations, allele frequencies in all 6 populations were largely consistent with Hardy-Weinberg expectations.

Genetic variation was similar in all 6 sampling areas, although somewhat lower for B-PNW and W-OK (Table 3). Linkage disequilibrium analysis followed by Bonferroni correction indicated that no loci were linked when each locus pair was tested across all populations.

TABLE 3.

Expected heterozygosity (H_E: calculated without bias; Nei 1978), observed heterozygosity (H_O), average number of alleles per locus (A), allelic richness (AR). SD is given in parentheses for H_E and H_O.

Color-locationa	HE	HO	A	AR
B-PNW	0.572 (0.213)	0.553 (0.248)	3.769	3.736
B-N	0.622 (0.232)	0.587 (0.232)	5.692	4.322
B-FP	0.667 (0.161)	0.628 (0.208)	5.462	4.652
B-FB	0.597 (0.240)	0.565 (0.252)	5.154	4.076
W-FB	0.576 (0.234)	0.539 (0.232)	6.000	3.995
W-OK	0.522 (0.255)	0.498 (0.247)	5.231	3.749

^a B = blue plumage, W = white plumage, PNW = Pacific Northwest, N = north-central United States, FP = Florida Peninsula, FB = Florida Bay, OK = outer Keys.

The FCA plot suggested loose population structuring: there were many individuals from a sampling location clustered together, but there were usually some individuals that clustered with individuals from a different sampling location (Figure 4).

FIGURE 4.

FCA for 6 heron sampling locations and 13 loci. Yellow = B-PNW (blue plumage-Pacific northwest); blue = B-N (blue plumage-north); white = B-FP (blue plumage-Florida Peninsula); gray = B-FB (blue plumage-Florida Bay); pink = W-FB (white plumage-Florida Bay); green = W-OK (white plumage-outer Keys).

In STRUCTURE, inspection of α, F, D and the likelihood plots indicated a sufficient burn-in and number of post-burn-in iterations. STRUCTURE HARVESTER indicated the best support for K = 2 (Figure 5) based on the method of Evanno et al. (2005). Cluster 1 (Figure 6) consisted of GBH from the Pacific Northwest (B-PNW), north-central region (B-N) and the Florida Peninsula (B-FP). Cluster 2 included blue and white individuals from Florida Bay (B-FB and W-FB) and GWH from the outer Florida Keys (W-OK). Mean Q-values (Table 4) indicated that individuals in cluster 1 had >80% ancestry from cluster 1, whereas W-FB and W-OK had >74% ancestry from cluster 2, and B-FB had ~42% ancestry from cluster 1 and ~58% ancestry from cluster 2. Some individuals in every sampling location had less ancestry in common with their own cluster than with the other cluster (Figure 7; Supplemental Material Table S7).

TABLE 4.

STRUCTUE mean ancestry (Q) values in each of 6 sampling locations for each cluster.

Color-locationa	Cluster 1	Cluster 2	n
B-PNW	0.809	0.191	11
B-N	0.873	0.127	30
B-FP	0.849	0.151	20
B-FB	0.418	0.582	35
W-FB	0.257	0.743	77
W-OK	0.144	0.856	37

^a B = blue plumage, W = white plumage, PNW = Pacific Northwest, N = north-central United States, FP = Florida Peninsula, FB = Florida Bay, OK = outer Keys.

FIGURE 5.

Plot of ΔK in STRUCTURE indicating K = 2 as the most likely number of clusters.

FIGURE 6.

STRUCTURE bar plot showing the proportion of ancestral alleles in each individual in 6 sampling locations. Purple = blue ancestry; gold = white ancestry. Two clusters were identified: the first consisted of B-PNW (blue plumage-Pacific northwest), B-N (blue-plumage north) and B-FP (blue plumage-Florida Peninsula), whereas the second consisted of B-FB (blue plumage-Florida Bay), W-FB (white plumage-Florida Bay) and W-OK (white plumage-outer Keys).

FIGURE 7.

Proportion of immigrants and SD of the marginal posterior distribution as estimated in BayesAss. Blue font: estimates of non-immigrants in each sampling location; black font: estimates of immigration between sampling locations. W-OK (white plumage-outer Keys), W-FB (white plumage-Florida Bay), B-FB (blue plumage-Florida Bay), B-FP (blue plumage-Florida Peninsula). (Top panel) Three-group scenario in which all herons in Florida Bay (FB) were combined into a single group. (Bottom panel) Four-group scenario in which blue and white plumage birds in Florida Bay were considered separately.

We analyzed data in ARLEQUIN by assigning heron sampling locations to 1 of 2 groups corresponding to the 2 clusters identified in STRUCTURE (group 1 = B-PNW, B-N, B-FP; group 2 = B-FB, W-FB, W-OK). Pairwise F_ST between group 1 and 2 was low but significant (0.069, P < 0.0001); similarly, pairwise F_ST among the 6 sampling locations was generally low (Table 5) but significant (P ≤ 0.027). The lowest F_ST value was observed between W-FB and B-FB (0.007), whereas the highest F_ST values were observed between B-PNW and W-OK (0.187) and between B-N and W-OK (0.131), the 2 most geographically distant pairs of sampling locations. AMOVA suggested that most variation was found within individuals (86.4%) with little variation explained by the 2 groups (5.72%) or among populations (i.e. sampling location) within groups (2.76%; Table 6). Three Mantel tests showed significant correlation between genetic (F_ST) and geographic (km) distance: (1) all 6 groups (r = 0.71, P = 0.003); (2) all 4 Florida groups (r = 0.95, P = 0.045); and (3) 5 groups in which the B-FB and W-FB birds were combined (r = 0.59, P = 0.007). The Mantel test that included the 3 Florida groups in which the B-FB and W-FB birds were combined was not significant (r = 0.99, P = 0.147).

TABLE 5.

Pairwise F_ST between each sampling location (all significant at P ≤ 0.027). Pairwise F_ST between group 1 and group 2 = 0.069 (P < 0.0001).

Color-locationa	B-N	B-FP	B-PNW	B-FB	W-FB	W-OK
B-N	0.000
B-FP	0.014	0.000
B-PNW	0.083	0.106	0.000
B-FB	0.066	0.041	0.114	0.000
W-FB	0.081	0.057	0.114	0.007	0.000
W-OK	0.131	0.107	0.187	0.025	0.025	0.000

^a B = blue plumage, W = white plumage, PNW = Pacific Northwest, N = north-central United states, FP = Florida Peninsula, FB = Florida Bay, OK = outer Keys.

TABLE 6.

Global AMOVA results as a weighted average over 13 loci. Two groups were compared based on STRUCTURE analysis. Group 1: B-PNW (blue plumage-Pacific northwest), B-N (blue plumage-north) and B-FP (blue plumage-Florida Peninsula). Group 2: B-FB (blue plumage-Florida Bay), W-FB (white plumage-Florida Bay) and W-OK (white plumage-outer Keys).

Source of variation	Sum of squares	Variance components	Percentage variation
Among groups	51.361	0.238	5.716
Among populations within groups	44.865	0.115	2.757
Among individuals within populations	791.553	0.214	5.156
Within individuals	729.000	3.590	86.371
Total	1616.779	4.157

STRUCTURE results indicated considerable admixture in the Florida Bay region particularly in blue herons; therefore, we examined gene flow (as estimated by BayesAss) in 2 ways (Figure 7). The first compared B-FP, Florida Bay (all herons combined) and W-OK. These results indicated that most individuals were derived from their sampling location (Figure 7). However, there were some interesting patterns. Most notably, ~29% of blue herons on the Florida Peninsula derive from the Florida Bay area, and ~32% of blue and white herons in Florida Bay derive from white herons in the outer Keys, indicating eastward gene flow (Figure 7). In contrast, only ~1–2% of white herons in the outer Keys derive from either B-FP or Florida Bay (Figure 7). Inbreeding coefficients were 0.40 (SD = 0.28) for B-FP, 0.15 (SD = 0.10) for Florida Bay herons and 0.09 (SD = 0.02) for W-OK, indicating that herons in the outer Keys are not particularly inbred despite a limited number of immigrants. Note that B-FP’s small sample, n = 20, and acquisition from a single area (although there were different nesting sites within that area), may account for its high inbreeding coefficient; it may also explain why this population is always involved in high F_ST values when compared with other populations. Our second analysis compared gene flow between sampling locations, but took into account plumage color in Florida Bay, resulting in 4 groups: B-FP, B-FB, W-FB and W-OK. Again, most birds in each sampling location consisted of non-immigrants, but the number of non-immigrants in W-OK was 96%, whereas in B-FP, B-FB and W-FB it was 68% (all 3 sampling locations; Figure 7). The largest proportion of immigrants in B-FP came from B-FB (0.25); in B-FB from W-OK (0.28); in W-FB from W-OK (0.28); and in W-OK from B-FB (0.02) and W-FB (0.01). These results suggest that blue herons in Florida Bay tend to disperse to the Florida Peninsula and white herons in the outer Keys tend to disperse to Florida Bay (Figure 7). Inbreeding coefficients were 0.43 (SD = 0.29) in B-FP, 0.04 (SD = 0.04) in B-FB, 0.38 (SD = 0.17) in W-FB and 0.1 (SD = 0.02) in W-OK.

DISCUSSION

Morphology

The pattern of plumage variation between GBH and GWH at the southern tip of Florida is consistent with 2 hypotheses: the intermediate plumage of Wurdemann’s Herons results from hybridization between B-FP and W-FB or it results from a cline. Until now, only one morphometric character (occipital plume length) had been compared between B-FP and W-FB/OK (Holt 1928). We compared 6 additional characters and found trends in size, but not a significant difference. Unfortunately, our analysis suffered from a lack of W-OK individuals and a dearth of B-FP individuals, and we measured only a small number of characters, several of which (culmen, depth of bill, tarsus and middle toe) are almost certainly correlated. Because of the small sample size, the power of our morphometric tests was low (1-β < 0.30 for all tests of H₀: μ_B-FB = μ_W-FB).

The size trends were signaled by mean values of 5 characters, with ${\overline{\mathrm{X}}}_{\mathrm{B}-\mathrm{F}\mathrm{B}}$ intermediate between ${\overline{\mathrm{X}}}_{\mathrm{B}-\mathrm{N}}$ and ${\overline{\mathrm{X}}}_{\mathrm{W}-\mathrm{F}\mathrm{B}}$ regardless of whether the character was larger in B-N (longest occipital plume) or W-FB (culmen, depth of bill, tarsus and middle toe). The dearth of specimens of W-OK and B-FP is particularly unfortunate because these 2 areas represent the extremes of plumage coloration in Florida (indeed, in the entire species), whereas Florida Bay—the source of our best sampling—represents mostly intermediate individuals. Zachow’s (1983) skeletal data suggest that B-FP are smaller than GWHs, but these differences are probably small compared with the substantial size differences we found between B-N and B-FP/W-FB. Overall, even if our sampling were better, morphometric comparisons would not resolve the hybridization-cline issue.

Nesting Behavior

Support for GWH as a subspecies of GBH stems largely from limited, usually anecdotal, reports of interbreeding between blue and white individuals across the Caribbean and the impression that mate choice is random with respect to plumage color (Holt 1928, Mayr 1956, Meyerriecks 1957, Bond 1961). However, no published data support a random mating hypothesis. Ours is the first study to conduct observations of nests in sufficient numbers and detail to confirm the phenotypes of both members of a mated pair and their offspring. These observations indicate more white-white and blue-blue pairs and fewer mixed pairs than expected in a randomly mating population. The findings are consistent with Robertson’s (1978) thesis of positive assortative mating. However, Robertson (1978) reported a 10-fold reduction in mixed pairing, whereas we observed half the mixed pairs expected in a randomly mating population. The difference suggests that either the proportion of blue herons breeding in Florida Bay fluctuates over time, or one (or both) of our estimates is (are) incorrect. Robertson’s (1978) surveys were done at 2-to-4-mo intervals from June 1959 through to May 1960. Aerial surveys (Gawlik 1998) show that peak nesting activity is essentially the same for white and blue herons in Florida Bay (approximately December to February), but that white heron breeding is more protracted, with some individuals breeding in virtually every month. Our estimate of relative proportions of white to blue herons in the Florida Bay breeding population (80% white, 20% blue) is consistent with surveys conducted at the time of our research (e.g., Gawlik 1998). If our data and Robertson’s are both accurate, then the discrepancy invites the following hypothesis. As blue herons become rare within the Florida Bay breeding population, the probability that a blue individual will pair with a white individual increases but remains below the level expected in a randomly mating population. This hypothesis is consistent with mate choice theory that predicts the degree of “choosiness” exhibited by courting animals will be influenced by the availability of potential mates (Crowley et al. 1991, Nuechterlein and Buitron 1998). Choosiness may also be influenced by the degree of backcrossing in the blue-colored Wurdemann’s Herons.

Although the pattern of mate choice in Florida Bay is assortative, the mechanism producing this pattern is unknown. The non-random pattern does not necessarily imply that the herons use plumage color as a criterion for mate choice. Other factors (habitat preference, timing of breeding, sex ratios, historical geographic distribution, etc.) may also influence patterns of mate choice. These factors may function at several spatial or temporal scales and could either inhibit or promote mixed pairing. For example, subtle differences in habitat preference (nest sites or foraging habitat) could segregate white and blue birds during the breeding season. Such differences might manifest themselves on small scales (within individual islands) or on larger scales (among islands). If white and blue herons occupy different habitats and mate choice is restricted primarily to birds within a preferred habitat, then the probability of encountering potential mates of the opposite phenotype is reduced and there will be fewer mixed pairs than expected in a randomly mating population.

Conversely, other factors may facilitate mixed pairing even if the herons prefer mates with similar phenotype. In Florida Bay, pairs nest singly or in loose association with other breeding birds, but rarely in the dense colonies typical of GBH populations on the mainland. Even large islands generally have <30 nesting pairs at any given time and most have <10 (Gawlik and Ogden 1996). Breeding is asynchronous. Not all birds acquire the visible signs of breeding condition (brilliant soft part coloration and elongated neck, back and occipital plumes) simultaneously, and nests at different stages of the breeding cycle (egg, nestling, fledgling) are commonly found in proximity (H. L. McGuire personal observation). Blue individuals are relatively rare in the Florida Bay breeding population. If mate choice occurs at a small spatial scale (within individual islands or among closely spaced islands), some blue herons may not encounter a suitable mate of the preferred phenotype. These birds may pair with the less desired phenotype rather than give up an opportunity to mate during a breeding season. Thus, some mixed pairs could occur even if white and blue herons have strong preferences for mates with similar phenotype.

Although positive assortative mating was demonstrated (Table 1), the pattern of mate choice with respect to plumage color suggests that prezygotic reproductive barriers are imperfect. For example, 17 of 114 nests comprised mixed pairs (15%), although this percentage would be an overestimate if the degree of backcrossing is a factor in mate choice. Without a more detailed understanding of some of the variables that influence mate choice, we can only say at this stage that blue and white herons prefer to mate with their own kind.

Blue–White Polymorphism

Our nesting data provide clues into the genetic basis of the plumage polymorphism observed in Florida Bay. Some ornithologists have suggested that 2 white adults might be capable of producing blue offspring (Mayr 1956, Meyerriecks 1957). Mayr (1956) proposed a model in which a dominant allele conferred white plumage and modifier genes were responsible for producing the intermediate plumage of Wurdemann’s Heron. However, we found only white offspring in nests where both parents were white. If plumage color is determined primarily at a single locus and white is dominant, then the probability of our finding 113 white and zero blue offspring in nests where both parents were white (Table 2) requires that most pairs (>88%) have at least one member homozygous for the dominant white allele. This is possible, given assortative mating and our estimate of the ratio of white to blue herons in the breeding population (4:1). However, when both parents were blue, we found both blue and white offspring in nests. If plumage color is controlled primarily at a single locus, 2 blue parents can produce white offspring only if white plumage is a recessive trait and both parents are heterozygotes. Furthermore, the proportion of white offspring found in nests of blue-blue pairs (19.4 ± 13.9%) is within the range expected under the hypothesis that white plumage is recessive (0–25%) and differs significantly from the expected value for the hypothesis that white plumage is dominant (0%).

Although nestling and adult phenotype data indicate that white plumage behaves as a recessive trait, the single locus hypothesis for inheritance of plumage color ignores the range of blue phenotypes found in the Florida Bay population. Several hypotheses can explain these intermediate plumages (e.g., incomplete dominance at a single locus or additive alleles at more than one locus). Testing any of these hypotheses would require genomic data or the examination of large numbers of offspring from known crosses and the ability to determine their adult phenotypes. Both are beyond the reach of our data. Moreover, we are assuming no extra-pair paternity, which is a tenuous assumption for colonially breeding herons (Ramo 1993, Krebs et al. 2004).

Regardless of whether the allele conferring white plumage is dominant or recessive, it appears to be expressed only in the Florida Keys and Caribbean populations. White individuals occasionally occur in other North American GBH populations, but these are generally considered vagrant, or possibly leucistic, individuals rather than locally breeding birds (Mitra and Fritz 2002). There is only one published observation of a white nestling outside of south Florida: in Texas (McHenry and Dyes 1983). The lack of likely breeding white individuals, and overwhelming numbers of blue individuals, in mainland GBH populations suggests that most North American GBHs do not carry an allele conferring white plumage. This implies that there is little emigration of either white or blue individuals from the predominantly white Florida Bay population to other GBH populations, or introgression of the white gene into mainland populations (but see the discussion of genetics below).

Inferring immigration from other GBH populations into the Florida Bay population based on mating preferences is more difficult. Although fewer mixed pairs occurred than expected from random mating, the number that we observed (15%) is not trivial. These matings could provide gene flow between white and blue herons in Florida Bay or between the Florida Bay breeding population and other GBH populations.

Genetics

STRUCTURE analysis of herons sampled from 6 A. herodias taxa revealed 2 groups (Figure 6), one consisting of A. h. fannini (B-PNW), A. h. herodias (B-N) and A. h. wardi (B-FP), and the other consisting of A. h. occidentalis (B-FB, W-FB and W-OK). However, IBD analysis across the entire GBH range suggests that most differentiation is attributable to geographic separation (for IBD within Florida, see below). Gene flow may have had a significant influence on genetic structure because all sampled taxa had some individuals with ancestry (Q) values more consistent with immigrant status, but incomplete lineage sorting may also be responsible for allelic similarities, especially between the most distant populations, which are least likely to experience gene flow (e.g., W-OK and B-PNW). FCA analyses showed that individuals from each of the 2 STRUCTURE groups largely clustered together, but some individuals in each group clustered most closely with individuals from a different sampling location (Figure 4). F_ST values, although significant in all pairwise comparisons, were fairly low except between W-OK and B-PNW and between W-OK and B-N (F_ST = 0.187 and 0.131, respectively), which is unsurprising given the distances that separate W-OK from B-PNW and B-N. All pairwise comparisons involving B-FP were at the high end of the F_ST range, but these may have been influenced by B-FP’s small sample size and single genetic sampling locality (Figure 2). Finally, ARLEQUIN results indicated that little genetic variation among the herons is attributable either to the 2 STRUCTURE groups or to the 6 geographic groups; most genetic variation is within individuals (Table 6).

Within Florida, IBD does not appear to play a large role in determining genetic structure; when all Florida groups were examined (B-FP, B-FB, W-FB and W-OK) the result was marginally significant (r = 0.95, P = 0.045), and when B-FB and W-FB birds were combined the result was not significant (r = 0.99, P = 0.147), although the correlation coefficients were large in both analyses. STRUCTURE analyses suggested that herons in Florida Bay, particularly birds with blue plumage, have a high level of admixture from each of the 2 clusters (Table 4). This finding suggests they are hybrids rather than intermediates responding to a clinal selection gradient (assuming the microsatellite data are selectively neutral). BAYESASS results indicated considerable gene flow from Florida Bay blue herons to the Florida Peninsula, and that the W-OK group receives few immigrants even though herons in the outer Keys do appear to disperse to Florida Bay (Figure 7). Given that BAYESASS may overestimate migration rates (Samarasin et al. 2017), gene flow to and from the outer Keys may be even lower than our results indicate. A distinct population in the outer Keys would support Meyerriecks (1957) proposal that the gap between Florida Bay and the outer Keys might split south Florida’s GWHs into 2 distinct breeding populations, although Robertson (1978) doubted the existence of this gap, and our results do not indicate such a difference.

Although our microsatellite data provide useful insight into the genetic distinctiveness between GWH and GBH, more precise information should be forthcoming from next generation sequencing of these taxa. Genomic analyses based on large sets of anonymous SNP loci (e.g., RAD-seq or ultra-conserved elements) or loci associated with plumage coloration and/or adaptation to distinct habitats (e.g., marine vs. freshwater) may reveal greater genetic differences than observed here with non-coding microsatellite loci, which are predominantly influenced by genetic drift. For example, Poelstra et al. (2014) showed that genetic differentiation between Carrion (Corvus corone) and Gray-coated Hooded Crows (C. cornix) only occurred at areas of the genome related to visual perception and plumage pigmentation, but because the 2 crow species mate assortatively with respect to plumage coloration, these small areas of genetic differentiation may drive reproductive isolation and therefore speciation. If the assortative mating of GWH is based on plumage coloration, as seems likely, population differentiation may be occurring in small islands of the genome, and although such loci might ultimately drive speciation, they would not be detected in our analyses.

Are GWH and GBH Distinct Species?

We can say with some assurance that the morphological trend from GWH in the outer Keys to GWH and Wurdemann’s Herons in Florida Bay to GBH on Florida Peninsula, although crossing a saltwater-freshwater ecotone, is probably not a cline driven by habitat selection on morphology. STRUCTURE analyses indicate that blue plumage birds in Florida Bay are highly admixed with genes coming from both parental populations (Table 4). Wurdemann’s Herons are hybrids and, thus, GWH and GBH are not simply a single evolutionary group separated by habitat and on their way to “ecological speciation” (Schluter 2009); they are most likely 2 distinct populations that have come together in a secondary contact zone in Florida Bay, possibly as a result of Pleistocene sea level changes influencing the shape and position of the Florida Peninsula (Mayr 1956, Lazell 1989).

With respect to assortative mating in the contact zone, blue individuals choose blue and white individuals choose white most of the time. Although 15% of the nesting pairs are mixed, this number may well be an overestimate given likely backcrossing and a range of “blue-ness” in the birds we observed. As for gene flow, BAYESASS analysis indicates that it is largely in one direction, from the outer Keys to Florida Bay, and via blue-colored herons from Florida Bay to the peninsula. There is no evidence yet that B-FP herons are moving into Florida Bay and breeding there to a substantial degree. Indeed, breeding data in addition to genetic data suggest the opposite: with respect to GWH and Wurdemann’s Heron, B-FP have a disjunct breeding distribution (Figure 2) and different peak breeding season (GWH: October–April, GBH: February–March).

Are GWH and GBH different species? Despite the importance of the species category as a way to organize our understanding of the biological world, some lines are difficult to draw—as evidenced by the proliferation of species concepts, longstanding debates among “lumpers” and “splitters”, and discussions among the authors of this paper. Highly differentiated populations that never interbreed or undifferentiated populations that interbreed freely rarely present problems for those attempting to delimit species boundaries. However, there are many cases where the level of actual or potential interbreeding is ambiguous, and the BSC offers no concrete rules to govern what level renders 2 populations conspecific. Many find this to be a weakness, but it is precisely those cases that challenge the BSC in which we may find the most interesting biological processes in progress. Given available data and a certain degree of subjectivity in species recognition, we believe the weight of evidence favors separating the GWH as a distinct species. Applying the logic of de Queiroz (1998, 2007), GWH is a metapopulation with its own evolutionary trajectory; it is separated by diagnostic plumage (white vs. blue), behavior (mate choice and timing of mating) and habitat (saltwater-island vs. fresh/brackish water mainland). The only remaining issue is distinctiveness in the face of a small amount of gene flow, which is not a taxonomically insurmountable problem, as we commonly recognize distinct species that occasionally exchange genes (Zarza et al. 2016).

Conservation Considerations

GWH is a small, distinct population in an endangered environment and, as such, more likely to go extinct than large populations spread over a wide area (Goodman 1987). Of particular concern is the steady decline of GWHs since the 1990s—even as other wading bird populations in the Everglades and Florida Bay have increased (Meyer and Kent 2011). GWH decline is linked to habitat loss and deterioration resulting from human activities (Hunter et al. 2006, Meyer and Kent 2011). Given the continued threats to south Florida’s ecosystems from human population growth, expanding agriculture, introduction of exotic species, and likely increase in the occurrence of natural disasters, the Keys’ GWH population will almost certainly face serious challenges to its survival in the future.

In comparison to other North American birds, the distribution and breeding behavior of GWH outside of Florida is remarkably poorly documented. This lapse probably results from birdwatchers and ornithologists in the greater Caribbean conflating (listing, recording, treating) breeding GWHs and wintering GBH as the same species—GBH—rather than 2 distinct entities. We predict that our understanding of GWH distribution and breeding will improve dramatically when it is treated as a full species by taxonomists.

ACKNOWLEDGMENTS

The Burke Museum, Bell Museum of Natural History, Field Museum (Chicago), and the Louisiana State University Museum of Natural Science (LSUMNS) provided heron tissues for genetic analysis. We thank the Academy of Natural Sciences of Drexel University, Carnegie Museum of Natural History, Field Museum, Cleveland Museum of Natural History, Museum of Comparative Zoology, LSUMNS and the National Museum of Natural History for access to their bird specimen collections. We are grateful to James Dame, Peter Frederick, Becky Hylton, Winslow McGuire, Melissa Powell, Martin Ruane, Kristina Skarin and Erin Waggoner for their help in the field. We thank the members of H.L.M.’s dissertation committee—Frank Rohwer, F.H.S., William Kelso, Cheryl Hedlund, J. V. Remsen and Mike Stine—for their guidance and support. Dr. Mohamed A. F. Noor and the people in his lab were extraordinarily generous with their time, resources, and expertise isolating and analyzing genetic markers. We are also grateful to J. V. Remsen, Dan Lane and the LSUMNS Vert Lunch group, and Kristen Ruegg, Kevin Winker and an anonymous reviewer for comments that improved the manuscript. Adam Vaccarella kindly assisted in preparing some of the figures and tables. Finally, we would like to extend special thanks to Laura Quinn and the Florida Keys Wild Bird Center for providing tissue samples and logistical support in the field.

Funding statement: This research project was funded by an Allan D. and Helen G. Cruickshank Award from the Florida Ornithological Society, a Research Award from the American Ornithologists’ Union, Florida Fish and Wildlife Conservation Commission grant NG96-025, and an LSU Graduate School Fellowship and Dissertation Fellowship to H.L.M. This material is based upon work that is supported by the National Institute of Food and Agriculture, U.S. Department of Agriculture, McIntire Stennis project LAB04066 and LAB94169.

Ethics statement: The research was conducted under permits from the Florida Fish and Wildlife Conservation Commission (WX96214), U.S. Fish and Wildlife Service (823445), Everglades National Park (1997-0099, 1998-0072 and 1999-0133), Key West and Great White Heron National Wildlife Refuges (96-027, 98-016, and 98-026), and was approved by the Institutional Animal Care and Use Committee at LSU (A96-44).

Author contributions: H.L.M. conceived and designed the project, conducted the fieldwork, produced and analyzed the data, and wrote the project as her Ph.D. dissertation (McGuire 2001). S.S.T. analyzed the genetic data using methods unavailable in 2001, and she edited the manuscript. F.H.S. converted the dissertation to manuscript and edited various drafts.