Population genetic structure can aide in developing conservation and management strategies by characterizing populations on local and regional scales. The Elegant Tern (Thalasseus elegans) has a restricted breeding range, with a majority of its nesting population historically found on Isla Rasa, Gulf of California, Mexico. Since the late 1950s, increased frequency of low marine productivity due to warm oceanographic anomalies, such as El Niño events, have caused northward expansion of Elegant Tern breeding colonies into southern California, USA. To test the hypothesis that high gene flow occurs between Gulf of California and southern California Elegant Tern breeding colonies, restriction-site associated DNA sequencing was used to analyze 5,510 single nucleotide polymorphisms from 69 Elegant Terns sampled across four known breeding sites: Isla Rasa, Mexico (n = 30), San Diego National Wildlife Refuge, USA (n = 17), Bolsa Chica Ecological Reserve, USA (n = 11), and Port of Los Angeles, USA (n = 11). Analyses revealed little population subdivision, with non-significant genetic differentiation (FST) among sites and no geographic association of individuals, but there was subtle clustering of individuals by breeding site. These results suggest a strong degree of gene flow among the Gulf of California and southern California nesting colonies and indicate that Elegant Terns have a fluid breeding distribution and move readily among nesting sites.
Changes in species distributional patterns can provide valuable information on local, regional, and global responses to climate change and anthropogenic impacts (Pacifici et al. 2015; Beaugrand and Kirby 2018). Natural colonization of new habitats that lead to range expansions are useful instances to examine potential factors and consequence of species' range changes (Lenoir and Svenning 2015). Range expansions can lead to novel evolutionary and ecological adaptations to newly encountered environments (Rehm et al. 2015), and understanding of the underlying patterns and processes in these instances will be useful in projecting future conservation scenarios. From a conservation perspective, understanding the source/sink dynamics of newly colonized versus ‘core’ habitats can be informative in predicting the future persistence of all populations for a species, especially given current and predicted impacts of global climate change.
As top predators, seabirds are important links in food webs and are considered good environmental indicators of the health of marine ecosystems and abundance of prey species (Einoder 2009; Horn and Whitcombe 2015; Paleczny et al. 2015; Velarde et al. 2015). Fitness of seabirds is directly tied to populations of fish prey species (Velarde et al. 2015; Gomez-Laich et al. 2015; Cecere et al. 2015). Throughout the annual cycle, seabirds exploit a heterogeneous landscape (Cecere et al. 2015) and must take advantage of pulses in food availability to successfully breed, molt, and migrate. Generally, seabirds are colonial, long-lived, and exhibit pair-bonding. Most migrate after the nesting season and demonstrate high fidelity to nesting sites (Levin and Parker 2012; Clucas et al. 2016), but see Coulson (2016) who suggests a range of fidelity exists within seabird populations to counter seasonal variability in resources.
Of all bird groups, seabirds are the most imperiled, with 30% of Charadriiformes threatened with extinction (Paleczyny et al. 2015; Mancini et al. 2016). In addition to habitat loss and degradation, seabird populations are threatened by invasive species (e.g., rats), egging, diverse types of fisheries (i.e. as bycatch), and plastics pollution. These environmental challenges cause changes in breeding phenology (Wanless et al. 2009), and loss of suitable breeding and foraging habitats. In response to environmental challenges, some seabirds can change their patterns of migration and dispersal by expanding their foraging range (McLeay et al. 2010; Cecere et al. 2015), shifting breeding phenology, colonizing new breeding locations, or a combination thereof (Monaghan 1996; Velarde et al. 2015; Dayton et al. 2017).
The Elegant Tern (Thalasseus elegans) is included in the Mexico category of Special Protection for species that have a very restricted breeding distribution and are under risk of decline due to decades of intensive egg harvesting on Isla Rasa, Mexico (SEMARNAT 2010). In the United States, the tern is a species protected under the Migratory Bird Treaty Act of 1918, and it is considered Near Threatened by the IUCN (IUCN 2020). Elegant Terns have historically nested mainly in the Gulf of California islands, with their most important nesting site in recent years being Isla Rasa (Velarde et al. 2015), as attested by the popular L.W. Walker's (1951) account. Following the breeding season, a portion of the population moves northward, some as far north as British Columbia, Canada. Generally, the majority of the population can be found from San Francisco Bay, USA south to Baja California, the Gulf of California, and within Mexico from Sonora south to Nayarit, spending several weeks along the trajectory of their post-breeding migration (during Northern Hemisphere Summer). In November, a portion of the population then moves southward, with individuals spreading out to winter along the coast of southern Mexico south to northern Chile. The spring, northward migration begins in late February-early March (during Northern Hemisphere late Winter and early Spring), when terns begin to return to breeding grounds in Mexico and southern California (Bent 1986; Harrison 1983; IUCN 2020). Therefore, the species range spans from approximately latitude 40°N to 22°S, or about 9,000 km along the eastern Pacific coast. In this sense, all the previously reported and present nesting sites fall within the annual distribution range of the species; therefore, new colonies are not a range extension per se, but a breeding range extension.
Over the last two decades (1991-2014), the majority (∼78%) of nesting Elegant Terns have been found on Isla Rasa, Mexico (Velarde et al. 2015), a known breeding sanctuary and migratory stopover for many seabirds. During this interval, numbers have fluctuated dramatically on Isla Rasa (0-150,000 pairs). Since the late 1950s, the Elegant Tern has expanded its breeding range northward into southern California: first in 1959 at the San Diego National Wildlife Refuge (SDB) (Gallup and Bailey 1960), in 1987 at the Bolsa Chica Ecological Reserve (BCER) (Collins et al. 1991; Collins 2006), and then in 1998 at the Port of Los Angeles (POLA) (Burness et al. 1999). These four colonies constitute the most important nesting colonies (based on nest numbers) of the species.
The establishment of new breeding colonies in Southern California may be associated with El Niño-Southern Oscillation (ENSO) events that cause food shortages of pelagic fish e.g., Pacific sardine (Sardinops sagax) and northern anchovy (Engraulis mordax) (Mellink 2003; Velarde et al. 2015). This shortage, coupled with overfishing of its main food source (small pelagic fish, especially Pacific sardine), was documented by Velarde et al. (2015) to influence reproductive success of the tern. McLeay et al. (2010) suggested that the restricted range of crested terns (such as Elegant Tern) during the breeding season may make them more sensitive to competition with fisheries that operate within the terns' foraging ranges, influencing terns to disperse away from the area of restriction. An additional possibility is that when a colony reaches a certain size (number of nesting pairs), increased competition for food among nesting individuals may produce the Ashmole Halo effect (Ashmole 1971) and influence a tendency among nesting individuals to search for alternative nesting sites that may cause spillover to alternative nesting areas. Coupled with the assumed presence of site fidelity, this dispersion of nesting individuals suggests segregation can exist between populations (Kirven 1969).
There is also evidence that wintering distributions and migration routes may play an important role in the structuring of breeding populations for many seabird species (Szczys et al. 2012; Friesen 2015; Dayton et al. 2017; Szczys et al. 2017). Friesen et al. (2007) proposed that isolation during the non-breeding season is crucial to the development of genetic differentiation in seabirds, therefore one would expect little population structure among breeding colonies of Elegant Terns, despite the fact that geographically discrete breeding colonies do exist and that a northward expansion of Elegant Tern breeding colonies appears to be underway.
We investigated population subdivision in the Elegant Tern across the four most important breeding colonies to determine the level of connectivity between nesting populations in southern California, USA, and Isla Rasa, Mexico. We expected that Elegant Terns lack population subdivision because of the high level of connectivity produced by the spillover effect, the tern's dispersal capabilities, and the absence of discrete non-breeding areas. An alternative hypothesis is that subtle population subdivision is present between Elegant Tern colonies due to the possibility of site fidelity during non-El Niño years. Determining the plasticity or rigidity in Elegant Tern population structure is crucial to establish long-term population trajectories for this species, given the increasing variability of oceanic conditions that is negatively impacting local and regional seabird populations (e.g., Wiley et al. 2013, Velarde et al. 2015, Gagne et al. 2018).
Table 1.
Location and number of individual blood samples taken for genetic analysis among four Elegant Tern (Thalasseus elegans) breeding colonies in the Gulf of California, Mexico (RASA) and southern California, USA (SDB, BCER, POLA).
Methods
Sample Collection and Processing
We collected blood samples from 69 Elegant Terns (from a combination of adults and chicks) at the four main breeding colonies in 2013 and 2014 (Table 1; Fig. 1). On Isla Rasa, adult terns were captured using mist-nets during the twilight and early evening hours. Nets were placed away from the colony area to avoid disturbance, and along a site where we had noticed the formation of a flight corridor of terns approaching the colony during the daytime hours. After carefully removing the terns from the mist-net, a hood was placed over their head to calm the birds. The individuals were taken to the periphery of the mist-netting area to take blood samples. At the southern California colonies, blood sample collection was spread over a course of three weeks in July, dependent upon stage of reproduction within the colony. Chicks were hand captured within the breeding colony, and blood samples were drawn from the brachial vein. Chicks were released within 10-15 minutes of capture back to the chick creche. Exact age of chicks sampled was unknown, but chicks were old enough to be part of a formed creche (greater than 15 days).
Blood samples were drawn from the brachial vein and stored in Queen's lysis buffer and transported to California State University Los Angeles for DNA extraction. The DNA was isolated using standard Qiagen DNeasy Blood and Tissue Kit (Qiagen Inc.) following the manufacturer's protocols for blood samples. DNA was quantified using the Qubit Broad Range reagents (Life Technologies). From each blood sample, 400 ng of DNA was shipped to Global Biologics (Columbia, Missouri, USA) for Double-digest Restriction Associated DNA sequencing (ddRADseq) to isolate single nucleotide polymoprhisms (SNPs) by creating a genome-wide sequencing library of loci. This sequencing method was carried out using the protocol of Peterson et al. (2012). Prior to library preparation, DNA quality was checked on an Agilent Bioanalyzer. The DNA samples were digested with two endonucleases, EcoR1 and MspI. The fragments generated were adaptor ligated, and fragments between 150-400 bp were selected for sequencing. Samples were pooled and sequenced on an Illumina NextSeq500 with PE100 sequencing. Unique five-base barcodes were used to allowing pooling of sample on a single lane of sequencing.
Bioinformatic Analyses
The STACKS v.1.29 pipeline (Catchen et al. 2011) was used to demultiplex the raw fastq files and identify loci through de novo assembly and define SNPs across individuals. This pipeline compiles loci in three steps: (1) ustacks identifies unique loci based on polymorphism and infers alleles and haplotypes to compile loci for each individual; (2) cstacks generates a catalog of loci based on ustacks output and merges homologous loci; (3) sstacks searches for haplotype combinations present in those catalogs to match loci to an individual. In ustacks, we used a minimum of three reads to form a stack (-m) and allowed for three mismatches (-M) between stacks. We also allowed for three mismatches in the formation of stacks (-n) in cstacks.
Once the SNPs were categorized, the populations program in the STACKS pipeline was used to analyze and compute genetic statistics and generate various output formats such as VCF and genepop. The final dataset kept individuals with SNPs present at a minimum of 80% (-r option) at all 4 breeding colonies (-p option) with a minimum coverage of 10 reads (-m option). This was done to ensure a large proportion of called SNPs were shared across sample sites, and that adequate sequence coverage was obtained for all genotypes. Bootstrap resampling of 1,000 permutations was set to measure accuracy or significance to sample estimates.
PLINK (Purcell et al. 2007) was used to recode and extract SNPs based on a set of parameters that aim to increase the percentage of the genotype data shared between individuals. The parameters include removing individuals with more than 20% missing genotypes (-- mind 0.2), then removing loci with a minor allele frequency < 0.05 (--maf 0.05) and genotyping rate < 20% (--geno 0.2). Lastly, we filtered out those loci not in Hardy-Weinberg equilibrium and linkage disequilibrium.
Genetic Variation
Observed heterozygosity (Ho), expected heterozygosity (He), inbreeding coefficient (Fis), and allelic diversity were estimated using the R package diveRsity v. 1.9.89 (Keenan et al. 2013) for each breeding colony. A one-factor ANOVA was used to test for any difference between means of heterozygosity and allelic diversity among breeding colonies.
Population Subdivision
Pairwise estimates of population differentiation (FST) (Weir and Cockerham 1984) were also calculated using diveRsity (Keenan et al. 2013) and included bootstrapping with 1,000 permutations to estimate 95% confidence intervals. Population subdivision was further analyzed using a Principal Component Analysis (PCA) and Discriminant Analysis of Principal Components (DAPC), and assigning ancestry coefficients with the program sNMF. The R package adegenet (Jombart and Ahmed 2011) was used for both PCA and DAPC analyses.
PCA removes variables based on maximum variance and finds covariance between variables, which are plotted to spatially represent those relationships. These PCA data and a priori group membership were then used to conduct a DAPC analysis. The purpose of this analysis was to determine clusters in genetic data to create discriminant functions that will further classify individuals into groups based on overall variance among breeding colonies. We employed the cross-validation approach to determine the appropriate number of principal component axes to retain for the DAPC.
Lastly, the program sNMF (Frichot et al. 2014) was used to determine the number of groups or clusters and assign ancestry coefficients (K) by using sparse non-negative matrix factorization algorithm (sNMF). The input values used in sNMF were K = 1-6. A cross-validation technique was used to validate cross-entropy criterion where lower values indicate valid predictive results or ancestry estimation. The cross-entropy values are expected to decrease exponentially if distinct populations are present. This method is similar to ADMIXTURE and STRUCTURE; however, it is more computationally efficient for large datasets (Frichot et al. 2014).
Results
Sequencing Data
The average number of reads per individual recorded were: 1,864,635 for RASA; 134,758 for SDB; 759,525 for BCER; and 1,146,631 for POLA; with the average proportion of reads retained per individual = 0.25 (Table 2). After the program populations in PLINK was run, the total SNPs retained was 5,521. Raw sequence reads have been deposited in the NCBI short read archive (PRJNA604787).
Genetic Diversity
The number of polymorphic loci ranged from 5153 to 5510 (Table 3). The range of average overall observed heterozygosity ranged from 0.279 to 0.296 (Table 3) among all breeding colonies. The coefficient of inbreeding calculated among all breeding colonies varied slightly across sites, but the values were close to zero (Table 3). Mean allelic richness ranged from 1.933 to 1.936 (Table 3) among all breeding colonies. Based on ANOVA, no significant difference was found in means for heterozygosity (F = 2.3125; P = 0.128) or allelic richness across breeding colonies (F = 0.3572; P = 0.550).
Table 2.
Summary of DNA sequencing reads for Elegant Terns (Thalasseus elegans) from four breeding colonies in the Gulf of California, Mexico (RASA) and southern California, USA (SDB, BCER, POLA).
Population Subdivision
Pairwise FST among the four breeding colonies is summarized in Table 4 and was used to explain the level of population subdivision. Our analysis showed consistent FST values close to or equal to zero for each pairwise comparison of colonies. Results from the PCA, DAPC, and sNMF revealed no distinct population clusters among breeding colonies. The PCA showed strong overlap among the four breeding colonies (Fig. 2). The optimal number of principal components retained from cross-validation were 25 for the DAPC. The DAPC analysis also showed overlap of individuals among the four breeding colonies (Fig. 3). Cross-entropy ancestry coefficients (K) obtained using the sNMF program fluctuated rather than decreasing exponentially, the latter would be expected with populations that are subdivided.
Table 3.
Descriptive statistics and averages at each depth of coverage for population genetics of four Elegant Tern (Thalasseus elegans) breeding colonies in the Gulf of California, Mexico (RASA) and southern California, USA (SDB, BCER, POLA). Measurements of variation include number of polymorphic loci, observed and expected heterozygosity, inbreeding coefficient and allelic richness (standard error in parentheses).
Discussion
Our analyses supported the hypothesis of low to no population genetic subdivision and high levels of connectivity among the four breeding colonies (RASA, SDB, BCER, POLA) of the Elegant Tern, contrary to the stereotype of strong philopatry for seabirds. The SNP dataset showed no evidence for distinct Elegant Tern breeding colonies based on several analyses. Pairwise FST values show low levels of genetic differentiation, suggesting no barriers to gene flow among the colonies. The four breeding colonies overlapped in both ordination analyses (PCA and DAPC). The cross-entropy results computed in the sNMF program fluctuated with an increasing number of clusters (K), further supporting patterns of high gene flow among the contemporary breeding colonies of the Elegant Tern.
While we suggest the elevated level of gene flow is attributed to high connectivity between Mexico and southern California populations, the lack of genetic differentiation observed here could also be a result of the relatively recent establishment of breeding colonies in southern California, and the fact that distinct non-breeding areas are absent for this species. As mentioned above, the first recorded instance of Elegant Terns breeding in San Diego (USA) was in 1959, with subsequent establishment of breeding colonies farther to the north at BCER and POLA. These colonies have grown considerably since 1959, but their numbers have varied widely between years (see Table S1 in Velarde et al. 2015). The population fluctuations indicate that the growth of these colonies is due to both new recruits from these nesting sites as well as to immigration of individuals from the fast-growing Isla Rasa nesting population. Also, the southern California colonies may be too new to reveal population subdivision.
Table 4.
Pairwise FST for three depth of coverage datasets among four Elegant Tern (Thalasseus elegans) breeding colonies in the Gulf of California, Mexico (RASA) and southern California, USA (SDB, BCER, POLA).
Given that the establishment of these colonies is recent, and documented population movements from the southern California nesting areas to Isla Rasa colonies during non-El Niño years (E. Velarde, unpubl. data), our genetic results suggest that ongoing gene flow (panmixia) is occurring among contemporary Elegant Tern populations. Between 2001-2007, 11 terns were captured on Isla Rasa that possessed bands. Of these 11, three were originally captured and banded on Isla Rasa (∼27%) (E. Velarde, unpubl. data). The other eight individuals were originally banded in southern California (C. Collins, pers. comm.). The band-resight information also suggests that, at least in certain years, there is movement of individuals from the southern California colonies into the Isla Rasa colony, resulting in a very dynamic flow of breeding terns between colonies, possibly in response to food availability conditions.
One explanation for the maintenance of gene flow between Elegant Tern populations stems from their ability to respond to warm ocean anomalies by establishing new colonies or adding to existing ones (Velarde et. al. 2015). The connectivity between breeding locations for the Elegant Tern may confer a reproductive advantage in seasons where oceanic anomalies are influencing prey resource distribution. Tern populations participate in group adherence, which causes entire colonies to disperse when habitat quality and prey abundance decreases (Palestis 2014). This response has been observed in Elegant Tern populations, where northward expansion of breeding colonies was associated with warm ocean anomalies and fishery decline (Velarde et al. 2015).
Also, the fact that all nesting colonies occur within the annual range of the species indicates that there has not actually been a species distribution expansion but, rather, that the new colonies have been initiated within the post-breeding distribution of the species. This may have been initiated as a result of Isla Rasa individuals failing to nest on this island and starting a post-breeding migration north, where they found adequate food abundance and resulted in a second nesting attempt. Nesting in the southern California colonies has been recorded to start a few weeks after the Isla Rasa breeding initiation, which also indicates that the optimal environmental and food conditions for Elegant Tern breeding along the California Current are possibly a few weeks delayed from those of the Gulf of California. So, failed breeders leaving the Gulf of California and migrating north may well still arrive in time to attempt re-nesting on the southern California coast. This may be another explanation for the panmixia found in this species.
Another possible explanation for the lack of population structure in Elegant Terns is that distinct non-breeding areas do not exist for this species. In addition to spillover during breeding events, mixture during the non-breeding season could promote gene-flow as suggested by Friesen et al. (2007). The lack of discrete non-breeding sites in the Arctic Tern (Sterna vittata) was proposed to contribute to a lack of genetic differentiation among subspecies (Connan et al. 2015). Similarly, Faria et al. (2010) found low levels of population structure in the South American Tern (Sterna hirundinacea), which may be attributed to mixing during the non-breeding season. Faria et al. (2010) also suggest that the seasonal/annual movements of food resources such as anchovies (E. anchoita) may contribute to the movement of adult terns and therefore promote gene flow among colonies.
There have been variable results in finding population subdivision in terns. Population subdivision has been observed in other North American species of terns, such as Caspian Tern (Hydroprogne caspia), Common Tern (Sterna hirundo) and Least Tern (Sternula antillarum), and Eurasian species like the Eurasian Black Tern (Chlidonias niger niger) and the Whiskered Tern (Chilodias hybrida). Whittier et al. (2006) found genetic differentiation among the three subspecies of Least Terns. Hierarchical population structure, with asymmetric gene flow was found in the eastern range of Common Terns (Szczys et al. 2012). A subsequent genetic analysis of Least Terns (Draheim et al. 2010) also detected genetic differences among subspecies, but little genetic subdivision was observed within subspecies. Similarly, significant genetic differentiation was reported across the breeding range of the Caspian Tern (Boutilier et al. 2014). Szczys et al. (2017) found significant genetic differentiation among breeding sites and limited contact at nuptial staging areas in Black Terns, and Dayton et al. (2017) found evidence for two distinct subpopulations of European Whiskered Terns. The geographic scale at which these previous studies uncovered genetic differentiation is much larger than the scope of our study and points towards larger spatial and temporal isolation in the generation of genetic differences among North American and Eurasian tern populations.
Instances of low or no population differentiation have also been observed in other seabird species. Recent studies on petrels (Quillfeldt et al. 2017), penguins (Clucas et al. 2016; Cristofari et al. 2016; Gorman et al. 2017), gulls (Yannic et al. 2016), and murres (Tigano et al. 2017) found little evidence for genetic structure across the geographic range of those species. This lack of structure is particularly interesting given that many seabirds are thought to be highly philopatric and therefore expected to show genetic structure, even over small geographic distances (Friesen et al. 2007). Nevertheless, site fidelity is highly variable in seabirds as reviewed by Coulson (2016). Perhaps having variability in site fidelity, or not having site fidelity at all, is an adaptive feature in seabirds that allows for the ability to shift nesting locations to mitigate impacts to reproductive success (Crawford 2009). Reproductive success of the Elegant Tern has been associated with habitat quality and prey abundance (Velarde et al. 2015).
Connectivity among seabird breeding colonies has been hypothesized to buffer against genetic drift (Ramirez et al. 2013) and provide resilience to human disturbances (Geary et al. 2017). The presence of high connectivity between breeding locations provides evidence that the Elegant Tern population seems to be moderately buffered from warming events and would be able to respond to future climate variation if there are suitable nesting locations and an adequate prey base within reach of the colonies. Breeding colonies are currently found in stable reserves (RASA, SDB, and BCER) or in locations off limits to development (POLA). What remains uncertain is the variation in prey base for the Elegant Tern and other seabird populations in the Gulf of California and the southern California Bight system and impacts to these systems in the face of climate change (and more so if it is combined with unsustainable fisheries) (Cheung et al. 2009; Velarde et al. 2015). Further analyses of movement and landscape patterns would be beneficial in developing an understanding of habitat characteristics that lead to reproductive success that can then be used in conservation and management plans to address changes in nesting distributions for mobile species similar to the Elegant Tern.
Acknowledgments
We would like to thank Drs. K. Fisher, E. Wood and two anonymous reviewers for comments on a previous version of this manuscript. This material is based upon work supported by the California State University Council on Ocean Affairs, Science & Technology (COAST) under Award No. (COAST-RR-2014-002) to AA and MH. GSP was supported by a CSU-LSAMP fellowship (NSF Grant # HRD-1463889). Field work for collection of samples in Isla Rasa was supported by the joint fund provided by Fondo Mexicano para la Conservación de la Naturaleza/Lindblad Expeditions/Packard Foundation grant number M-FG-E-VA16-12-05-II and M-FG-E-VA16-12-05-III, and by the UCMEXUS program to EV. Collections were made possible by permits and permissions from the USFWS San Diego Bay National Wildlife Refuge Bird Banding Station Permit 2278 to KG, BBL Bird Banding Permit 20047 to RP, and SGPA/ DGVS/01954/13 and SGPA/DGVS/02131/14 from SEMARNAT, as well as SATI/PC/017/13 and SRPAP/ PC/016/14 from SEGOB to EV. Banding permits to EV were kindly provided by Dirección General de Vida Silvestre and USFWS: Master Bander Permit number 21380. All methods meet the ethical guidelines for the use of birds in Mexico. All applicable ethical guidelines for the use of birds in research have been followed, including those presented in the Ornithological Council's “Guidelines for Use of Wild Birds in Research” (Fair et al. 2010).