DNA evidence for historic population size and past ecosystem impacts of gray whales

Genetic Variation and Mutation Rates at Multiple Loci.

We sequenced amplified gene segments for seven autosomal introns, two X-linked introns, and the mitochondrial marker cytochrome b from up to 42 individuals from the eastern Pacific gray whale population, and estimated substitution rates for these markers [see Materials and Methods and supporting information (SI) Methods, including SI Table 2]. The average rate of substitutions across autosomal nuclear introns was 4.8 × 10⁻¹⁰ substitutions per base pair per year⁻¹ (ranging from 1.5 × 10⁻¹⁰ to 10 × 10⁻¹⁰) (Table 1).

We used the coalescent analysis program LAMARC (12) to estimate genealogies from individual sequences and calculate the genetic diversity parameter, ϴ = 4N_eμ, where N_e is the effective population size and μ is the average mutation rate. We combined data from all loci into a joint-likelihood analysis. The overall maximum-likelihood point estimate of ϴ was 0.001021, with 95% confidence intervals ranging from 0.000925 to 0.001130.

Variance of Effective and Census Population Sizes.

To calculate the effective population size N_e, we divided the joint maximum likelihood estimate of ϴ by estimates of generation time ranging from 15.5 to 22.28 years (13, 14), and by four times the average autosomal substitution rate (μ), after applying scaling factors described in Materials and Methods. Generation time range was calculated as the median age of 54 sexually mature females (13) and as the mean period elapsing between the birth of a parent and the birth of offspring (14).

We calculated the average long-term effective population size of gray whales to be on average 34,410 with 95% confidence limits of 31,175 and 38,084. However, census size of animal populations is typically higher than effective size because not all adults successfully breed. We converted effective size N_e into total census estimates (N) by multiplying by a conservative 2:1 ratio of total adults to breeding adults (15) and the ratio of total population size to total adults, estimated between 1.58 and 1.78 (13, 16) based on census and fisheries data.

By using the relationship ϴ = 4N_eμ, and the conversion factors above, we computed 95% confidence intervals on census size N by randomly sampling 10,000 times from uniform distributions of ϴ (0.000925–0.001130), generation time (15.5–22.28 years), and the ratio of census population size to effective population size (3.2–3.6). This procedure gives 95% confidence limits of 78,500–117,700 with a mean of 96,400 (Fig. 1), or 3.5–5.3 times today's census population size. This range of values incorporates uncertainty in measures of genetic diversity, uncertainty in mutation rates among loci, and uncertainties in generation time and juvenile abundance.

Testing for a Prewhaling Population Bottleneck.

