On 24 March 1989, the T/V Exxon Valdez ran aground on Bligh Reef in Prince William Sound, Alaska, spilling approximately 41 million liters of crude oil (Piatt and Lensink 1989, Piatt et al. 1990). Subsequent wind and ocean currents spread the oil southwest through Prince William Sound (PWS), along the Kenai and Alaska Peninsulas, and along the Kodiak Archipelago. Much of the oil was deposited in nearshore intertidal and subtidal habitats (ADEC 1992, Neff et al. 1995), which are of importance to a large number of vertebrates including molting and wintering waterfowl. Indeed, the nearshore environment of Prince William Sound received about 40% of the oil spilled (Galt et al. 1991). The effect of this oil spill to resident fish and wildlife has been dramatic (Piatt et al. 1990, ECI 1991) and the subject of extensive investigations. Harlequin Ducks (Histrionicus histrionicus) are year-round inhabitants of nearshore environments within the oil spill zone (Isleib and Kessel 1973, Agler et al. 1994). Two hundred and twelve Harlequin Ducks were recovered from beaches after the oil spill (J. Piatt pers. comm.). The actual number of Harlequin Ducks that died from oiling is suspected to be much higher because the recovery of waterbirds was probably very low due to birds being scavenged, sinking or not being found along shorelines (Ford et al. 1987, ECI 1991). Adjusted estimates which correct for these recovery problems range from 1,298 to 2,650 Harlequin Ducks killed from oiling (ECI 1991, Piatt and Ford 1996, J. Piatt pers. comm.). Postspill studies suggest continuing constraints to recovery from the spill, based on differences in winter survival between oiled and unoiled areas (D. Esler unpubl. data), declines in numbers of molting birds within the spill zone (D. Rosenberg unpubl. data), and detectable levels of hydrocarbons in Harlequin Ducks and their prey from 1989 to 1993 (Patten 1995). To understand the process of post-spill recovery of Harlequin Ducks, and to identify impediments to recovery, it is critical to determine whether aggregations of individuals within local areas of the marine environment are discrete and demographically independent, i.e. whether the population is structured within the oil spill region. For example, if birds located in oiled and unoiled areas belong to demographically distinct population segments (i.e. management units; Moritz 1994), then certain segments of the population may have been (and continue to be) impacted disproportionately by oil. Alternatively, lack of spatial structure would imply that oil effects are distributed throughout a larger, panmictic population. Unfortunately, for Harlequin Ducks, very little is known about the level of movements (or gene flow) that occur among areas of the marine environment. In the absence of direct observational data, population discreteness can also be inferred using genetic markers (Slatkin 1985, 1995, Slatkin and Barton 1989). In this study we used three classes of molecular genetic markers which differ in their mode of inheritance (bi-parental, maternal and sex-linked) to determine whether there was population structuring among winter aggregations in coastal marine habitats of PWS and APKA. Information on population structure was then used to assess the potential constraints of the Exxon Valdez oil spill on the recovery of Harlequin Ducks.
Methods. –We compared genetic characteristics of four aggregations of molting birds (referred to as populations) captured within each of two regions (PWS and APKA, Fig. 1). Because band resightings and radio telemetry data indicate Harlequin Ducks almost always winter at or near their molting sites (Robertson 1997, D. Esler unpubl. data), we considered these molting aggregates to represent discrete winter populations. Sampling sites within PWS region included Bay of Isles, the Green Island area, the Foul Bay area and Montague Island, whereas sampling sites within APKA region included Uganik and Uyak bays and Afognak Island of the Kodiak Archipelago, and southern portions of the Alaska Peninsula. Molting Harlequin Ducks were captured at each of the eight sampling sites during July to Sept of 1995, 1996 and 1997 by herding the flightless birds into pens (after Clarkson and Goudie 1994). The sex of each duck was determined by plumage characteristics and aged by bursal probing (Mather and Esler in press). Blood was extracted from the jugular or tarsal vein and preserved with lysis buffer (Longmire et al. 1988). DNA was extracted using Puregene DNA extraction kits (Gentra Systems, Inc., Minneapolis, MN) and standard proteinase-K, phenol-chlorophorm methods (Sambrook et al. 1991). Thirty-three microsatellite loci (Fields and Scribner 1997, Cathy et al. 1998, Buchholtz et al. 1998; Scribner, unpubl. data) were examined for variation by screening two individuals each from four populations on the west coast of North America. Of these, four bi-parentally inherited loci (Sfi?4, Hhi?2, Hhi?5, and Bca?10) and two sex- linked loci (Sfi?1, Bca?4) were variable. The sex-linked (Z-specific) loci provide an estimate of male-mediated gene flow from the preceding generation if sampling is conducting using females. Microsatellite loci were assayed using the polymerase chain reaction, employing end-labeled (32P-?ATP) primers. Specific conditions for each locus are available from R. Lanctot. Products were visualized on 6% denaturing polyacrylamide sequencing gel after autoradiography. A M13 sequencing reaction (Amersham Life Sciences Sequenase Kit, Arlington Heights, IL) and individual standards of known genotypes were run adjacent to the samples to provide an absolute-size marker for determining the size of microsatellite alleles. Mitochondrial DNA (mtDNA) specific primers were designed so that mtDNA sequences, and not nuclear DNA sequences originating from transposed mtDNA (numt, Sorenson and Fleischer 1996), would be amplified. This was verified by sequencing mitochondrial rich heart and mitochondrial poor blood DNA extracts from the same individual as described in Sorenson and Quinn (1998). Using the mtDNA specific primers, we amplified a ?385 base-pair fragment of the 5' end of the control region which is comparable to the hypervariable region 1 of the mammalian mtDNA control region (Vigilant 1990, Wakely 1993). Mitochondrial DNA specific primers included HADUM1L (L16744, ref. Chicken, Desjardins and Morais 1990): 5' TGC CCG AGA CCT ACG GCT C 3'; HADUM2L (L12, ref. Chicken, Desjardins and Morais 1990): 5' TCT AAA ATG ACT CAA CAG TGC C 3'; and HADUMITH (H737, ref. Chicken, Desjardins and Morais 1990): 5' TGA GTA ATG GTG TAG ATA TCG 3'. Nuclear specific primers included HADUN1L (L16744, ref. Chicken, Desjardins and Morais 1990): 5' TAC CCG AGA CCT ACA GCT T 3' and HADUNUCH (H737, ref. Chicken, Desjardins and Morais 1990): 5' TGA GTT ATG GTG TAG ATA CTA 3'. MtDNA sequences (Genbank accession numbers AF101372-81) were obtained by matching either the HADUM1L and HADUM2L primers with the HADUMITH primer. Mitochondrial DNA was PCR-amplified, purified and sequenced using Sequitherm's EXCEL DNA Sequencing Kits (Epicentre Technologies, Madison, WI) and 1.5 mM of HADUMITH primer (contact S. Talbot for specific information). Sequences were visualized using autoradiography, manually scored and aligned. Genetic analyses were restricted to birds ?1 year of age. Data from males and females were used for the bi-parentally inherited loci, but were restricted to females for the maternally inherited mtDNA and sex-linked microsatellite loci. For each population, all bi-parentally inherited loci were tested for linkage (all two-locus comparisons) and Hardy-Weinberg disequilibrium using the Fisher's Exact Test in the Genetics Data Analysis (GDA) program (Lewis and Zaykin 1998). P-values were adjusted for the number of statistical tests. Mean number of alleles per locus were calculated using BIOSYS (Swofford and Selander 1989). We estimated observed (Ho) and expected (He) heterozygosity under Hardy-Weinberg assumptions for each locus (BIOSYS) and for each population and locus (GDA). These estimates of Ho and He were used to generate inbreeding coefficients (F = 1-[Ho/He]) combined across loci for each population (Wright 1951) and tested for significance as described in Li and Horivitz (1953). Statistical analyses of spatial heterogeneity in gene frequency between and within regions for each locus were assessed using hierarchical F-statistics (Weir 1996) in the GDA program at three levels: (1) among individuals within populations, (2) among populations within regions, and (3) between regions. We also calculated RST, an analogue of FST (Michalakis and Excoffier 1996) using the Analysis of Molecular Variance (AMOVA) program (Excoffier et al. 1992). Significance of FST values was based on 95% confidence intervals determined by bootstrapping across loci. Confidence intervals that overlap zero were considered non-significant. We used allele frequencies to calculate Cavalli-Sforza and Edwards (1967) chord distances among the populations and constructed boot-strapped population trees using subroutines within the PHYLIP program (Version 3.572c, Felsenstein 1993). This distance metric has been described as producing robust tree topologies for groups separated evolutionarily over time periods comparable to these populations (Takezaki and Nei 1996). We used birds sampled at Shemya, on the outer Aleutian Islands of Alaska, as an outgroup in this analysis. For the Sfi?1 and Bca?4 sex-linked loci, we calculated allele frequencies and a measure of genetic diversity (D, equation 8.3, Nei 1987). Estimates of variance among individuals within populations (?SC), among populations within regions (?CT), and between regions (?ST) were derived using the AMOVA program. Mitochondrial DNA sequence haplotypes were assigned based on at least a single base-pair substitution or insertion/deletion across the 163 base-pair segment. Pairwise haplotype distances were calculated by the MEGA program (Kumar et al. 1991) and used to estimate haplotype (h) and nucleotide diversity (?) indices for nonselfing populations (equations 8.1 and 10.4, respectively of Nei 1987) using the REAP program (McElroy et al. 1991). We tested heterogeneity of genotype distribution among samples using Monte- Carlo resampling and the chi-square test of Roff and Bentzen (1989). This approach is suitable for genetic data matrices in which many or most elements are very small (<5) or zero. Estimates of regional, population and individual variance (F-statistics) in haplotype frequency were assessed using the AMOVA program as discussed by Excoffier et al. (1992). We considered haplotype frequencies alone and by weighting the number of base-pair substitutions among haplotypes. Evolutionary relationships among haplotypes were estimated using distances generated under the Tamura-Nei model of sequence evolution (Tamura-Nei 1992). The neighbor-joining method of Saitou and Nei (1987) was used to generate the evolutionary tree. We next determined relationships among populations by constructing neighbor-joining phylogeographic trees using coancestry coefficient distances (Reynolds et al. 1983) generated using the AMOVA and PHYLIP programs. Harlequin ducks (n=16) sampled in Shemya, Alaska, were used as an outgroup for this tree.
Results. – Each population of Harlequin Ducks had three to eight alleles at each bi- parentally inherited locus (Table 1). Mean number of alleles per locus, and observed and expected heterozygosities were moderately high and concordant across all populations (Table 2). The bi-parental loci did not deviate significantly from Hardy-Weinberg equilibrium (i.e. there was no heterozygote deficit) and no evidence of linkage was observed in any population. This lack of heterozygote deficiency corresponds with the low and non-significant values of f (alleles within individuals variation) present for each locus and across all loci (Table 1). Inbreeding coefficients (F) were low when averaged across loci (Table 2) but were more variable when analyzed for each population and locus combination (range: –0.266 to 0.262; mean = 0.008; data not shown). These values were not significant in all cases except for the Katmai population for one locus (Sfi?4, P < 0.05). Estimates of spatial variation for the four bi-parental microsatellite loci were low (Table 1). No locus showed a significant difference in allele frequencies among populations within regions or between regions (all P>0.05). When averaged across all loci, allele frequencies also did not differ significantly among populations within regions or between regions (Table 1). Estimates of RST did not differ appreciably from FST spatial variation values, and also showed no significant difference in allele frequencies among populations within regions or between regions (data not shown). The lack of spatial variation among populations was further evident based on the low levels of bootstrap support for most nodes of the consensus neighbor-joining tree using the bi-parental loci (Fig. 2). Only three nodes within the microsatellite tree received >50% bootstrap support. The fact that clusters were composed of populations from both regions, coupled with the lack of consistent clustering at most tree nodes, suggests weak population structuring. However, our ability to place accurate and precise bootstrap values on the various nodes of the tree may be compromised by the relatively small number of bi- parental loci used (Takezaki and Nei 1996). Genetic diversity was high for the Bca?4 sex-linked locus where the number of alleles in any one population frequently exceeded ten (Table 2). The Sfi?1 locus had much lower genetic diversity scores. Like the bi-parentally inherited loci, there were no significant differences in allele frequencies for sex-linked loci among populations within regions or between regions (Table 1). The absence of population structure using the Bca?4 locus was potentially misleading, however, given the large number of rare alleles and the relatively low sample sizes for each population. Allele frequencies were generally very similar among populations for the Sfi?1 locus. One hundred twenty seven Harlequin Ducks from the PWS (n = 63) and APKA (n = 64) regions (Table 1) were sequenced to test for variation in the mtDNA control region. Six of the 163 base-pair sites were variable across all individuals; only two were present in more than one haplotypes, and all occurred as transition substitutions (data not shown). There were 53 nucleotide substitution differences between the numt and the mtDNA haplotypes. Forty-nine sites occurred as transition substitutions and the remainder were inversion/deletions (data not shown). This variation resulted in a total of seven unique haplotypes. Six of the seven haplotypes were found in both the PWS and APKA regions, and only one haplotype ("F" in Table 1) was detected in a single individual whereas the most common haplotype ("G" in Table 1) accounted for 52% of the samples. Haplotype diversities ranged from 0.427 to 0.745 (Table 2), with the highest diversity located in the Montague Island and Afognak populations. Haplotype diversity was not significantly lower in PWS (0.574) than in the APKA (0.702) region (Monte Carlo simulation, P = 0.08). Haplotype diversity values were highly concordant across populations, and were not significant among populations or regions based on Monte Carlo simulations (all P > 0.05). Nucleotide diversity corresponded closely to haplotype diversity. Neighbor- joining phylogenetic trees based on distances between haplotypes failed to support phylogeographic relationships (i.e., particular haplotypes were not restricted to specific sampling sites). Indeed, phylogenetic relationships (tree not shown) determined by the neighbor-joining analysis were only weakly supported by the bootstrapping analyses. Hierarchical analyses conducted without information on haplotype evolutionary relationships (i.e. number of nucleotide substitutions) revealed no significant differences in mtDNA haplotype frequencies (Table 1). When haplotype evolutionary information was included, similar results were produced, however a marginally significant difference in frequency of the two most common haplotypes ("A" and "G") was found among regions (P=0.049, see Table 1). Haplotypes differed in sequence by 1-2 base-pair substitutions (data not shown). The lack of structuring among populations was also evident when we constructed a neighbor-joining tree using the mtDNA haplotype frequency information (Fig. 3 – maybe combined with 2). Branch lengths were extremely short and it was not unusual for populations from the PWS and APKA to be clustered together. A comparison of the population consensus tree generated from the four bi-parental microsatellite loci (Fig. 2) and the branch length tree generated with mtDNA sequence information (Fig. 3 - maybe combined with 2) revealed several similarities. First, the Bay of Isles and Katmai populations, and the Uyak and Uganik populations clustered together in both trees. The Afognak population also clustered near the Uyak and Uganik populations in both trees, although this cluster differed in its relative position within each of the trees.
Discussion – Our analyses show that Harlequin Ducks molting and wintering within PWS and APKA share population genetic characteristics consistent with those expected for organisms belonging to one panmictic population (i.e. there is little or no spatial population structuring). Indeed, within regions, all three classes of molecular genetic markers indicated no significant differences in allele and haplotype frequencies among populations. Further, the topologies of neighbor-joining consensus trees and the short lengths of tree branches revealed few consistent population groupings. There was, however, a significant difference between PWS and APKA regions with the mtDNA marker (Table 1). Given the lack of information on evolutionary rates on microsatellite loci relative to the mtDNA, it is difficult to interpret these findings. Nevertheless, this regional difference was found only when the frequencies of mtDNA haplotypes were compared, i.e. there was a notable paucity of private alleles within populations. Consequently, we interpret these results as a general lack of historic structuring within and across regions. Our finding that Harlequin Ducks have little to no population structuring within PWS and APKA was surprising given that life history characteristics suggest that discrete, reproductively isolated populations exist. Indeed, males and females are highly philopatric to molting and wintering areas (Robertson 1997, D. Esler unpubl. data) and pair formation occurs on the wintering grounds during early to mid-winter (Gowans et al. 1997, Robertson et al. 1998). This disparity may be explained in four ways. First, Harlequin Ducks in PWS and APKA may represent a recent range expansion from a refugial population (Wenink et al. 1994). Second, barriers to gene flow may have only recently become established and insufficient time has elapsed for genetic differences to have evolved. Third, habitat alterations (e.g. Katmai Eruption, Rigg 1914; Great Alaskan Earthquake of 1964, Jacob 1986) may result in episodic dispersal and possible gene flow among populations which are otherwise presently reproductively isolated. And finally, low levels of adult or juvenile movement may occur between populations and regions. Recent band-resighting data within APKA and PWS regions support the idea that individuals occasionally move between populations (although apparently not between regions; D. Esler and D. Zwiefelhofer unpubl. data). If these emigrants mate with local individuals, they may provide sufficient gene flow to homogenize gene frequencies (Wright 1931). Indeed, subadults and juveniles are the most likely avenue for gene flow because they may remain in the marine coastal environment for several years during maturation (xxxxxxx). Lack of natal philopatry to molt and/or wintering sites may, in essence, nullify the effects of high adult fidelity to these same sites. Although the specific reason for the absence of genetic structuring among these populations is not known, our results indicate genetic variation was not lost as a result of the Exxon Valdez Oil Spill. This assumes, of course, that our genetic markers are representative of the overall genetic characteristics of the populations. The poor recovery of Harlequin Ducks in oiled areas does not appear to be related to Harlequin Duck population genetic characteristics, but may instead be associated with unfavorable local environmental conditions. Further, our inability to determine whether or not there is currently gene flow among regions or populations, makes it impossible to know if immigration will bolster population numbers in damaged areas within the spill zone. Even though genetic characteristics of the populations appear to be similar, translocations of Harlequin Ducks to augment populations which have not recovered seems unwise given the paucity of information on adult movements, breeding philopatry, mating strategies, and other basic natural history data on this species. Nevertheless, these genetic data do serve as an important baseline from which to assess the impact of future environment perturbations on genetic characteristics.
Acknowledgments – Funding for this project was provided by the Exxon Valdez Oil Spill Trustee Council (project no. 96161), the Alaska Biological Science Center of the U. S. Geological Survey, Biological Resources Division, Kodiak National Wildlife Refuge, Katmai National Park and Preserve, and the U.S. Fish and Wildlife Service Ecological Services Lafayette Field Office. Support during the final stages of data analysis and writing were provided by the Belgium Fund for Scientific Research, the Michigan Department of Natural Resources, and the Department of Fisheries and Wildlife at Michigan State University through the Partnership for Research and Management (PERM) program. J. Pearce helped test microsatellites for amplification and variability. J. Piatt provided unpublished information on Harlequin Duck mortality resulting from the Exxon Valdez oil spill. Office and field assistance was provided by G. Johnson, H. Brokate, B. Hobbins, T. Arensburg, B. Rice, C. Berg, L. Bennett, R. Leatherman, and A. Goatcher (APKA); and B. Baetsle, R. Ballas, B. Benter, T. Bowman, K. Burek, J. DeGroot, B. Jarvis, D. Mather, D. Monson, J. Morse, D. Mulcahy, D. Ruthrauff, D. Schaeffer, and K. Trust (PWS). J. Pearce, B. Greene, G. Mittelhauser, J. Rhymer, and G. Shields provided useful comments on analyses or earlier drafts of this paper.