Climate change during the late Quaternary period (LQP), in both temperate and tropical zones, has been a major driver in the shaping of species distributions and abundances. Our understanding of the various effects of climate change on population dynamics of marine species in temperate zones (such as the Gulf of California) is growing. However, studies on the demographic history of seabirds are rare and there is no description of how regional climate change has affected high-trophic-level marine species such as Heermann's Gull (Larus heermanni; Charadriiformes: Laridae). We investigated whether the demographic history of Heermann's Gull reflects population change consistent with past changes in climate in the Gulf of California during the LQP. We also explored whether those past changes affected the demographic history of codistributed marine organisms in a similar way as found for Heermann's Gulls. To test our hypotheses, we sequenced a fragment of the mitochondrial DNA cytochrome b gene (1,033 base pairs) and performed tests to investigate whether demographic change occurred within the LQP. The results of the 3 approaches used were consistent with a historical demographic expansion during the LQP. All analyses (Fu's FS test, Tajima's D neutrality test, mismatch distribution analysis [MDA] and associated demographic parameters, and Bayesian skyline plots [BSP]) were consistent with a model of population expansion in Heermann's Gulls. The MDA estimated the expansion event at ∼48,000 yr before present (yr BP; 95% confidence interval: 34,000–72,000), whereas the BSP showed that population growth began ∼100,000 yr BP and lasted until ∼45,000 yr BP. We discuss possible associations between the demographic expansion of this seabird species and large-scale ecological shifts or demographic expansions of other marine species.
INTRODUCTION
In order to understand contemporary biological diversity and species distributions, it is necessary to recognize the influence of recent and historical climate change on these factors, at both global and local scales. The glacial–interglacial cycles of the late Quaternary period (LQP) have played a major role in shaping the distribution and abundance of contemporary species and populations, and in the extinction of others (Hewitt 2004). The study of these factors has proved useful in recovering information of postglacial demographic histories, speciation, and divergence (Hewitt 1996, 2000, Ramírez et al. 2013). Some studies have compared intraspecific genetic patterns among several taxa within the same geographic area, in order to search for congruent patterns of genetic variation (Brunsfeld et al. 2007). When such congruence has been found, it is considered an indicator of the influence of common historical factors (Soltis et al. 1997, Petit et al. 2003, Brunsfeld et al. 2007).
Changes in climate during the LQP in the central region of the Gulf of California (particularly the Guaymas Basin) have been studied. Major changes in ocean productivity off both coasts of the Baja California peninsula, on the northwestern Mexico mainland, and in the southwestern United States have been documented for periods ranging from 140,000 years before present (yr BP) to contemporary times (Cheshire et al. 2005, Pérez-Cruz 2013, Price et al. 2013). Environmental oscillations in the Gulf of California have likely promoted changes in the demographic histories, genetic structuring, and genetic diversity of individual species as well as in the composition of marine communities over the past 200,000 yr to present times. Some studies have set the time of population expansion at 190,000–130,000 yr BP (Pfeiler et al. 2008, Liu et al. 2011), others at 71,000–50,000 yr BP (Pfeiler et al. 2005, Hurtado et al. 2007, Díaz-Viloria et al. 2012). These studies have focused almost exclusively on species that occur strictly in marine environments, such as fishes, mollusks, and crustaceans; one study considered a nonaquatic vertebrate (the fish-eating bat Myotis vivesi; Mejía et al. 2011), which, though not inhabiting the ocean, lives on oceanic islands and feeds almost exclusively on marine fishes and crustaceans (Maya 1968, Otálora-Ardila et al. 2013). Although studies of seabird genetics exist for other geographic regions (Milot et al. 2007, Dantas et al. 2012, Welch et al. 2012), little is known about the effects of regional climate change on seabird demographic history in the Gulf of California. By focusing on Heermann's Gull (Larus heermanni; Charadriiformes: Laridae), the present study will shed light on how climate change has affected high-trophic-level species in this region.
Several studies have examined the biology and ecological relationships of Heermann's Gull, including adaptation to the hot, arid environment in which it nests (Bartholomew and Dawson 1979, Bennett and Dawson 1979, Rahn and Dawson 1979); interactions with predators and breeding biology (Velarde 1992, 1993, 1999); effects of parental age and food availability on breeding success (Vieyra et al. 2009); and how its diet composition can be used to predict forage fish abundances (Velarde et al. 2004, 2013, 2015a). However, the effects of past climate change on this species' population size through time, and whether changes in its population size are at all related to population changes in other marine species in the region, have been unknown.
Here, we investigate whether the demographic history of Heermann's Gull reflects population change consistent with past changes in climate in the Gulf of California during the LQP and with population change in other species that occur in the area. We hypothesized that if past changes in climate affected the primary productivity and, consequently, the demographic histories of both Heermann's Gull and other codistributed marine organisms in similar ways, then the periods when this occurred would also be similar.
METHODS
Study Population
Heermann's Gull has been placed under special protection (Pr) by Mexican federal environmental law (SEMARNAT 2010) because most of the population (∼95%) nests on a single island, Isla Rasa (28.82°N, 112.98°W; Mellink 2001) (Figure 1). Recently the species has been reported to nest in 14 sites, with Isla Rasa harboring ∼130,000 pairs; Isla Cardonosa (4 nautical miles north of Isla Rasa), with ∼4,000 pairs, is the second-largest colony (Velarde 1999, Mellink 2001, Velarde et al. 2005). The species feeds on small pelagic fishes and is an important part of the regional trophic web (Velarde et al. 2004, 2013, 2015a).
Collection of Blood Samples, DNA Purification, Amplification of Cytochrome b, and Sequencing
Samples were collected from Isla Rasa and Isla Cardonosa (28.89°N, 113.03°W). Gulls were captured using walk-in traps in 2011 and 2012, and 60 μL of blood was drawn from the brachial vein and placed in Longmire buffer (Longmire et al. 1988). Samples were stored at –20°C until DNA was extracted using the DNeasy Tissue Kit (Qiagen, Valencia, California, USA) following the manufacturer's protocol. Polymerase chain reaction (PCR) amplification of mitochondrial DNA cytochrome b gene (Cyt-b, 1,033 base pairs) was performed with the primers L14841 (5′ AAA AAG CTT CCA TCC AAC ATC TCA GCA TGA TGA AA 3′) and H16065 (5′ GGA GTC TTC AGT CTC TGG TTT ACA AGA C 3′) (Crochet et al. 2003), using an ABI 2720 Thermal Cycler (Applied Biosystems, Foster City, California, USA). The PCR reaction contained the following: 1× PCR buffer (10 mM Tris–HCl, 50 mM KCl), 1.5 mM MgCl2, 10 mg mL−1 BSA, 0.25 mM of each dNTP, 2U of Amplitaq (ABI), and 20 ng of genomic DNA. Cycling parameters were as follows: initial denaturation for 5 min at 94°C; 35 cycles of 1 min at 94°C, 1 min at 54°C, and 1 min at 72°C; and a final extension of 5 min at 72°C. Amplified products were purified using the Qiaquick PCR Purification Kit (Qiagen) and sequenced at the University of California Berkeley Core Facility on an ABI 3730 DNA Analyzer (Thermo Fisher Scientific, Waltham, Massachusetts, USA). Chromatogram files were edited manually using Finch TV 1.4.0 (Geospiza, Seattle, Washington, USA; http://www.geospiza.com) and assembled with Seaview 4.5.4 (Gouy et al. 2010). Sequences were deposited in GenBank (accession nos. KT275317–KT275602).
Analyses of Genetic Variation and Demographic History
We estimated the haplotype richness (h), number of segregating sites (S), average number of nucleotide differences (k), haplotype diversity (H), and nucleotide diversity (π) (Tajima 1983, Nei 1987) using Arlequin 3.11 (Excoffier et al. 2005). A parsimony network was produced using TCS 1.21 (Clement et al. 2000), and no insertions or deletions were considered as a fifth character state, given that they were not found in the dataset. In order to detect departures from the mutation–drift equilibrium, Fu's FS (Fu 1997) and Tajima's D (Tajima 1989) were estimated with DnaSP 5.10.00 (Rozas et al. 2003) and Arlequin.
To examine the demographic history of Heermann's Gull, 3 different methods were used: (1) neutrality against population growth, (2) mismatch distribution analysis (MDA), and (3) Bayesian skyline plots (BSP). We examined the mismatch distributions (Harpending 1994, Schneider and Excoffier 1999) using Arelquin. Estimated parameters were as follows: expansion time τ = 2ut (where u is the cumulative substitution rate and t is the number of generations since the expansion), θ0 = 2 μN0, and θ1 = 2 μN1 (where θ is the scaled mutation rate, μ is the mutation rate, and N0 and N1 are the effective population sizes before and after the expansion). To test the validity of the demographic model, we calculated Harpending's (1994) raggedness index and the distribution of the sum of squared deviations (SSD) between the observed and an estimated mismatch distribution. For all these parameters, we also obtained 95% confidence intervals by means of 10,000 randomizations of the data. We investigated the time since the presumed expansion event from τ, using substitution rates of 1.8% Ma−1 (Nabholz et al. 2009) and 2.0% Ma−1 (Quinn 1992). Time-since-expansion values were estimated using the spreadsheet tool for the application of the mismatch analysis (Schenekar and Weiss 2011).
Finally, historical population dynamics were also estimated through the BSP approach (Drummond et al. 2005), implemented in BEAST 1.8.0 (Drummond and Rambaut 2007). We employed a strict molecular clock using the aforementioned substitution rates and identified HKY (Hasegawa-Kishino-Yano) as the best-fitting model of nucleotide substitution for each codon position using Akaike's Information Criterion with Modeltest 3.7 (Posada and Crandall 1998). Markov chain Monte Carlo sampling was run for 50 million generations and sampled every 1,000 generations, with the first 5,000 sampled events discarded as burn-in. All estimators of effective sample size were >1,000. The results of 3 independent chains were combined in LogCombiner 1.5.4 (Drummond and Rambaut 2007), and the BSPs were generated in Tracer 1.5.0 (Rambaut and Drummond 2009).
RESULTS
Genetic Diversity and Neutrality
We obtained sequences of the mtDNA cytochrome b gene (aligned and edited length of 1,033 base pairs) from 266 individuals from Isla Rasa and 20 from Isla Cardonosa, yielding a total sampling size of 286. Sequences were highly conserved, with only 41 variable sites across 35 haplotypes (Table 1). Two of the haplotypes were common (H1 was observed in 165 individuals and H3 was observed in 70 individuals), 6 were more rare (H9, H11, H13, H14, H27, and H29), and the remaining 27 haplotypes were found in only one individual. The parsimony network (Figure 2) revealed that most haplotypes are separated from the most common haplotypes (H1 and H3) by a single mutation; 4 (H2, H6, H16, and H19) are separated by 2 mutations; and 3 (H4, H12, and H31) are separated by 3 mutations. The most common haplotypes (H1 and H3) are separated from each other by a single mutation.
TABLE 1.
Nucleotide polymorphism and results of neutrality tests for Larus heermanni mtDNA cytocrome b (1,033 base pairs) haplotypes.
Notes: n = sample size; h = number of haplotypes; S = number of segregating sites; k = average number of nucleotide differences; H = haplotype diversity; π = nucleotide diversity; FS = Fu's FS statistic; D = Tajima's D; SSD = sum of squared deviations; raggedness = Harpending's raggedness index; *P < 0.001, ns = nonsignificant.
FIGURE 2.
Parsimony network showing the relationship among Heermann's Gull mtDNA cytochrome b haplotypes. The size of the circles is proportional to the frequency in the dataset. The number of individuals sharing a haplotype is shown in parentheses; when only one individual was found with the haplotype, no number is shown. Empty circles represent inferred (but not sampled) haplotypes.
Haplotype richness, haplotype diversity, and nucleotide diversity values are shown in Table 1. Fu's FS and Tajima's D-tests for neutrality were both negative and highly significant (Table 1).
Demographic History of Heermann's Gull from Late Quaternary to Present
The mismatch distribution for all samples was unimodal (Figure 3). The demographic parameters estimated by the MDA were consistent with a null model of population expansion and were supported by nonsignificant SSD and Harpending's raggedness index (Table 1). The estimated time of population expansion was dated to 48,000 yr BP (95% confidence interval: 34,000–72,000) or 43,000 yr BP (31,000–65,000), using substitution rates of 1.8% or 2.0%, respectively. Finally, the BSP using the 1.8% substitution rate showed that the Heermann's Gull population appears to have experienced continuous population growth beginning during 100,000–45,000 yr BP (1.8%) or during 95,000–35,000 yr BP (2.0%) (Figure 4). In both cases, a relative constant population size can be seen from 25,000 yr BP to the present, even with the small differences between substitution rates.
FIGURE 3.
Demographic history of Heermann's Gull, estimated from mitochondrial cytochrome b gene. The plot depicts the distribution of pairwise numbers of nucleotide differences between individuals in the study population. The observed distributions (bars) are compared for their goodness-of-fit to a Poisson distribution under a sudden-population-expansion model illustrated by the curved line.
FIGURE 4.
Bayesian skyline plot of Heermann's Gull population, detected with mitochondrial cytochrome b gene. The fluctuation rate of effective population size over time is shown. The middle line (black) is the median estimate; upper and lower limits (95% highest posterior density) are also given (gray lines). The x-axis represents a population fluctuation rate of 1.8% Ma−1.
DISCUSSION
The present study represents an important effort to describe the demographic history of a seabird species (sampling size = 286). In this context, assessing Heermann's Gull genetic diversity is a first step toward understanding future trends for population conservation. For example, our estimated cytochrome b nucleotide diversity (0.00083) may be useful for comparison with similar studies (also using the same gene) that consider conservation status (see Lawrence et al. 2008 and references therein). On one hand, both the Magenta Island Petrel (Pterodroma magentae; Lawrence et al. 2008) and the Great Bustard (Otis tarda; Pitra et al. 2000)—two critically endangered or vulnerable species, as recognized by the World Conservation Union (IUCN, 2016)—have been characterized with a nucleotide diversity of 0.0013. On the other hand, species with a conservation status of “least concern” (IUCN, 2016), such as the European Storm-Petrel (Hydrobates pelagicus; Cagnon et al. 2004) and the Ferruginous Pygmy-Owl (Glaucidium brasilianum; Proudfoot et al. 2006), showed a nucleotide diversity of 0.004 (again, the same for both species). Although there seems to be no clear tendency in number of haplotypes in relation to IUCN conservation status (given that its estimation is highly dependent on sampling size), in the case of the value of nucleotide diversity, Heermann's Gull is more similar to species with a higher (more endangered) IUCN status. Heermann's Gull has been given special protection under Mexican federal norms but is categorized as “least concern” by IUCN. No conclusions are intended here regarding the relationship between conservation status and factors such as number of haplotypes or nucleotide diversity; however, our current dataset does not support a recent population bottleneck, in that we found evidence of a historical population expansion. In this sense, additional analyses with variable nuclear markers should shed light on more contemporary processes.
The highly significant negative Tajima's D values may represent evidence of purifying selection or evidence of a historical population expansion in Heermann's Gull. The population-expansion interpretation is also supported by the highly significant negative value of Fu's FS, a more sensitive indicator than Tajima's D. Evidence from MDA and BSP indicates that Heermann's Gull experienced a demographic expansion between 100,000 and 45,000 yr BP, a period that roughly coincided with the end of the Pleistocene glacial maxima to the last glacial retreat. The signature of major population expansion during the late Pleistocene for Heermann's Gulls is likely related to large-scale ecological shifts in the North Pacific and Gulf of California due to historical climate change and the subsequent demographic expansions of associated marine species, including species that are presently important components in the diet of Heermann's Gulls, such as the northern anchovy (Engraulis mordax; Velarde et al. 1994, 2013, 2015a).
During the late Pleistocene, the North Pacific developed into a cold, rocky shore ecosystem that provided complex and productive environments for many marine species, and dynamic ocean currents and upwelling that led to higher plankton productivity (Paine 2002, Graham et al. 2003, Lyle et al. 2010). Coastal areas in the North Pacific are highly variable environments, where changes in current patterns, ocean circulation, and upwelling have reversed in different time scales, ranging from 10,000 to 100,000 yr BP (Kennett and Ingram 1995, Kiefer et al. 2001, Hewitt 2004), affecting productivity (Webb and Bartlein 1992). Similar patterns of abrupt variations in climate and major changes in ocean productivity have been documented for the Gulf of California (Cheshire et al. 2005, Pérez-Cruz 2013, Price et al. 2013).
Numerous studies indicate that Quaternary climate change had a significant impact on the demographic history of marine species in both the northeastern Pacific and the Gulf of California (Table 2). A survey of 10 studies focused on this geographic region indicated 2 general periods of expansions, one before 130,000 yr BP and the other after 100,000 yr BP. Studies that included Pacific sardine and northern anchovy (main prey items of Heermann's Gull) found differential evidence for early vs. late population expansions. Lecomte et al. (2004) and García-Rodríguez et al. (2011) reported more ancient expansions of these 2 species, while Díaz-Viloria et al. (2012) suggested that northern anchovies expanded more recently. The latter result is more in line with the timing of expansion we have reported here for Heermann's Gull. While not the focus of our study, the differences reported for these species across studies could be attributed to different molecular markers and molecular clocks employed by various authors.
TABLE 2.
Summary of species that display population expansions from the eastern Pacific–Gulf of California. Studies are ordered chronologically, starting with more recent expansions.
Notes: cyt b = cytochrome b; CR = control region; COI = cytochrome oxidase I; yr BP = yr before present.
The fact that studies on 5 additional species, aside from Heermann's Gull, show evidence of expansions during the late Pleistocene (100,000–30,000 yr BP) indicates that regional oceanographic change during that period affected multiple marine species in a similar manner. The results of Díaz-Viloria et al. (2012) on the northern anchovy are most relevant to our results, given that long-term studies on the diet of Heermann's Gull have shown that northern anchovy is an important prey species (Velarde et al. 1994, 2013, 2015a). The population expansion of northern anchovy (Díaz-Viloria et al. 2012) could have been a key factor in the increased breeding success of Heermann's Gull and its own demographic expansion. Even though these results are correlative, they suggest a role for long-term predator–prey dynamics in the regulation of seabird populations from the Gulf of California. Population expansions of sympatric species may follow one another quite rapidly, particularly in the case of predator–prey systems, as suggested by the results of the present study, and future work that utilizes genome-wide markers may provide additional insights into these studies.
Conclusion
Although based on rough estimates, our results may indicate that predator–prey population dynamics have resulted in demographic expansions occurring in a relatively close time frame, particularly in the highly productive marine region of the Gulf of California. Furthermore, recent local and regional weather and oceanographic anomalies have been shown to affect both the marine productivity and the breeding distribution and success of several seabird species (Humphries et al. 2015, Velarde et al. 2015b). Studies of the demographic responses of seabirds to similar events in the past can help us have an early understanding of possible effects of such anomalies on species that may be particularly vulnerable. This information has important implications not only in the study of the genetics and evolutionary history of 2 very different species (both taxonomically and ecologically), but also for present and future conservation and management decisions regarding these seabirds and their prey species, as well as other predator–prey systems.
ACKNOWLEDGMENTS
We acknowledge support from the Secretaría de Marina-Armada de México for marine transportation. We thank Comisión Nacional de Áreas Naturales Protegidas, Baja California and Sonora regional offices, for logistic support; and T. Weiss for kindly providing the spreadsheet tool for the application of the mismatch analysis. J. Arce-Smith, R. Olachea, J. Aguilar, and T. Morales provided field support and helped with sample collecting. We thank three anonymous reviewers for valuable comments on earlier versions of the manuscript.
Funding statement: UCMEXUS-CONACYT provided support to E. Velarde and A. Aguilar through the Collaborative Research Grants Program for 2012, and a fellowship to E.A.R. (grant FE-10-53). Support was also provided to E.V. by the joint fund by Fondo Mexicano para la Conservación de la Naturaleza–Lindblad Expeditions/National Geographic–Packard Foundation.
Ethics statement: This study complies with Mexican regulations regarding the ethical treatment of research subjects. The Dirección General de Vida Silvestre, Subsecretaría de Gestión para la Protección Ambiental, Secretaría de Medio Ambiente y Recursos Naturales, granted the corresponding research and sample collection permits (SGPA/DGVS/02313/10, SGPA/DGVS/02017/11, SGPA/DGVS/01410/12, SGPA/DGVS/01522/12, and SGPA/DGVS/03382/12) and exportation permit (28758). The USFWS granted the importation permit (MB40037A) to A.A. Secretaría de Gobernación granted permits to stay and carry out research activities in Isla Rasa.
Author contributions: E.A.R., A.A., and E.V. conceived the study. E.A.R. and A.A. conducted and supervised the research. E.A.R. and E.V. collected the samples and data. E.A.R. performed the experiments, analyzed the data, and wrote the first draft of the manuscript. All authors proofread and edited the subsequent and final versions. A.A. and E.V. contributed materials, resources, and funding.