Genetic data provide population size estimates potentially averaged over thousands of generations or more. Thus, the genetically estimated gray whale population size might be higher than at the start of commercial whaling if a large population decline occurred before the mid-1800s. Tests for departure from neutral distributions of alleles within populations and segregating sites (Tajima's D, Fu and Li's D*, Fu's F) (SI Table 3) showed no evidence for loci under selection, or for significant population growth or decline. A separate way to estimate past population dynamics is to use coalescent analysis (12, 17). However, the low average mutation rate across whale introns and the large gray whale population size lengthen the time period over which our data provide useful views of population changes. Coalescent analyses performed on our intron data set show no long-term growth or decline, but lack the power to detect relatively recent population swings.

However, loci with relatively high mutation rates and more sensitivity to genetic drift are more useful for testing for changes in population size, so we used current mitochondrial DNA haplotype diversity in gray whales to examine the effects of different bottleneck scenarios. We carried out coalescent simulations of bottleneck events at a variety of times before the onset of commercial whaling, to test whether such an event could explain differences between genetic and historic estimates (see Materials and Methods).

We simulated a prewhaling census size change from 96,000 to 22,000, approximately the difference between genetic estimates and average census size today, varying the bottleneck time from 15 generations to 10,000 generations ago. We then tested for significant reductions in mitochondrial haplotype diversity in simulated vs. current populations. These tests show that haplotype diversity (Hd) significantly declines for all bottlenecks occurring longer ago than 73 generations, or ≈1,100–1,600 years ago (Fig. 2). These analyses are therefore inconsistent with the hypothesis that a bottleneck from 96,000 to 22,000 animals further back in time than 1,100–1,600 years can explain the difference between our genetic estimate and prior estimates of historic abundance: such a bottleneck would have eroded haplotype diversity well below the current value.

Quantifying the Effects of Population Structure.

A final consideration is population structure, because genetic diversity among subpopulations can inflate N_e. In general, estimates of effective population size exceed real population size in this situation by a factor of ≈1/(1 − F_ST), where F_ST is the proportion of genetic variance distributed spatially (ref. 18, equation 2.7 based on ref. 19). Available data show no observable structure within the eastern Pacific gray whale population (reviewed in ref. 20). To further confirm this result, we collected data from six microsatellite loci for the same whales used in this study and estimated the probabilities of one population (K = 1) versus multiple populations (K > 1) by using STRUCTURE, version 2 (21). No population subdivision is apparent in our data, suggesting that our estimate of genetic diversity is not inflated by current population substructure.

However, two other populations may have contributed genetic diversity to the eastern Pacific population in the past: the western Pacific population (16) and an extinct Atlantic population (22). A recent study comparing mitochondrial control region data between the eastern and western Pacific gray whale populations indicates an F_ST value of 0.087 between these two populations (23), potentially inflating our estimate of population size in the eastern population by about a factor of ≈9.5% [1/(1 − 0.087) = 1.095]. However, this formulation assumes the populations are at migration–drift equilibrium, which is unlikely to be true for gray whales. To capture the complex contribution of potentially nonequilibrium migration scenarios to current-day genetic diversity, we constructed a series of population simulations to estimate the impact of migration from partially isolated western and Atlantic gray whales on genetic diversity in eastern gray whales. Our basic question was whether periodic migration from small satellite populations could significantly increase genetic diversity in a central population.

We simulated a set of scenarios likely to be realistic given the appearance and disappearance of biogeographic barriers over the past 200,000 years (SI Fig. 4). Migration between the eastern and western Pacific populations is allowed during the last glaciation (18–70 kya), and the Little Ice Age (400–750 ya) when sea ice or lower sea level may have blocked migration through the Bering Sea (24). Arctic migration between the eastern Pacific and the Atlantic is assumed to have been possible only during interglacial, warm periods, most recently during the Sangamonian Interglacial period (114–131 kya).

Simulations with a variety of migration rates were conducted by using effective population sizes as input values, but in the following, we converted these effective sizes to estimated census size by using the conversion factors described above. Simulations of an eastern population with an effective size of 11,500 (corresponding to a census size of 40,000) and small populations in the western Pacific and Atlantic (census size, 3,000–6,000 each) show no impact on the genetic diversity of the eastern Pacific (Fig. 3). Large satellite populations, however, can increase the genetic diversity of the eastern population. Of course, at these large sizes (approximate census size of 30,000) (Fig. 3), the western and Atlantic populations contain a substantial fraction of the global population. Simulations also indicate that increased diversity in the eastern population is much more likely to have been caused by immigration from the western Pacific population than the Atlantic. Gene flow from a large Atlantic population alone increases genetic diversity values in the present eastern Pacific population by <10% (data not shown).

Sources of Uncertainty.

These estimates involve two major levels of analysis, the estimate of effective size from genetic diversity and the estimate of census size from effective size. In the use of genetic diversity to measure effective size, genetic variation and subsequent estimates of long-term effective population size differ among loci, and are affected by many evolutionary forces, such as selective sweeps, population expansions and bottlenecks, variance in reproductive success, stochastic retention of genetic diversity, population structure, and mutation rates. Although most of these factors are likely to decrease estimates of effective population size, stochastic effects, mutation and population structure can potentially increase it. As in other genetic studies of wild populations, such as chimpanzees (25), we have attempted to control for these three factors by measuring effective size across independent loci, estimating mutation rate independently for each locus, and by accounting for population structure.

Because data from multiple unlinked loci independently assay past population patterns, analyzing large numbers of loci greatly increases the accuracy of estimates of ϴ, even when the number of individuals sampled is moderate (26). Among the nine polymorphic loci we examined, genetically based population census estimates vary from 44,000 to 272,000, but exceed current census estimates in all cases. Only in the X-linked intron G6PD, in which no variation was observed, is the genetic estimate of long-term abundance less than the current population size (Table 1).

Different loci potentially have different mutation rates, so we independently measured mutation rates for each locus. Prior genetic measurement of whale populations relied on variation in the mitochondrial control region, which has a complex pattern of molecular evolution (4). In contrast, for whale introns the relationship between divergence time and genetic distance is largely linear (SI Figs. 5 and 6). The most critical fossil date used to calibrate substitution rates, the divergence between right whales and other mysticetes, is well supported. This time point represents the earliest well dated split within mysticetes and thus has a large impact on the relationship between genetic distance and divergence. Fossil dates (26–30 Mya) for this split are concordant with molecular clock estimates (27.3 ± 1.9 Mya) (ref. 27 and ††), suggesting this time point is a robust anchor for our mutation calibrations. Our calculated rate of intron evolution is similar to that of other mammals after taking into consideration the large size and low metabolic rate of whales (28).

Because the amount of genetic data varies between markers, our average rate might differ slightly if we calculated a rate by using data from all loci in the same analysis. This approach generates a slightly lower rate [4.15 × 10⁻¹⁰ ± 0.3 × 10⁻¹⁰ (SE) bp⁻¹ year⁻¹] (SI Fig. 6). Because a lower mutation rate will result in higher estimates of N_e, we used the higher rate obtained from averaging across loci in our analysis (4.8 × 10⁻¹⁰ bp⁻¹ year⁻¹). In our LAMARC analysis, the estimate of ϴ for each genetic marker is scaled by that marker's relative mutation rate. This allows us to use an average mutation rate to calculate overall census size estimates from LAMARC's output.

Ecological counts include nonreproductive adults and juveniles observed from shore, so to compare our results to published census data we must correct for these parts of the population. The ratio of adult population to effective population has been estimated as >2 for many mammals (15): models of the relationship between population size and extinction risk, including those for whales, often assume ratios of 2–10 (29). In the absence of a direct estimate of this parameter from lifetime measures of reproductive success in males and females, we assumed the most conservative ratio of 2. Juvenile abundance can be more confidently estimated from census data (13, 16).

We estimated overall variance by incorporating uncertainty in both the generation of effective size from diversity and the generation of census size from effective size through a Monte Carlo sampling method. This approach gave us confidence intervals that ranged from 76,000 to 118,000. Other sources of uncertainty exist, such as (i) variance in reproductive success among adults, (ii) impact of selection at individual loci, and (iii) swings in population size over time. However, each of these factors, if known in detail and taken into account, would tend to increase our estimates of the ratio of census to effective population size, and would increase our calculated long-term census population size by some unknown amount.

Long-Term Nature of Population Estimates.

Genetic data provide long-term estimates of average effective population size, and there are many population trajectories throughout the past that could give rise to a particular effective population size today. Until reconstructions of whale population size in different past time periods become available, the single value of long-term size we provide should be treated cautiously. In particular, it is possible that the population sizes of whales just before whaling might have been lower than their long-term sizes. Although no direct evidence is available to evaluate this hypothesis, such a population trajectory could explain the discrepancies between historical and genetic views of whale abundance (4).

We explored this possibility with simulations of gray whale population bottlenecks from our estimated census size of 96,000 down to 22,000 individuals. Simulations show that such a decline would rapidly erode mitochondrial haplotype diversity. Current data would be able to detect this erosion if the prewhaling bottleneck had happened more than 1,100–1,600 years ago (73 generations).

Two possibilities remain. First, a decline in gray whale numbers from 96,000 to 22,000 might have occurred over the past 1,100–1,600 years. If this happened before western whaling began in the 19th century, then both long-term genetic estimates and traditional estimates of gray whale abundance just before whaling could be correct. A second possibility is that a population much larger than 96,000 collapsed to 22,000 further back in time. Analysis of gray whale genetic diversity from ancient samples may clarify these possibilities, and allow the exploration of other possible population trajectories.

The Influence of Population Structure on the Results.

Although there are many demographic events that reduce genetic variability, there are a few that increase it. Balancing selection can increase diversity, but is unlikely to be operating across all of the multiple loci we analyzed, and there is no signature of excess high-frequency alleles in our data (SI Table 3). More likely might be that undetected, nonequilibrium population migration could increase the apparent genetic diversity of whale populations. Accounting for past population structure is complex, given a potentially infinite number of past migration scenarios. We simulated a set of scenarios in which genetic diversity might have been injected into the eastern Pacific population from the western Pacific population when it was larger, or from the Atlantic before that population became extinct. Our simulations suggest that small populations in the western Pacific or in the Atlantic would not have appreciably increased genetic diversity in the eastern population. However, historical migration from a large western Pacific population could partially explain high genetic diversity in the eastern Pacific.

As a result, our long-term estimate of 96,000 gray whales was probably distributed between the eastern and western populations. Very few data are currently available on historical abundance and whaling effort in the western Pacific stock (20), but it is possible that this population was substantially larger than it is now. Future analysis of these two populations using multiple nuclear loci and the IM approach used by Won and Hey (25) to investigate ancestral and effective sizes in chimpanzees might be a profitable way to address these issues.

Implications for Management.

Conservation efforts have resulted in a steep decline in extinction risk for eastern Pacific gray whales, but our results suggest that full demographic recovery has not occurred and that this population remains depleted. Mean census counts of ≈22,000 (8) are 28–56% of the minimum genetic long-term estimates, assuming up to half the population have typically been western. Because marine mammal populations <50–70% of historic population size are considered depleted (7), a cautionary response to these data would be to consider the eastern Pacific gray whale as a depleted stock. Such a designation could change the allowable human-caused mortality of this stock under the U.S. Marine Mammal Protection Act from ≈417 animals a year to ≈208, assuming the recovery factor (7) used in calculating potential biological removal (7, 8) changes from 1.0 (8) to 0.5 (7, 30). Currently, an aboriginal take of 125 animals per year is allowed by the International Whaling Commission (8, 30).

For gray whales, historic population data help clarify population trends. Models that incorporate census data from 1999–2000 indicate that gray whales have reached their carrying capacity (7), but other interpretations suggest that high mortality in 1999–2000 represented an ephemeral episode due to climatic shifts (9). Genetically determined long-term population estimates support the hypothesis that the current population has not fully recovered and will continue growing if ecological conditions permit.

An important question for current whales is whether carrying capacity has declined over time. If this is the case, gray whales may be depleted relative to historical numbers but may have reached carrying capacity today. Ecological surveys of gray whale feeding areas on the Bering Sea shelf suggest that this area alone could support ≈90,000 whales annually (see SI Methods). However, recent evidence suggests that gray whale feeding habitat may be declining as Arctic benthic prey populations are reduced because of changing climate in the Bering Sea (9). Although additional survey data will be critical to determining whether carrying capacity has been reached, our estimate of typical gray whale abundance suggests that recent problems in gray whale feeding, including reports of thin adults or high calf mortality, may result from changing conditions in northern feeding grounds.

Ecosystem-Wide Impacts of Gray Whale Depletion.

In addition to implications for management, these data are a first step toward quantification of the ecosystem effects of whale population depletion in the North Pacific. Gray whales are important ecological structuring agents in Bering Sea benthic marine communities (31, 32). Because they are bottom feeders that suck up mouthfuls of sediment, study of feeding gray whales shows ≈1.2 × 10⁸ m³ of sediment were annually resuspended by the eastern Pacific gray whale population of the early 1980s (33, 34). Assuming a population size then of ≈16,000 individuals, a population of 96,000 gray whales would rework ≈7.2 × 10⁸ m³ in a summer, ≈12 times larger than the sediment transport load of the largest river emptying into the Bering Sea, the Yukon River (35, 36). Decreased sediment reworking could dramatically change nutrient recycling, and create shifts in benthic species dominance (32).

Similarly, feeding by gray whales provides nutrient subsidies from benthic marine communities to terrestrial ones, including food subsidies for at least four species of seabirds that feed on benthic crustaceans brought to the surface by gray whale feeding (37). The number of birds attending the foraging activities of a single gray whale averages ≈11, although the number of plumes used per bird and the use of plumes from different whales are not known. If this average is representative across the summer season, we calculate that a population of 96,000 whales could provide food subsidies to 1.03 million birds. The quantitative impact of these subsidies on sea bird reproduction or fitness is currently unexplored. In addition, gray whales may have provided an important food source for predators and scavengers such as orcas (38) and California condors (39).

More numerous gray whales in the past may not have fed only on the Bering and Chukchi shelves, the areas in which the majority of individuals feed today. As gray whale populations have increased, more and more are observed feeding in other coastal locations including Oregon, Washington, British Columbia, and southeast Alaska (e.g., ref. 41). If feeding in other areas was common, the ecological impact of diminished gray whale populations would not be restricted to the Arctic.

Overall, these simple calculations of gray whale impact do not fully document the likely ecological impact of whales; they merely suggest the order of magnitude of effects of gray whales on their environment. Further research on ecological subsidies of whales to seabirds, sediment and nutrient effects, and ecological shifts of whales during population expansion are critically needed to better understand how marine ecosystems have been impacted by reduced whale populations.

Estimating Genetic Population Parameters.

The diversity parameter ϴ was calculated by using LAMARC (12), by jointly estimating ϴ across multiple loci while incorporating relative μ and N_e values for each locus. Among genomic regions (seven autosomal nuclear introns, two X-linked introns, and cytochrome b), we applied scaling factors of 4 and 4/3 to our mitochondrial and X chromosome N_e values, respectively, and computed ϴ by using a scale in which the average relative mutation rate among autosomal loci is unity. These two sets of scaling factors allowed us to calculate the population's ϴ from a multidimensional likelihood surface constructed from all loci. To ensure model convergence, we performed each LAMARC analysis 15 times, using three different random number seeds and five different trial values for ϴ. For each gene, we performed 10 MCMC searches of 31,000 iterations each followed by two searches of 1,001,000 iterations each, sampling every 20th genealogy estimate after discarding the first 1,000 estimates of each search.

Substitution rate was estimated by comparing pairwise genetic distance between 4 and 14 species of baleen whales and their respective divergence times (27). We used the program MODELTEST (42) to determine the appropriate mutational model, and PAUP* to generate pairwise genetic distances (43).

Bottleneck Simulations.

We used SIMCOAL (44) to simulate a sequence of 523 bp for 42 individuals, the length of the control region segment we obtained, and used a mutation rate of 5 × 10⁻⁸ bp⁻¹ year⁻¹, which is the mutation rate necessary to obtain the same average effective size of eastern North Pacific gray whales for control region data as we found for intron data. We used an empirically determined transition bias of 0.96 and a rate heterogeneity parameter of 0.8. The empirical likelihood of haplotype diversity (Hd) under each demographic model was calculated as a one-tailed test [see Belle et al. (45)] and decreases as the time of bottleneck becomes more ancient (Fig. 2). Linear regression shows that a likelihood of 0.01 is obtained under a bottleneck occurring 73 generations ago.

Migration Between Subpopulations.

To test the effect of immigration from satellite populations on effective population size of females, we used SIMCOAL (44) to simulate a central population with migration of varying levels from satellite populations of varying sizes, and calculated the ratio between N_e(f) of the central population computed with and without migration. The central population was assigned an N_e(f) of 5,000 individuals, corresponding to the estimated effective size of the female population today. We simulated immigration rates of 0.001, 0.01, and 0.1 (proportion of immigrants per generation) and satellite populations of sizes 500, 1,000, 2,500, and 5,000.