Geolocator tagging links distributions in the non-breeding season to population genetic structure in a sentinel North Pacific seabird.

We tested the hypothesis that segregation in wintering areas is associated with population differentiation in a sentinel North Pacific seabird, the rhinoceros auklet (Cerorhinca monocerata). We collected tissue samples for genetic analyses on five breeding colonies in the western Pacific Ocean (Japan) and on 13 colonies in the eastern Pacific Ocean (California to Alaska), and deployed light-level geolocator tags on 12 eastern Pacific colonies to delineate wintering areas. Geolocator tags were deployed previously on one colony in Japan. There was strong genetic differentiation between populations in the eastern vs. western Pacific Ocean, likely due to two factors. First, glaciation over the North Pacific in the late Pleistocene might have forced a southward range shift that historically isolated the eastern and western populations. And second, deep-ocean habitat along the northern continental shelf appears to act as a barrier to movement; abundant on both sides of the North Pacific, the rhinoceros auklet is virtually absent as a breeder in the Aleutian Islands and Bering Sea, and no tagged birds crossed the North Pacific in the non-breeding season. While genetic differentiation was strongest between the eastern vs. western Pacific, there was also extensive differentiation within both regional groups. In pairwise comparisons among the eastern Pacific colonies, the standardized measure of genetic differentiation (FꞌST) was negatively correlated with the extent of spatial overlap in wintering areas. That result supports the hypothesis that segregation in the non-breeding season is linked to genetic structure. Philopatry and a neritic foraging habit probably also contribute to the structuring. Widely distributed, vulnerable to anthropogenic stressors, and exhibiting extensive genetic structure, the rhinoceros auklet is fully indicative of the scope of the conservation challenges posed by seabirds.

Introduction

Barriers to dispersal are generally thought to be less rigid and absolute in marine environments than in terrestrial environments, implying that genetic connections in marine organisms will span longer distances, and population genetic structuring will be less prevalent [1-3]. Many marine organisms also exist at very high densities, which reduces the potential for genetic differentiation to arise from drift [4]. Nonetheless, a wide variety of effectively unbridgeable barriers to dispersal and gene flow, many of them cryptic, are present in the ocean. These barriers include spatial gradients in salinity [5] and temperature [6], eddies and gyres [7], fronts [8], ocean currents [9], irregular coastlines [10], and weak connections between ocean basins [11]. Consequently, population genetic structure exists in a broad range of marine taxa, including phytoplankton [12], vascular plants [13], crustacean zooplankton [14], molluscs [15], cephalopods [16], and vertebrates including marine fish [17], reptiles [18], and mammals [19]. Interpreting the complex patterns of genetic structuring displayed by marine organisms is a task that continues to challenge conservation researchers and managers [4, 20], suggesting the utility of targeted, hypothesis-driven approaches [21].

Among marine vertebrates, seabirds possess a suite of ecological and behavioural traits that are particularly relevant to investigations of population differentiation [22, 23]. A taxonomically diverse group usually considered to include ~360 species in six Orders, seabirds inhabit all of the world’s oceans and share a core suite of “slow” life-history traits: long pre-breeding periods, small clutch sizes, and high adult survival rates [24, 25]. Among the traits common to seabirds that are especially germane to population genetics are extreme mobility, philopatry, and high-density breeding in discrete groups (colonies). Most seabirds are strong fliers capable of rapidly traversing long stretches of ocean. However, the variation in dispersal capability across the group is considerable, with the flightless penguins (all species) and cormorants (one species) at one extreme, and the albatrosses, wide-ranging ocean wanderers [26], at the other. Such high capacity for dispersal violates the assumptions of many models of population divergence and speciation [22], and may render some of the smaller-scale oceanographic barriers that limit gene flow in other marine taxa irrelevant to seabirds. Conversely, philopatry restricts gene flow [27], thereby promoting population structure for species that return annually to the same breeding colonies. From a logistical perspective–spatial and temporal predictability, ease of access, sample sizes–the habit of colonial breeding generally makes seabirds good subjects for study [23].

Friesen, Burg, & McCoy [22] reviewed the literature on the extent and causes of population genetic structuring in seabirds based on mitochondrial DNA. They found that among all species studied, the extent of structuring varied from virtual panmixia to reciprocal monophyly, but structure was present in most. More specifically, genetic or genetic-plus-phylogeographic structure was found in all 12 species with breeding distributions fragmented by rigid physical barriers (land and/or ice), and in 26 of 37 species (70%) with continuous distributions. From among a suite of seven factors proposed to explain structuring in the latter group (geographic distance, the pattern of colony dispersion, distribution outside the breeding season, foraging range around colonies, population bottlenecks, retained ancestral variation, and cryptic physical barriers), the presence or absence of population genetic structure appeared to be best explained by distributions outside the breeding season. Structure was associated with the behavioural tendency to occupy multiple population-specific non-breeding areas, and/or to reside year-round at or near breeding colonies.

Two more recent reviews have further assessed the factors underlying genetic structuring in seabirds based on mitochondrial and nuclear DNA. The first of these reviews [27] concluded that geographic barriers, natal philopatry, and the occupation of multiple oceanic regions (e.g., basins, current systems) were important factors restricting gene flow. The second review [28] examined seabirds of the Southern Ocean, a marine region relatively free of physical barriers to dispersal [29], and identified ocean currents (especially for penguins) and philopatry as primary factors promoting structure. Of note, both of the newer studies confirmed one of the major conclusions of the original [22], specifically, that spatial separation in the non-breeding season was associated with restricted gene flow. Surprisingly, given its prominence in the reviews, the hypothesis is rarely tested directly. In two investigations to date, one found no support in populations of three species of seabirds that exhibited spatio-temporal segregation at sea, but little genetic structure [30]; and the second found support for the hypothesis, but concluded that gene flow between the two populations being investigated was restricted not by spatial segregation in the non-breeding season per se, but by population-level differences in the timing of breeding linked to habitat specialization at that time [31].

Here, we test the hypothesis that segregation in wintering areas is associated with genetic differentiation in the rhinoceros auklet, a widely-distributed North Pacific seabird that exhibits extensive and complex population genetic structuring [32]. The migratory movements of the species in the eastern Pacific are poorly known, however, beyond a putative southward post-breeding trajectory, e.g., “the bulk of the population appears to winter off California” [33]; “winters mostly offshore and along coasts, in North America mainly from southern British Columbia (casually from SE Alaska) south to Baja California” [34]. Our approach was to expand on the study of Abbott et al. [32] by increasing sample sizes for microsatellites, individuals, and colonies, while concurrently deploying light-level geolocator (GLS) tags to track the migration of rhinoceros auklets from 12 colonies distributed more-or-less continuously along the Pacific Ocean coastline from California to the Alaska Peninsula. For all colony-pairs, we then correlated the standardized measure of genetic differentiation (FꞌST) with the extent of spatial overlap in wintering areas. We also tested for isolation by distance (IBD) among the eastern Pacific Ocean colonies, to evaluate whether average dispersal distance accounts for genetic structure, assuming all populations have equal mean dispersal distances. Finally, we collected blood samples on five colonies in Japan, including one site where geolocator tags were deployed previously [35]. Our sampling thus spanned the global range of the species, the only exception being the Russian Far East, which was found to support large numbers in recent surveys [36].

Materials and methods

Study species, study locations, and oceanographic setting

The rhinoceros auklet is an abundant, colonial, burrow-nesting seabird with a wide distribution across the temperate North Pacific Ocean. The species is misnamed, in that it does not belong to the auklet clade (Aethiini) within the family Alcidae. Rather, the genus Cerorhinca is either basal within the puffin clade (Fraterculini), or it forms a sister clade with Fratercula [37].

The main part of our study occurred on 13 rhinoceros auklet breeding colonies in the eastern North Pacific Ocean (Fig 1). These consisted of three colonies in Alaska (Chowiet, Middleton, and St. Lazaria islands), six colonies in British Columbia, which supports the majority of the North American breeding population (Lucy, Moore, Pine, Triangle, and Cleland islands, plus S’Gang Gwaay), two colonies in Washington State (Protection and Destruction islands), and two colonies in California (Southeast Farallon and Año Nuevo islands). The North Pacific Current flows east from Asia, bifurcating roughly off southern British Columbia into the southward-flowing California Current, an upwelling system, and the northward-flowing Alaska Current, a downwelling system [38]. These two systems exert the major oceanographic influences around the 13 study colonies. From a historical perspective, many of the colonies north of California would have been covered by the Cordilleran Ice Sheet at its maximum extent ca. 19–14.5 kya near the end of the Pleistocene [39, 40].

Fig 1

Locations of rhinoceros auklet breeding colonies along the Pacific Coast of North America where samples were collected for genetics analysis and light-level geolocator tags were deployed.

For the genetics we also included five western North Pacific colonies, all in Japan: three in the Japan Sea off Hokkaido (Todojima, Teuri, Matsumae-Kojima Islands), one along the Pacific coast of Hokkaido (Daikoku Island), and one in Mutsu Bay in Aomori (Taijima Island). Geolocator tags were deployed previously on rhinoceros auklets on Teuri Island [35]. The location and population size of all 18 colonies are provided in Table 1.

Table 1

Locations of rhinoceros auklet breeding colonies included in the genetics study, and the sizes of their breeding populations.

Breeding colony	Location	Population size (pairs)	Latitude (°N)	Longitude (°W)
Daikoku	Japan	78,000 [41]	42.954744	144.866057
Taijima	Japan	150 [41]	41.263284	140.345585
Matsumae-Kojima	Japan	40,000 [41]	41.360155	139.818113
Teuri	Japan	379,000 [41]	44.417646	141.312281
Todojima	Japan	35,000 [41]	45.366808	141.035141
Chowiet	Alaska	400 [42]	56.015513	-156.740272
Middleton	Alaska	10,000 (S. Hatch, unpublished data)	59.415254	-146.345472
St. Lazaria	Alaska	2,000 [43]	56.986502	-135.710838
Lucy	British Columbia	25,000 [44]	54.294418	-130.621907
S'Gang Gwaay	British Columbia	14,000 [44]	52.092634	-131.225633
Moore	British Columbia	40,000 [44]	52.678344	-129.418847
Triangle	British Columbia	42,000 [44]	50.851023	-129.066292
Pine	British Columbia	90,000 [44]	50.976062	-127.729909
Cleland	British Columbia	1,000 [44]	49.171516	-126.091075
Destruction	Washington	6,500 [45]	47.674599	-124.484817
Protection	Washington	36,000 [45]	48.126341	-122.930289
Southeast Farallon	California	4,500 [46]	37.695357	-123.000752
Año Nuevo	California	330 [46]	37.107584	-122.337026

Genetics methods

Research protocols employed in this study were approved by Simon Fraser University Animal Care Services (#974B-94), the Western and Northern Animal Care Committee of Environment and Climate Change Canada’s Canadian Wildlife Service (15MH01 and 16MH01), US Geological Survey Federal Bird Banding Permit (#09316, #20570), the Animal Ethics Committee of Hokkaido University, and the Aomoru prefecture (#3021). Access to field sites was provided by the Año Nuevo State Park, the California Department of Parks and Recreation, the United States Fish and Wildlife Service, the Washington Department of Fish and Wildlife, British Columbia Parks, the Archipelago Management Board of Gwaii Haanas National Park Reserve, and the Japanese Agency for Cultural Affairs.

Rhinoceros auklets involved in the genetics (2010–2016) and geolocator tagging (2014–2016) studies were captured using a variety of methods. Most were removed from breeding burrows, others blocked while departing the colony in knock-down or mist nets, others caught in purse nets or noose mats set in the entrances of burrows, and others captured by hand on the surface.

Blood samples (1.0 ml) and/or feather samples (5–10, from the breast) were collected from seven to 80 individuals on each of the 18 colonies, for a total of 704 individuals. All samples collected for the earlier study [32] were included, and supplemented with additional samples from Chowiet and St. Lazaria islands, and from 10 new sites. Samples were stored in ethanol and then at -20°C upon return to the lab. DNA was extracted from samples collected at four sites in Japan (Daikoku, Taijima, Matsumae-Kojima, Todojima) using the DNeasy® Blood and Tissue Kit (Qiagen), and at one of our California sites (Año Nuevo, where muscle tissue from nestlings found dead in burrows was used) using Macherey-Nagel DNA extraction kit. For all remaining samples, DNA was extracted using a modified Chelex protocol [47, 48].

A small set of samples (three to six) was genotyped to check for amplification and polymorphism with microsatellite loci from the genomes of rhinoceros auklet, crested auklet (Aethia cristatella), and whiskered auklet (Aethia pygmaea). Of the 31 microsatellite loci tested, 12 were polymorphic: CMms2, CMms3, CMms4, CMms9, CMms14, CMms22, CMms23, and CMms26 [49]; Pal11 and Pal26 [50]; Apy06 and Apy09 [51]. The eight CMms loci were the same ones used in the earlier study [32]. Due to inconsistent amplification, Pal26 and Apy09 were removed from the study and individuals were genotyped at the remaining ten microsatellite loci.

DNA was amplified in a 10 μL polymerase chain reaction (PCR) containing colorless GoTaq® Flexi (Promega) or Truin buffer, 0.2 mM dNTP, 2.5 mM MgCl2, 0.5 μM forward and 1 μM reverse primers, 0.05 μM M13 fluorescently labelled primer and 0.5 U GoTaq® Flexi or 1 U Truin Taq polymerase. The lone exception was CMms4, for which 3 mM MgCl2 was used. One percent formamide was added to CMms4, CMms9, and CMms14 reactions. A M13 sequence was added to the 5’ end of the forward primer to allow incorporation of a fluorescently-labelled M13 primer for visualization of the PCR product. The Apy09 reverse primer was pig-tailed (gtttctt) at the 5’ end to improve scoring [52]. Nine loci were amplified using the following thermal cycling profile: 120 s at 94°C, 45 s at 55°C, and 60 s at 72°C; seven cycles of 60 s at 94°C, 30 s at 55°C, and 45 s at 72°C; 31 cycles of 30 s at 94°C, 30 s at 57°C and 45 s at 72°C; and a final cycle at 72°C for 300 s. The tenth locus, Pal11, was amplified using the same program except the steps annealing Ta1 and Ta2 were increased to 60°C and 62°C respectively.

PCR products were run on a 6% polyacrylamide gel using a LI-COR 4300 DNA Analyzer (LI-COR Inc.). Alleles were scored via visual inspection with all genotypes confirmed by a second person. To maintain consistent scoring, three positive controls of known allele size were present on each load.

Genetic diversity analyses

Seven to 53 individuals remained from each breeding colony after excluding individuals missing three or more loci for a total of 424 individuals from 18 breeding colonies genotyped with ten microsatellites loci (Table 2). For the seven individuals from the Southeast Farallon colony, CMms3 failed to amplify. CMms26 had a high percentage of missing data (> 35%). Thus, analyses sensitive to missing data excluded the colony at Southeast Farallon and the CMms26 locus. GENEPOP v4.2 [53, 54] was used to check colonies and loci for linkage disequilibrium and deviations from Hardy-Weinberg equilibrium (HWE) using Markov chain parameters of 1000 iterations, 300 batches, and 2,000 dememorization steps. MICRO-CHECKER v2.2.3 [55] was used to detect scoring error due to stutter, null alleles, and drop out of large alleles.

Table 2

Rhinoceros auklet breeding colonies (populations) included in the genetics analyses, sample size (n), number of alleles (Na), expected heterozygosity (HE), observed heterozygosity (HO), private alleles (PA), and allelic richness (AR).

Breeding colony	n	Na	HE	HO	PA	AR
Daikoku	20	5.3	0.65	0.68	1	3.89
Taijima	10	4.8	0.65	0.77	0	3.97
Matsumae-Kojima	25	5.3	0.64	0.71	0	3.76
Teuri	21	5.4	0.63	0.66	1	3.61
Todojima	12	4.6	0.64	0.67	0	3.75
Chowiet	18	4.9	0.59	0.59	0	3.58
Middleton	50	6.3	0.64	0.60	1	3.71
St. Lazaria	22	5.9	0.65	0.65	1	3.87
Lucy	26	5.3	0.60	0.63	0	3.53
S'Gang Gwaay	53	5.5	0.62	0.61	1	3.56
Moore	18	4.8	0.62	0.62	0	3.55
Triangle	26	4.7	0.61	0.62	0	3.44
Pine	27	5.6	0.63	0.60	0	3.64
Cleland	30	5.8	0.66	0.64	1	3.85
Destruction	18	4.9	0.66	0.71	0	3.94
Protection	19	5.2	0.61	0.63	0	3.71
Southeast Farallon	7	3.3	-	-	-	-
Año Nuevo	22	5.5	0.66	0.68	1	3.83

For some measures, values are excluded for Southeast Farallon Island because of small sample size and data missing from locus CMms3.

Levels of population genetic diversity, observed and expected heterozygosities, private alleles, and number of alleles per locus were calculated in GenAlEx v6.5 [56, 57]. Because smaller sample sizes are expected to have fewer alleles, allelic richness was calculated using statistical rarefaction in HP-Rare v1.1 [58] which standardizes measurements to account for differences in sample size.

GenAlEx v6.5 [56, 57] was used to calculate both global and pairwise FST and FꞌST values to characterise population genetic structure with 999 permutations used to test significance. To address the multi-allelic nature of microsatellites, FꞌST provides a standardized value by dividing each FST with the maximum possible FST for the data [59]. For all tests with multiple comparisons, statistical levels of significance were adjusted using the modified False Discovery Rate (FDR) correction [60].

STRUCTURE v2.3.4 [61, 62], a non-spatial Bayesian method, was used to determine the number of genetic clusters present among populations. STRUCTURE captures the underlying population structure of the data without overestimating it. However, the program struggles to cluster individuals when genetic differentiation is low (FST ≤ 0.03) [63] and substructure might be present among populations. Therefore, we applied a hierarchical STRUCTURE approach (i.e., separate multistep runs with admixed individuals) that might detect genetic structure that is not apparent when all populations are run together [64].

STRUCTURE v2.3.4 was run using the recommended admixture ancestry model and correlated allele frequencies with sampling locations as locpriors. The locpriors option can be informative when population structure is weak [65]. For the 18 breeding colonies, ten independent runs were completed for each value of K from 1–6 with a burn-in of 100,000 and 120,000 Markov chain Monte Carlo (McMC) repetitions. The most appropriate number of clusters (K) was determined using several methods including lnPr(X|K) values, STRUCTURE HARVESTER v0.6.94 [66] and Bayes factor [62]. To further measure whether substructure was present within the populations, a hierarchical analysis for both the western (n = 5) and eastern (n = 13) breeding groups was completed using a burn-in of 60,000 and 70,000 McMC repetitions. K ranged from K = 1–4 for the western group and K = 1–5 for the eastern group. Five additional runs were completed at the optimal K to ensure convergence.

GenAlEx v6.5 [56, 57] was used to perform a principal coordinates analysis (PCoA) using standardized covariance from the FꞌST pairwise matrix. Three separate analyses were performed: all colonies, only western Pacific colonies, and only eastern Pacific colonies.

Geolocator deployments

We deployed light-level based Global Location Sensing (GLS) tags (Intigeo-C65, Migrate Technology Ltd., 1.0 g, range 4, mode 6) on breeding rhinoceros auklets on 13 eastern Pacific colonies in 2014 and 2015 (Fig 1). Most tags were deployed during the second half of the offspring provisioning period, as close to fledging as possible. Each tag was pre-attached to a Darvic leg band using self-amalgamating tape, epoxy, and an ultraviolet-resistant zip tie with a stainless-steel barb threaded through two custom holes in the band such that the whole assembly (< 2.0 g) could be quickly attached to the tarsus. Tags were groundtruthed (run at a known location in the colony) for ≥ 3 days prior to deployment. At deployment, the morphometrics of each bird were recorded, a numbered USGS metal band and a geolocator tag assembly were applied to the right and left tarsi, respectively, and the bird was released. Efforts were made to recapture birds and retrieve the tags beginning early in the following breeding season. Retrieved tags were groundtruthed again for ≥4 days. It was not possible to test for tag-induced changes to non-breeding behavior, and because birds can change burrows or nesting areas among years, failure to recapture a bird during this study does not necessarily indicate mortality.

Geolocation data analyses

We used the TwilightFree package to derive daily locations from the light level data [67]. This approach uses a Hidden Markov Model in which the hidden states are the daily geographic locations and the measured response is the observed pattern of light and dark over a 24 h period, and identifies the most likely tracks (based on the centre coordinates of a pre-defined grid) predicted to have generated the light record. The user specifies the spatial grid, minimum light threshold, zenith, and parameters related to the probability of light sensor shading and of movements between grid cells [67].

Based on the ground-truthing periods, a preliminary calibration was conducted on each light record to determine the minimum light threshold and zenith, which were used for all runs (zenith = 96°, threshold = 10 lux). We defined a grid that spanned from -180 to -100° longitude, and 10 to 65° latitude, and used a sea mask to restrict movements to the ocean. Shading likelihood parameters were selected to represent moderate amounts of shading (alpha = 5), and movement probabilities based on other seabirds (beta = 5) [67]. Each lux file was then processed wherein the model was initially fit with a grid of 4°×4° cell sizes to identify days with missing light data, which were then excluded from subsequent analyses. This reduced light record was then refit on a finer grid of 1°×1° cell sizes, and the track files derived from this second fit used in further analyses.

Home ranges in winter

The breeding seasons of rhinoceros auklets on colonies in the eastern Pacific Ocean span the period from April to August, generally starting and ending about one month earlier in the south than in the north [33]. In defining a non-breeding season for this study, we restricted tracks to the period from 1 November to 28/29 February. These dates roughly span the period from the end of fall migration (post-breeding dispersal from the colony) to the start of spring migration (pre-breeding return to the colony). These dates also avoid errors in geolocation due to the confounding of movements and light levels that occurs at the equinox.

We characterized the non-breeding distribution of each bird by its Utilization Distribution (UD), the bivariate function giving the probability density that an animal is found at a point in space according to its geographical coordinates [68]. The daily spatial coordinates were used to calculate the UD for each individual bird, as well as an aggregate from all the tags deployed in each colony. The two years of deployment (2014, 2015) were initially considered separately, and then were consolidated following the observation that tags from each colony had very high overlaps across years. We used function kernelUD from package adehabitatHR in R [69] to calculate the UDs, with the ad hoc method to derive the smoothing parameter and estimated over a grid of 1°×1° cell sizes. The function getverticeshr was used to extract the 50 and 90 percent home-range contours for each bird or colony aggregate. These contours were projected to a Lambert Equal‐Area Azimuthal projection, and the area of the polygons was calculated in km2.

Spatial overlap

The function kerneloverlaphr was used to calculate the UD overlap index (UDOI) between each pairwise combination of colony UDs [70]. The UDOI is 1.0 when both UDs are uniformly distributed and share 100% overlap, 0 when UDs have no overlap, and > 1.0 when the two UDs are non-uniformly distributed and share high overlap.

Correlations with Fꞌst

To test whether genetic distance between rhinoceros auklets from different breeding colonies was related to (i) the distance between breeding colonies, or (ii) the spatial overlap among colonies in winter, we correlated these two measures with the standardized measure of genetic differentiation (Fꞌst) values based on Wright’s fixation indices. We calculated Spearman’s correlation (r) between Fꞌst values and (i) distance between breeding colonies (IBD) or (ii) the UDOIs. Given the lack of independence between pairwise values, we calculated the statistical significance of r-values using a randomization procedure similar to a Mantel test. We built a sampling distribution of r through a permutation of holding Fꞌst values and randomly assigning it a pair from the available UDOI or IBD values. This was repeated 5000 times, and the proportion of times when the randomized r was larger than the observed r was used as an approximate P-value. Given the variability in the underlying data, we used an alpha-level of 0.10 to minimize the probability of a Type II statistical error.

Results

Over all colonies and loci, the number of alleles ranged from 2 to 16 with seven colonies having one private allele. Excluding the Southeast Farallon colony due to small sample size and missing data from locus CMms3, overall mean expected heterozygosity across all loci and samples was 0.62. Expected heterozygosity ranged from 0.59 (Chowiet) to 0.66 (Cleland, Destruction and Año Nuevo) with observed heterozygosity 0.59 (Chowiet) to 0.77 (Taijima). Allelic richness (corrected to a sample size of 10) was similar among colonies, ranging from 3.44–3.97 (Table 2). After FDR correction, two loci (CMms2, CMms22 at St. Lazaria) showed deviations from HWE (Table 2). There was no evidence for linkage disequilibrium between any of the loci. MICRO-CHECKER found no evidence of null alleles, large allele dropout or scoring error due to stutter.

FST statistics excluded the Southeast Farallon colony and CMms26 locus because of small sample sizes and missing data, respectively. Using 17 colonies and nine microsatellite loci, global FST was 0.039 (p < 0.001). Pairwise FST values ranged from 0.000 (Pine and Chowiet) to 0.112 (Lucy and Teuri), and FꞌST values ranged from 0.000 to 0.307 (Table 3). After FDR correction, 106 out of 136 tests were significant indicating a high level of genetic differentiation among all 17 remaining colonies.

Table 3

Pairwise FST values (below diagonal) and FꞌST values (above diagonal) for 17 rhinoceros auklet breeding colonies based on nine microsatellite loci.

		West					East
West		DAI	TAI	MAT	TEU	TOD	CH	MID	STL	LU	SGG	MO	TRI	PI	CL	DE	PR	AN
	DAI		0.044	0.026	0.117	0.046	0.264	0.254	0.252	0.301	0.288	0.214	0.254	0.258	0.196	0.166	0.213	0.174
	TAI	0.013		0.075	0.228	0.107	0.202	0.215	0.254	0.254	0.284	0.188	0.251	0.221	0.138	0.152	0.192	0.166
	MAT	0.008	0.022		0.122	0.034	0.205	0.196	0.212	0.246	0.245	0.175	0.166	0.237	0.168	0.158	0.188	0.174
	TEU	0.036	0.071	0.038		0.034	0.302	0.272	0.284	0.307	0.275	0.269	0.223	0.248	0.293	0.305	0.283	0.244
	TOD	0.014	0.032	0.010	0.011		0.234	0.196	0.218	0.274	0.250	0.191	0.175	0.217	0.191	0.166	0.203	0.200
East	CH	0.086	0.067	0.067	0.105	0.078		0.009	0.010	0.008	0.013	0.045	0.037	0.000	0.034	0.070	0.019	0.051
	MID	0.083	0.071	0.064	0.093	0.065	0.003		0.037	0.005	0.009	0.062	0.037	0.019	0.053	0.080	0.076	0.093
	STL	0.070	0.075	0.064	0.090	0.066	0.003	0.012		0.045	0.041	0.052	0.029	0.027	0.031	0.054	0.062	0.039
	LU	0.104	0.090	0.086	0.112	0.098	0.003	0.002	0.016		0.000	0.098	0.064	0.012	0.066	0.124	0.097	0.136
	SGG	0.098	0.098	0.084	0.098	0.087	0.005	0.003	0.014	0.000		0.077	0.036	0.005	0.075	0.104	0.061	0.094
	MO	0.070	0.062	0.057	0.093	0.064	0.016	0.022	0.017	0.037	0.028		0.064	0.062	0.026	0.103	0.061	0.052
	TRI	0.081	0.081	0.055	0.076	0.058	0.013	0.013	0.010	0.024	0.013	0.023		0.033	0.081	0.098	0.099	0.103
	PI	0.084	0.073	0.078	0.085	0.072	0.000	0.007	0.009	0.004	0.002	0.022	0.012		0.043	0.103	0.069	0.074
	CL	0.056	0.039	0.048	0.089	0.055	0.011	0.017	0.009	0.022	0.025	0.008	0.025	0.014		0.030	0.033	0.061
	DE	0.048	0.044	0.047	0.095	0.049	0.023	0.026	0.016	0.043	0.035	0.034	0.032	0.034	0.009		0.032	0.024
	PR	0.066	0.060	0.059	0.094	0.064	0.007	0.026	0.020	0.035	0.022	0.021	0.034	0.024	0.010	0.010		0.037
	AN	0.052	0.049	0.052	0.077	0.060	0.017	0.031	0.012	0.048	0.032	0.017	0.033	0.024	0.018	0.007	0.012

One population and one locus were removed because of small sample size and missing data (Southeast Farallon Island and CMms26). Bold values indicate statistical significant at p < 0.05; with underlined values at p < 0.01 following corrections for multiple tests. Acronyms are DAI = Daikoku, TAI = Taijima, MAT = Matsumae, TEU = Teuri, TOD = Todojima (all in the western Pacific); and CH = Chowiet, MID = Middleton, STL = St. Lazaria, LU = Lucy, SGG = S’Gang Gwaay, MO = Moore, TRI = Triangle, PI = Pine, CL = Cleland, DE = Destruction, PR = Protection, AN = Año Nuevo (all in the eastern Pacific).

The five western Pacific Ocean colonies were significantly differentiated at p < 0.001 from the 13 eastern Pacific Ocean colonies based on FST values. The only exception was that Taijima in the western Pacific was significantly different at p < 0.05 from Destruction and Año Nuevo islands in the eastern Pacific, possibly due to the small sample size at Taijima (n = 10).

Among the western Pacific colonies, FST significance varied between colonies. All five significant values included Teuri and Taijima, both of which were significantly different from three of the other four colonies (Table 3). Among the eastern Pacific breeding colonies, 41 of the 66 pairwise comparisons were significant (Table 3). Interestingly, the three Alaskan colonies (Chowiet, Middleton, St. Lazaria) accounted for most (16 of 25) of the non-significant values. Excluding those three sites in Alaska, 8 of the remaining 9 non-significant values were between nearby colonies, either in the central (Lucy-S’Gang Gwaay, Lucy-Pine, S’Gang Gwaay-Pine) or the southern (Cleland-Protection, Cleland-Destruction, Protection-Destruction, Protection-Año Nuevo, Destruction-Año Nuevo) regions of our study. The other non-significant pairwise comparison was between Cleland and Moore islands in British Columbia.

Visual inspection of the STRUCTURE and ΔK plots showed two distinct clusters, one including all five western Pacific Ocean colonies and one including all 13 eastern Pacific Ocean colonies. Most individuals had ancestry coefficient Q > 80% for one of the two clusters (Fig 2A). Comparison between the average lnPr(X|K) values at K = 2 (-10381) and K = 3 (-10352) indicated three genetically distinct populations. At K = 3, most individuals from the five western Pacific colonies had Q > 70% for the same cluster, but the eastern colonies split into two clusters: Q > 60% for one cluster consisting of Chowiet, Middleton, St. Lazaria, Lucy, S’Gang Gwaay, Triangle, Pine, Destruction, Protection, and Southeast Farallon; and Q > 60% for a second cluster consisting of Moore and Cleland, with 60% of individuals from the Año Nuevo colony showing Q > 50% for cluster two (Fig 2B).

Fig 2

Genetic groups for 18 North Pacific Ocean rhinoceros auklet colonies for ten microsatellite loci as inferred by STRUCTURE v2.3.

Histogram plot shows ancestry coefficient (Q) on the y-axis with individuals plotted on the x-axis. Breeding populations are listed from the western North Pacific (left) to eastern North Pacific (right). Genetic clusters at a) K = 2 and b) K = 3. (a) At K = 2 two genetic clusters occur between the western Pacific (Q > 70% yellow; Daikoku (DAI), Tajima (TAJ), Matsumae (MAT), Teuri (TEU), Todojima (TOD)) from the eastern Pacific (Q > 80% purple; Chowiet (CH), Middleton (MID), St. Lazaria (STL), Lucy (LU), S’Gang Gwaay (SGG), Triangle (TRI), Moore (MO), Pine (PI), Cleland (CL), Destruction (DE), Protection (PR), Southeast Farallon (SEF), Año Nuevo (AN)). (b) At K = 3 one cluster occurs for the western Pacific (Q > 70% yellow). For the eastern Pacific there are two clusters: Q > 60% for Chowiet, Middleton, St. Lazaria, Lucy, S’Gang Gwaay, Triangle, Pine, Destruction, Protection, Southeast Farallone (in purple); and Q > 60% for Moore and Cleland, with 60% of individuals from the Año Nuevo colony showing Q > 50% for cluster two (in light blue).

The hierarchical STRUCTURE analysis for the western Pacific colonies showed K = 1 (highest ln(Pr(X|K) = -2157). For the eastern Pacific colonies, the highest lnPr(X|K) value occurred at K = 2 (-8177). The eastern Pacific colonies split into two clusters at Q > 60%. The first cluster included Chowiet, Middleton, St. Lazaria, Lucy, S’Gang Gwaay, Triangle, Pine, Destruction, and Southeast Farallon, while the second cluster included Moore, Cleland, Protection and Año Nuevo (Fig 3).

Fig 3

Hierarchical analysis of 13 eastern Pacific rhinoceros auklet colonies for ten microsatellite loci as inferred by STRUCTURE v2.3.4.

Histogram plot shows ancestry coefficient (Q) on the y-axis with individuals plotted on the x-axis. Populations are listed as they move along the coastline from Alaska to California. Substructure was found for the 13 eastern breeding colonies, and supported by mean (ln Pr(X|K)) value and delta K (ΔK) at K = 2. Substructure includes two clusters at Q > 60%, one including Chowiet, Middleton, St. Lazaria, Lucy, S’Gang Gwaay, Triangle, Pine, Destruction, and Southeast Farallon (purple), and the second including Moore, Cleland, Protection, and Año Nuevo (light blue). No additional substructure was found among either the western or eastern breeding groups.

The PCoA with ten loci and 17 breeding colonies showed clear separation between the western Pacific vs. eastern Pacific colonies. The first two axes explained 55.7% and 17.5% of the variation (third axis 8.4%; Fig 4) and the clustering is concordant with the STRUCTURE results. When the western Pacific colonies were examined alone, the first two axes explained 77.5% and 15.4% of the variation (third axis 6.8%). Three clusters were evident: Teuri with Todojima; Matsumae-Kojima with Daikoku, and Taijima. When the eastern Pacific breeding colonies were examined alone, the first two axes explained 47.3% and 16.5% (third axis 13.2%) of the variation. The majority of colonies in Alaska and British Columbia clustered together (Chowiet, Middleton, Lucy, S’Gang Gwaay, and Pine) with Triangle forming a separate cluster and the two WA colonies (Destruction and Protection) forming a third cluster. The remaining colonies of St. Lazaria, Moore, Cleland, and Año Nuevo showed some degree of separation from all of the other colonies.

Fig 4

PCoA analysis conducted in GenAlEx v6.5 based on pairwise FꞌST values for the western and eastern North Pacific rhinoceros auklet breeding colonies.

Coordinates 1 and 2 explain 55.7% and 17.5% of the variation (not shown: coordinate 3 at 8.4%). The colours correspond to the three groups detected using STRUCTURE (Figs 2 and 3): western Pacific (yellow), the larger eastern Pacific group (purple) and smaller eastern Pacific group (blue).

Geolocator tagging results

Of the total of 370 geolocator tags deployed on rhinoceros auklets, 150 tags were retrieved (Table 4). Due to logistical issues we were unable to retrieve tags on Moore Island in 2016. Retrieval rates varied widely by colony and year ( = 40.3%, range = [0.0, 80.0]). We were able to obtain 142 lux files and 141 tracks following geolocation analyses of data from these devices, with some losses due to file corruptions or device failures.

Table 4

Sample sizes and retrieval rates of GLS devices deployed on rhinoceros auklets at breeding colonies along the Pacific Coast of North America, 2014–2015.

Colony	Year	Deployed	Retrieved	Retrieved (%)	Track files
Chowiet	2014	14	5	35.7	4
Middleton	2014	20	16	75.0	16
St. Lazaria	2014	17	8	47.1	8
Lucy	2014	25	12	48.0	12
Lucy	2015	30	12	40.0	12
Moore	2015	5	0	0.0	0
S'Gang Gwaay	2014	30	11	36.7	10
S'Gang Gwaay	2015	31	14	45.2	14
Triangle	2014	30	13	43.3	12
Triangle	2015	31	16	51.6	16
Pine	2014	30	4	13.3	4
Pine	2015	22	5	22.7	4
Cleland	2014	20	2	10.0	2
Protection	2015	23	9	39.1	7
Destruction	2015	7	2	28.6	2
Farallones	2014	20	9	45.0	9
Año Nuevo	2014	15	12	80.0	9
Total		370	150	40.3	141

Home range estimates from the geolocator tags varied widely by colony, but in general rhinoceros auklets remained in continental shelf waters during the winter months (Fig 5). The mean area of colony aggregates was 570,408 km2 for 50 percent UD, and 2,216,996 km2 for the 90 percent UDs (approximately 3.5x larger, a pattern that was consistent for individual-level UDs). For individual birds, the mean 50 percent UD had an area of 331,847 km2, and 1,112,287 km2 for the 90 percent UDs. In general, there was a latitudinal cline in home range areas, with markedly smaller home ranges for birds breeding at the two colonies in California, Southeast Farallon and Año Nuevo (Table 5, Fig 5).

Fig 5

Home ranges of rhinoceros auklets during the wintering period (Nov-Feb) from breeding colonies on the Pacific coast of North America, 2014 to 2016.

Light green dots indicate the locations of breeding colonies. Light gray marks are point locations on a 1°×1° grid, as derived from TwilightFree geolocation analyses of light recorders deployed on all birds at each colony. Red polygons indicate the 50% Utilization Distribution (‘core wintering area’), and dark red polygons indicate the 90% Utilization Distribution.

Table 5

Home-range area estimates (km2) during wintering months (Nov-Feb) of rhinoceros auklets breeding at 12 colonies on the Pacific Coast of North America, 2014 to 2016, as calculated from the Utilization Distribution (UD) derived from geolocation of light loggers.

Colony	Colony (50p)	Colony (90p)	Individual Mean (50p)	Individual SD (50p)	Individual Mean (90p)	Individual SD (90p)
Chowiet	1,415,726	5,120,609	866,352	875,273	2,875,308	2,643,150
Middleton	627,880	3,106,851	338,206	266,197	1,075,572	825,712
St. Lazaria	408,793	2,428,101	420,982	355,115	1,399,460	1,147,132
Lucy	650,667	2,311,713	502,806	486,398	1,610,695	1,550,292
S'Gang Gwaay	637,480	2,313,936	423,452	394,763	1,409,576	1,303,438
Triangle	837,116	2,687,246	473,992	378,839	1,496,879	1,150,994
Pine	404,399	1,505,507	136,205	125,384	468,058	377,293
Cleland	490,341	1,746,645	275,473	328,310	989,485	1,107,242
Protection	1,044,610	3,948,580	322,455	271,926	1,091,980	839,139
Destruction	121,136	429,640	87,846	7,867	350,103	65,891
Farallones	69,851	324,594	66,589	30,521	288,058	118,205
Año Nuevo	136,904	680,534	67,809	120,074	292,267	433,373

Colony refers to the UD derived from all birds at colony, and Individual refers to UDs of individual birds. Colonies are sorted from north to south (Fig 1).

Overall, the UDs from different colonies showed moderate spatial overlap during the winter months, with a mean UDOI of 0.45 across all pairwise values (Table 6). Spatial overlap during the winter months had a strong regional component, with the three Alaskan colonies all sharing much of the same area (UODIs ranging from 0.65 to 0.92). Similarly, the two colonies in California had an UDOI of 1.10, indicating they shared the same high use areas. Colonies from central and northern British Columbia (Lucy, S’Gang Gwaay, Triangle, Pine) all had high overlap with each other (UDOIs ranging from 0.83 to 1.34), and somewhat lower with Cleland Island along the west coast of Vancouver Island (UDOIs ranging from 0.54 to 0.72). Birds from the two colonies in Washington State, despite being in close proximity while breeding, appeared to have little overlap in winter (UDOI of 0.17). The wintering range of birds from Protection was similar to that of birds from British Columbia, whereas the two birds from Destruction were more strongly associated with the Pacific coast of the United States (Fig 5).

Table 6

Spatial overlap between home ranges of during wintering months (Nov-Feb) of rhinoceros auklets breeding at 13 colonies on the Pacific Coast of North America, 2014 to 2016, as calculated from the Utilization Distribution (UD) derived from geolocation of light loggers.

	Middleton	St. Lazaria	Lucy	S’Gang Gwaay	Triangle	Pine	Cleland	Protection	Destruction	Farallones	Año Nuevo
Chowiet	0.92	0.65	0.37	0.38	0.40	0.21	0.10	0.54	0.05	0.01	0.04
Middleton		0.74	0.37	0.37	0.39	0.12	0.07	0.46	0.03	0.02	0.04
St. Lazaria			0.88	0.99	0.99	0.47	0.16	1.10	0.12	0.02	0.05
Lucy				1.34	1.22	0.92	0.54	0.90	0.13	0.10	0.22
S’Gang Gwaay					1.26	0.92	0.56	0.96	0.18	0.10	0.21
Triangle						0.83	0.62	0.88	0.17	0.13	0.27
Pine							0.72	0.70	0.34	0.15	0.33
Cleland								0.31	0.34	0.46	0.78
Protection									0.17	0.03	0.12
Destruction										0.15	0.24
Farallones											1.10

Overlap measures are calculated as a UD overlap index (UDOI) between each pairwise combination of colony UDs (Fieberg and Kochanny 2005). The UDOI is 1.0 when both UDs are uniformly distributed and share 100% overlap, 0 when UDs have no overlap, and > 1.0 when the two UDs are non-uniformly distributed and share high overlap. Colonies are sorted from north to south (Fig 1).

Relationship between winter distribution and genetic structure

The genetic distance of rhinoceros auklets from the eastern Pacific colonies as measured by Fꞌst was not correlated with the distance between the breeding colonies (IBD), with a Spearman correlation coefficient of 0.08 (p = 0.290 from randomization; Fig 6A). In contrast, Fꞌst was negatively correlated to the spatial overlap among colonies during the winter, with a Spearman’s correlation coefficient of -0.22 (p = 0.056 from randomization; Fig 6B).

Fig 6

Correlations between genetic distance, as measured by standardized measure of genetic differentiation (Fꞌst), and (a) distance between breeding colonies (IBD), or (b) spatial overlap during the non-breeding season (Nov-Feb). Spatial overlap was measured by the Utilization Distribution overlap index (UDOI) between each pairwise combination of colonies of rhinoceros auklets at breeding colonies on the Pacific Coast of North America, 2014 to 2016.

Discussion

We combined light-level geolocator tracking with range-wide genetics analyses to directly test a key hypothesis often invoked to explain population genetic structuring in seabirds; specifically, that genetic differentiation is associated with spatial segregation in wintering areas [22, 28]. Our results support and build on a previous analysis [32] in showing that contemporary genetic structuring is extensive and complex in the rhinoceros auklet, a widely-distributed seabird of the temperate North Pacific Ocean. Genetic differentiation was strongest between the eastern vs. western Pacific populations, but structure also existed within both regional groups. Neither physical barriers to movement nor isolation by distance within the eastern group accounted for the pattern of structure. As predicted, geolocator tagging revealed a negative association between the extent of population differentiation between colony pairs (Fꞌst) and the extent of overlap in non-breeding season distributions.

Genetic structure and its relationship to winter distribution

The genetics component of our study, which spanned 18 breeding colonies and 424 individuals, further supports the major conclusions of the Abbott et al. study [32]. All of FST and FꞌST values, population-based PCoA, and individual-based STRUCTURE plots indicated that there is a high level of genetic differentiation between rhinoceros auklets breeding in the western vs. eastern Pacific Ocean.

Morphometric variation has been linked to genetic differentiation in marine vertebrates [71, 72], including seabirds [73, 74]. In the rhinoceros auklet, individuals breeding on colonies in the western Pacific are on average larger in linear dimensions and body mass than those in the eastern Pacific [75, 76]. On colonies involved in this study, the mean (SD) mass of 278 birds from four western Pacific colonies (all but Taijima) was 558 g (37.1), with a range of 460 to 695 g; and 502 g (32.5) for 392 birds from all 13 eastern Pacific colonies, with a range of 425 to 620 g. There was little variation in mass within either region. No subspecies are currently recognized in the rhinoceros auklet [33], but given the marked genetic and morphological differences, the eastern and western Pacific Ocean populations might meet criteria for subspecific designation.

Abundant on both sides of the temperate North Pacific Ocean, the rhinoceros auklet is virtually absent as a breeder in the Aleutian Islands and Bering Sea, where two close relatives, the tufted puffin (Fratercula cirrhata) and horned puffin (F. corniculata), are abundant [77]. There are no obvious physical barriers to dispersal over that long gap in the distribution of the rhinoceros auklet. However, both a species distribution model (S1 Table, S1 Fig) and a resistance surface map (S2 Fig) developed in conjunction with this study [78] indicate that the deep ocean habitat along the northern continental shelf (Alaska to Asia) inhibits movement between populations in the eastern and western Pacific Ocean. In agreement with those analyses, no geolocator-tagged rhinoceros auklets crossed from the western to the eastern Pacific [35], nor from the eastern to the western Pacific (this study), in the non-breeding season. Cryptic barriers to dispersal limit gene flow in other seabirds as well [79, 80]. However, in another auk species, the ancient murrelet (Synthliboramphus antiquus), many individuals do cross from the eastern to the western North Pacific over their annual life cycle [81]. And as predicted from the hypothesis that use of common wintering areas is one factor that links with gene flow among populations, no genetic differentiation was found between ancient murrelets breeding in the eastern vs. western Pacific based on mitochondrial DNA [82].

There is historical context to consider as well. In other marine vertebrates, ice cover during the late Pleistocene appears to have promoted an east-west genetic split across the North Pacific Ocean [83-85]. For the rhinoceros auklet, it is plausible that extensive ice cover over the North Pacific during the last glacial maximum ~19–14.6 kya [39] forced a southward range shift that isolated populations into eastern and western refugia, thereby leading to genetic divergence [32].

While genetic differentiation in rhinoceros auklets was strongest between the eastern and western North Pacific Ocean, there was a striking amount of structure within both regional groups as well. Based on FꞌST values, five of 10 comparisons between pairs of colonies were significant in the western Pacific Ocean, as were 41 of 66 pairwise comparisons in the eastern Pacific Ocean. That smaller-scale differentiation exists in the absence of obvious physical barriers to dispersal, and could not be attributed to isolation by distance. Such a complex pattern of genetic differentiation is unusual among seabirds [47], and there are undoubtedly many factors involved [27]. Here, we will consider the role of three primary factors that could be associated with genetic structuring in rhinoceros auklets: (1) philopatry; (2) a mainly coastal (neritic) foraging habitat; and (3) the extent of overlap in wintering areas.

Philopatry

Philopatry to breeding sites is a common behavioural trait in marine vertebrates [86, 87], including colonial seabirds [88]. Like other colonial auks, the rhinoceros auklet appears to be highly philopatric [34]. Between 1984 and 2019, over 4000 individuals, including both adults and nestlings, were banded on the colonies in British Columbia included in this study. Many hundreds have been re-encountered, but only on the colonies where they were banded [89, 90].

Nonetheless, some movement among colonies does occur. A number of rhinoceros auklet breeding colonies were newly established, or re-established, following extirpation due to human influence in Oregon and California from the 1960s to the 1990s [91-93], including the Farallon Islands [94] and Año Nuevo Island [95]. There were clues from Fst values about the scale and nature of movement in rhinoceros auklets. Only 25 of 66 pairwise comparisons among the eastern Pacific colonies were non-significant, and 16 of those 25 involved the three breeding colonies in Alaska. Of the nine remaining non-significant comparisons, eight involved nearby pairs of sites: five at the southern half of the range (California to southern British Columbia, all situated within the California Current Marine System), and three further north (central and northern British Columbia). Those observations imply, first, that individuals breeding in Alaska are more likely to move away from their natal colonies than individuals breeding elsewhere, a notion consistent with the generally large wintering ranges of Alaskan birds; and second, that movements of individuals on breeding colonies south of Alaska, when they occur, tend to be regional.

Neritic foraging habit

The majority of rhinoceros auklets breeding in the eastern North Pacific Ocean inhabit colonies in British Columbia, and the largest colonies are located along the mainland coast [44]. That distribution presumably reflects the importance of Pacific sand lance (Ammodytes personatus) to these seabirds [96-98]: the availability of that specific prey species is closely tied to the annual productivity of rhinoceros auklet colonies in British Columbia [99, 100]. The Pacific sand lance, like other species of sand lances, inhabits shallow, coastal environments [101]. The fact that deep ocean habitat across the North Pacific Ocean appears to act as a barrier to movement between eastern and western rhinoceros auklet populations provides further evidence of the neritic nature of the species.

There is evidence in several seabird taxa, including albatrosses [102], boobies [103] and penguins [79], and in some marine mammals [104], that species that forage predominantly in coastal, neritic environments tend to exhibit more extensive genetic structuring than species that forage in open-ocean, pelagic environments. That appears to be true among the auks of the North Pacific Ocean as well. Genetic structure is prominent in three coastal species, the pigeon guillemot (Cepphus columba) [105], marbled murrelet (Brachyramphus marmoratus) [106] and Kittlitz’s murrelet (B. brevirostris) [107]. But structure is less prominent in the more oceanic common murre (Uria aalge) [108], thick-billed murre (Uria lomvia) [109], ancient murrelet [82], Cassin’s auklet (Ptychoramphus aleuticus) [50], and crested auklet [110]. A potentially confounding factor is that all of the species in the latter group breed in large colonies, whereas all of the species in the former group breed solitarily or in small aggregations [33]. But the rhinoceros auklet nests in large colonies, is mainly neritic, and exhibits extensive genetic structuring. Those same traits are also characteristic of the whiskered auklet [111], a species that forages in tide-rips close to shore [111], and a species in which birds outfitted with geolocator tags attended the breeding colony year-round [112]. Neritic seabirds might tend to structure extensively because they conform to a one-dimensional, stepping-stone model of dispersal [105], and/or because they are unlikely to come into contact at sea with individuals from populations other than their own [79].

Overlap in wintering areas

Two lines of evidence in our results supported the hypothesis of a link between segregation in wintering areas and population genetic structuring. First, none of the tagged rhinoceros auklets breeding on one side of the North Pacific Ocean crossed to the other side [35, this study], and there was strong differentiation between populations in the eastern vs. the western Pacific. And second, in comparisons between pairs of colonies in the eastern Pacific, there was a negative correlation between the standardized measure of genetic differentiation (Fꞌst) and the extent of spatial overlap in wintering areas. The wintering range of rhinoceros auklets breeding on colonies in North America encompasses nearly all continental shelf waters of the eastern North Pacific Ocean. On average, spatial overlap in winter was moderate among all colonies (UDOI = 0.45), and there was some regional structure, especially for the two colonies in California (UDOI values > 1) that had similar high use areas in winter, not just broad overlap.

It is most often suggested that segregation in wintering areas is associated with genetic differentiation because segregation reduces the chance that an individual from one population will encounter foreign breeding colonies and/or individuals from other populations [22, 27]. For rhinoceros auklets, the strong differentiation between western vs. eastern Pacific populations appears to be caused in large part by a more-or-less impermeable barrier to dispersal in the form of deep ocean habitat in the Aleutian Islands and Bering Sea–perhaps equivalent in effect to differentiation in other seabirds that inhabit separate ocean basins and in which genetic exchange is extremely rare [113] or non-existent [108]. Historical isolation during the late Pleistocene also could have led to differentiation between the eastern and western Pacific Ocean populations. Not surprisingly, the differentiation within the two populations of rhinoceros auklet was more subtle, and in the eastern population at least, linked to the degree of segregation in wintering areas. That, we propose, is consistent with the idea that the probability of encountering foreign colonies and/or individuals at sea is a mechanism linking distribution to genetic structure.

Other mechanisms are also plausible. In Cook’s petrel (Pterodroma cookii), genetic differentiation was promoted by population-level differences in the timing of breeding, driven by differences in migration timing associated with habitat specialization in wintering areas [31]. But allochrony seems an unlikely primary mechanism for rhinoceros auklets. On colonies in the eastern Pacific, timing of breeding becomes later with increasing latitude, e.g., the date of laying of the first egg in the year ranged from 8 April to 11 May in 12 years on Año Nuevo Island, CA (Oikonos, unpublished data); 30 April to 7 May in 5 colony-years on Protection and Smith islands, WA [114]; 22 April to 4 May in 5 years on Triangle Island, BC [115]; and 4 to 17 May in 3 years on the Semidi Islands, AK [116]. But while there is a latitudinal trend in laying dates, there was no strict latitudinal pattern of population differentiation.

Although the Friesen, Burg, & McCoy [22] hypothesis has not been directly tested in other auks, separate investigations of winter distribution and genetic structure provide further support, and suggest areas for future research. First, the ancient murrelet is panmictic across the North Pacific Ocean [82], and geolocator tagging revealed that some individuals crossed from the eastern into the western Pacific in the winter [81]. Second, population genetic structure exists in the whiskered auklet [111], and geolocator tagging showed that individuals in this species overwinter close to breeding colonies [112]. Third, the winter distribution of most Cassin’s auklets breeding on two very large colonies in British Columbia (Triangle, Frederick islands) overlapped with that of individuals from Southeast Farallon Island in California [117, 118], and the two groups are genetically homogeneous [50]. And fourth, analysis of highly differentiated loci in thick-billed murres from five colonies in Canada’s eastern Arctic showed that individuals from the three colonies showing highest overlap in wintering distributions also were those that were the most genetically similar [119].

But there is also contradictory evidence among the auks. Genetic structure was more extensive in Atlantic Ocean populations of common murres [108] than in thick-billed murres [109], whereas geolocator tagging in the western Atlantic found that population segregation in winter was stricter in thick-billed than in common murres [120]. That apparent contradiction could signal the existence of real variation among auks, even between close relatives, or differences in the scale of the genetics vs. the tagging projects. Finally, a range-wide geolocator tagging study of Atlantic puffins (Fratercula arctica) [121] revealed a complex pattern of overlap in wintering areas, whereas a much smaller-scale study found little genetic differentiation based on allozyme patterns [122]. An investigation of population genetic structure in the Atlantic puffin to match the scale of the geolocator tagging could provide much insight.

Conservation implications

Large-scale tagging programs are providing valuable new information on the year-round habitat use of seabirds [121], with direct application to marine conservation [123]. Seabirds are particularly wide-ranging organisms of high conservation concern, they face a wide array of anthropogenic threats both on the land and at sea, and many species are experiencing population declines [124]. Geolocator tagging revealed that the rhinoceros auklet is very widely distributed throughout the year in continental shelf waters of the temperate North Pacific Ocean. Within this realm, rhinoceros auklets are vulnerable to a wide variety of anthropogenic stressors including oiling at sea [125, 126], bycatch in fishing gear [127-129], and chemical contamination of food webs [130-133]. They are also affected by oceanographic change [100, 134]. Widespread and vulnerable, the rhinoceros auklet also displays an extensive and complex pattern of population genetic structure. The maintenance of genetic diversity is critical to ensuring that species remain resilient to natural and anthropogenic stressors over the long term [135, 136], and is a particularly challenging component of biodiversity conservation [137, 138]. Thus, many issues uncovered here, for one North Pacific species, indicate the scope and complexity of the conservation challenges posed by seabirds.

Environmental layers used to develop the Species distribution model.

Contributions to the model are determined using a heuristic approach that depends on the path of the Maxent code. Permutation importance is determined by values randomly permutated along training points and measurements for the decrease in training AUC. Variables with a higher influence have a larger percent value.

(DOCX)

Click here for additional data file.

Species distribution model for rhinoceros auklet created from GBIF breeding season occurrences (May-July).

Map was produced using the SDM toolbox (Brown, 2014; Brown et al., 2017), Maxent (Phillips et al., 2006; Phillips & Dudík, 2008), and ArcMap. The final map was visualised using ArcMap 10.2 (ESRI®) in the Azimuthal Equidistant (180° meridian) projection. Most suitable habitat for is shown in cool (blue) unsuitable habitat in warm (orange to red).

(TIF)

Click here for additional data file.

Resistance surfaces for 18 rhinoceros auklet breeding colonies located throughout the North Pacific Ocean.

Areas of low resistance are in blue with gradients of orange to red indicating area of higher resistance. Breeding colonies denoted by the pink circled star. Resistance surfaces were developed using friction surfaces from an inverted SDM for rhinoceros auklet and least-cost corridors function within the SDM toolbox (Brown, 2014; Brown et al., 2017). The final map was visualised using ArcMap 10.2 (ESRI®) in the Azimuthal Equidistant (180° meridian) projection.

(TIF)

Click here for additional data file.

21 May 2020

PONE-D-20-09231

Geolocator tagging links distributions in the non-breeding season to population genetic structure in a sentinel North Pacific seabird

PLOS ONE

Dear Dr. Hipfner,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

We would appreciate receiving your revised manuscript by Jul 05 2020 11:59PM. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Janice L. Bossart

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. In your Methods section, please state the volume of the blood samples collected for use in your study.

3. In your Methods section, please provide additional information regarding the permits you obtained for the work. Please ensure you have included the full name of the authority that approved the field sites access and, if no permits were required, a brief statement explaining why.

4. We note that you have stated that you will provide repository information for your data at acceptance. Should your manuscript be accepted for publication, we will hold it until you provide the relevant accession numbers or DOIs necessary to access your data. If you wish to make changes to your Data Availability statement, please describe these changes in your cover letter and we will update your Data Availability statement to reflect the information you provide.

Additional Editor Comments (if provided):

Please pay close attention to and address concerns raised by both reviewers. Both make legit points. Major concerns that need to be addressed are:

1) General issues regarding IBD, e.g. whether appropriate or not across the entire range, and needing to be introduced and conceptually developed within the Introduction.

2) Correlation vs. causation.

3) Clarity regarding breeding vs. non-breeding/winter areas, timing of such, etc... e.g. as identified by Reviewer 2. To the extent that genetic structure is a consequence of gene flow, which of course is tied directly to mating, there needs to be better description of these colonies and where/when breeding takes place versus where/when they are relative to this when not breeding. Like Reviewer 2, I was confused. The reader is provided insufficient information to understand the situation.

In addition to the other issues the two reviewers identified, below is a list of minor problems I found.

Throughout, replace the apostrophe in F'ST with a straight apostrophe (e.g. using insert symbol in Word).

Figure 2a & 2b. Double check the ordering/naming of the colonies. To the extent, n is larger for MAT than TEU, these two, at least, apparently switched order between a & b. I also wondered about DE & PR. Be 100% sure that all names are associated with the correct genetic info in the STRUCTURE diagram.

Also, Figure 2a & 2b. Both should look identical structurally given they are essentially part of the same figure: 2a currently has more space and an additional solid black line on the left. The region names and brackets also need to better match the associated breeding colonies. Currently, TOD looks to be part of Alaska, especially in 2b. The brackets and region names in general appear to be less carefully aligned/applied in 2b versus 2a. There's also no reason the bracket for Japan should extend well beyond the colony info. Making things exact is easy using PowerPoint by inserting figures, then adding text/brackets and aligning using grid lines. Seems to me 2a and b could have easily been pre-assembled into a single jpeg.

Seems to me given their geographic location and their genetic similarity, SGG and TRI should be ordered together in all STRUCTURE figures/diagrams.

In general, bracket lines on all STRUCTURE figures should line up with the vertical bracket lines in the STRUCTURE diagram.

Figure 4. Relevance of colors and legend. Change label for Protection to PR.

Figure 6. Correlation coefficients are generally designated using 'r' not 'p'. Until I read the text, I thought the figure included p-values of non significance. For clarity, replace 'p' with 'r'.

Line 393 and 394. Replace Table 4 with Table 3.

Line 419. Replace 'Moore, Cleland' with 'Moore and Cleland'

Line 420. Identify which cluster

Line 436-440. Given that you're dealing with Q>60% throughout, use it once only... or at most twice. E.g. The eastern Pacific colonies split into two clusters at Q>60%. Cluster one included... etc. Cluster two included... In cluster one, delete 'Protection'. In both cases, precede the last colony in the list by 'and'. Why were SEF and MID singled out given all member of cluster 1 and cluster 2 were associated with the respective clusters at Q>60%?

Line 447. Delete Protection. The text says 60%. Which is it 50 or 60%?

Line 448. Delete 'and' before Protection. Replace '&' with 'and' before Ano Nuevo. The text says 60%. Which is it 50 or 60%?

Line 520. Figure 6b says -0.21. Which is it?

Lines 641-643. Given the amount of spread, i.e. low correlation, and fact the correlation was only significant at 0.05, I think this is a bit overstated.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: This paper contributes to our understanding of the factors promoting population differentiation in seabirds and other species. I really enjoyed reading this paper. It's written very well and and it is clear and easy to follow. The science is sound and the methods are appropriate to test the hypotheses presented in the introduction. The conclusions are supported by the data and are well presented and discussed in a comparative framework as well. I have only very minor edits to suggest.

371-372: negative FST values are an artifact of some programs and should be replaced by zeros.

389-390: A significant Mantel’s test indicates correlation between genetic and physical distance but this doesn’t translate necessarily into isolation by distance. The fact that you get a significant Mantel’s test when you include populations from both sides of the Pacific but not when you analyze eastern Pacific colonies alone (lines 516-518) indicate that in the former the pattern is driven by strong pop structure between eastern and western Pacific colonies rather than by IBD. See Meirmans, P.G. 2012. The trouble with isolation by distance. Molecular Ecology 21: 2839-2846.

516-518: I would say that all results from Mantel’s should be presented together but testing for IBD between eastern and western Pacific colonies may be inappropriate (see above).

547-548: See comments above. This is structure between two differentiated groups rather than IBD.

Line 682: Analysis of highly differentiated loci in thick-billed murres (Tigano et al. 2017) has revealed a similar pattern, although not quantitively. The three Atlantic colonies that were more genetically similar at these outlier loci were also the colonies showing the highest overlap in wintering distributions. However, this pattern could have also been explained by colonization history following the retreat of the Pleistocene ice sheets or selection during the non-breeding time.

Figure 1. Please show all colonies included in the study and highlight those that have also GLS data

Figure 2-4 Structure plot image quality is low and pop labels are hard to read. Make sure that images are high quality.

Fig. 4 What do colors represent? Add legend.

Fig. 6 Add p-values on graphs and/or figure caption.

Reviewer #2: This paper makes an important contribution to our understanding of population structure, distribution, and space use of the Rhinoceros Auklet. The study is particularly valuable because it aggregates data from breeding colonies across nearly the entire distribution of the species, which will be critical baseline information for management and conservation of the species. The data are robust and well analyzed. Overall this is a well researched and well written paper!

My only real issue with the paper is that the abstract, intro, and discussion use language that makes it appear the authors found a causal relationship when they only tested a correlation. The study is framed as a test of the hypothesis that segregation in the wintering distributions drives genetic differentiation among breeding colonies. However, the paper does not test a causal mechanisms by which nonbreeding segregation would drive genetic differentiation, and as a paper the authors cite (Rayner et al 2011) points out, “Spatial segregation during the non-breeding season appears to provide an intrinsic barrier to gene flow among seabird populations that otherwise occupy nearby or overlapping regions during breeding, but how this is achieved remains unclear.” What the study does do is test if there is a correlation between genetic differentiation and degree of difference in space use during the nonbreeding season. I think the language can easily be revised to clarify that this distinction by removing words like “promotes” and “drives” that imply causation and replacing them with words like “correlates with” that make it clear what the authors are testing.

The authors only set up one hypothesis test in the intro (winter segregation is correlated with genetic structure), but they test an additional hypothesis in the paper (distance between breeding colonies is correlated with genetic structure). I assume the authors test this to evaluate if average dispersal difference accounts for genetic structure, assuming all populations have equal mean dispersal differences (which is not a given). Since the authors present results of both tests and use them to say that overwintering overlap better accounts for genetic difference than mean isolating distance, both hypothesis tests should be set up in the intro. The authors also discuss two additional hypotheses in the discussion: genetic structure correlates with egg lay date or morphological differences like body size or mass (potentially indicating differences in foraging niche). These are also testable hypotheses, and it appears the authors have the data to test them instead of just discussing them. I think this could potentially strengthen the paper, though it’s not essential.

A limitation of the study is that the authors did not consider temporal segregation, only spatial segregation. Because of that, it is unclear to what degree these breeding colonies are potentially in contact during the nonbreeding season. The geolocator study cited from Japan (Takahashi et al. 2015) did find that space use varied with month, so I’m curious why the authors did not consider or account for that. I think the authors should address this in the discussion at a minimum.

Authors use inconsistent terminology for geolocators. I found the following terms used in the ms: geolocator, geologger, logger, and GLS device. Bird carrying geolocators were also referred to as loggered and geologgered. Best to pick one and use it consistently to avoid confusion.

325: When do the colonies initiate and end breeding? How does that relate to the “winter” dates selected?

365, comma missing

631-633: I found this sentence confusing. I recommend “in which the population showed the unusual behavior of remaining in close proximity to a breeding colony year-round and foraging…”

640-643: change word from “promotes” as you don’t have evidence of a mechanistic relationship.

658: not “driven,” correlated with. Again, no causation has been shown.

669-670: These are alternative hypotheses that are testable.

Table 4: I think there’s an error in Middleton, which shows more track files recovered than geolocators recovered.

Fig 1: As this is the only full study area map, the Western Pacific sites should be included and differentiated with a different color or icon.

Fig 4: Unclear what the colors represent. Please include a color legend. It would be logical for colors to be grouped by the five regions identified in figures 2 and 3, but there are more than five colors presented here. Please clarify.

Figure 5: The resolution of the file attached for figure 5 is low, so it’s hard to figure out what the gray lines are. The caption says they represent the geolocator tracks, but the look like dashed, horizontal lines, whereas geolocator tracks are typically presented as a single connected line. The gray lines here look more like shading, but it’s unclear how that shading was generated. Does this incorporate some metric of uncertainty around each point in the track or aggregate all tracks in some way? Please clarify.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Anna Tigano

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

1 Sep 2020

Journal Requirements:

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

We have carefully reviewed the document to ensure complaince with PLoS ONE format.

2. In your Methods section, please state the volume of the blood samples collected for use in your study.

We added the following information: Blood samples (1.0 ml) and/or feather samples (5-10, from the breast) were collected from seven to 80 individuals on each of the 18 colonies, for a total of 704 individuals [line 225].

3. In your Methods section, please provide additional information regarding the permits you obtained for the work. Please ensure you have included the full name of the authority that approved the field sites access and, if no permits were required, a brief statement explaining why.

We have provided all information on permits to access field sites and conduct the research [lines 211-219].

4. We note that you have stated that you will provide repository information for your data at acceptance. Should your manuscript be accepted for publication, we will hold it until you provide the relevant accession numbers or DOIs necessary to access your data. If you wish to make changes to your Data Availability statement, please describe these changes in your cover letter and we will update your Data Availability statement to reflect the information you provide.

We have deposited the microsatellite data in Dryad and the GLS tagging data in Movebank.

Comments from the Editor

1) General issues regarding IBD, e.g. whether appropriate or not across the entire range, and needing to be introduced and conceptually developed within the Introduction.

We agree with Reviewer #1 that our test for IBD across the full range of the rhinoceros auklet was not appropriate, and have removed this. We have retained the test for IBD among the eatsern Pacific Ocean colonies (finding no support), and have briefly introduced the topic in the Introduction using the helpful phrase from Reviewer #2 [lines 165-167].

2) Correlation vs. causation.

We made this change throughout the manuscript, replacing “promotes” and “drives” with “is linked to” and “is associated with”.

3) Clarity regarding breeding vs. non-breeding/winter areas, timing of such, etc... e.g. as identified by Reviewer 2. To the extent that genetic structure is a consequence of gene flow, which of course is tied directly to mating, there needs to be better description of these colonies and where/when breeding takes place versus where/when they are relative to this when not breeding. Like Reviewer 2, I was confused. The reader is provided insufficient information to understand the situation.

We regret that we were not as clear as we could have been on this point. We have added a few sentences to the Methods section to clarify why we made the decision to examine the period from 1 November to 28/29 February [lines 355-358].

Throughout, replace the apostrophe in F'ST with a straight apostrophe (e.g. using insert symbol in Word).

We made this change throughout.

Figure 2a & 2b. Double check the ordering/naming of the colonies. To the extent, n is larger for MAT than TEU, these two, at least, apparently switched order between a & b. I also wondered about DE & PR. Be 100% sure that all names are associated with the correct genetic info in the STRUCTURE diagram.

We checked all of this carefully, and believe all of these issues have been addressed.

Also, Figure 2a & 2b. Both should look identical structurally given they are essentially part of the same figure: 2a currently has more space and an additional solid black line on the left. The region names and brackets also need to better match the associated breeding colonies. Currently, TOD looks to be part of Alaska, especially in 2b. The brackets and region names in general appear to be less carefully aligned/applied in 2b versus 2a. There's also no reason the bracket for Japan should extend well beyond the colony info. Making things exact is easy using PowerPoint by inserting figures, then adding text/brackets and aligning using grid lines. Seems to me 2a and b could have easily been pre-assembled into a single jpeg.

We have made extensive changes to the figures, following the recommendations of all reviewers. We believe they are now substantially improved.

Seems to me given their geographic location and their genetic similarity, SGG and TRI should be ordered together in all STRUCTURE figures/diagrams.

We respectfully disagree with this suggestion. The colony at TRI is closer to PI than it is to SGG, and similar to the distances between several other BC colonies. And while physically close, the two colonies are in quite different oceanographic zones – TRI in the transition zone between the California and Alaska Currents, but SGG in the Alaska Current. Finally, the PCoA provided some (weak) evidence that the colony at TRI is genetically distinct from all other colonies, as we concluded in our earlier (and less comprehensive) analysis (Abbott et al. 2014).

In general, bracket lines on all STRUCTURE figures should line up with the vertical bracket lines in the STRUCTURE diagram.

We believe these issues have now been addressed.

Figure 4. Relevance of colors and legend. Change label for Protection to PR.

We changed the label for PR, and added this sentence to the Figure caption: “The colours correspond to the three groups detected using STRUCTURE (Figs 2 and 3): western Pacific (yellow), the larger eastern Pacific group (purple), and smaller eastern Pacific group (blue).”

Figure 6. Correlation coefficients are generally designated using 'r' not 'p'. Until I read the text, I thought the figure included p-values of non significance. For clarity, replace 'p' with 'r'.

Those symbols were meant to be the lower case Greek letters rho, the symbol for Spearman’s correlation coefficient. But we have learned that this symbol is used as the coefficient for the entire population, and not for a sample. We have made this change in Fig. 6.

Line 393 and 394. Replace Table 4 with Table 3.

We made this change [lines 421-422].

Line 419. Replace 'Moore, Cleland' with 'Moore and Cleland'

We made this change [line 449].

Line 420. Identify which cluster

That should have read “cluster two” rather than “one cluster”. We made the change [line 450].

Line 436-440. Given that you're dealing with Q>60% throughout, use it once only... or at most twice. E.g. The eastern Pacific colonies split into two clusters at Q>60%. Cluster one included... etc. Cluster two included... In cluster one, delete 'Protection'. In both cases, precede the last colony in the list by 'and'. Why were SEF and MID singled out given all member of cluster 1 and cluster 2 were associated with the respective clusters at Q>60%?

We went through all of the STRUCTURE results very carefully, and think the text and Figure captions are now accurate and succinct [lines 441-481].

Line 447. Delete Protection. The text says 60%. Which is it 50 or 60%?

See comment above.

Line 448. Delete 'and' before Protection. Replace '&' with 'and' before Ano Nuevo. The text says 60%. Which is it 50 or 60%?

See comment above.

Line 520. Figure 6b says -0.21. Which is it?

We made this change - -0.22 is correct [line 554].

Lines 641-643. Given the amount of spread, i.e. low correlation, and fact the correlation was only significant at 0.05, I think this is a bit overstated.

We thought carefully about how to deal with the marginal correlation. In general, we believe a correlation at (two-tailed) P = 0.056 in the direction predicted by the hypothesis being tested should be considered statistically significant. Our resident biostatistician (M. Drever) also noted that given dispersion in the data, we are justified in using P = 0.1 to avoid Type II error. We added this to the Methods section [lines 387-388].

Reviewer #1

371-372: negative FST values are an artifact of some programs and should be replaced by zeros.

We made this change to text and tables.

389-390: A significant Mantel’s test indicates correlation between genetic and physical distance but this doesn’t translate necessarily into isolation by distance. The fact that you get a significant Mantel’s test when you include populations from both sides of the Pacific but not when you analyze eastern Pacific colonies alone (lines 516-518) indicate that in the former the pattern is driven by strong pop structure between eastern and western Pacific colonies rather than by IBD. See Meirmans, P.G. 2012. The trouble with isolation by distance. Molecular Ecology 21: 2839-2846.

We agree with the reviewer, and thank her or him for the advice. We removed the test for IBD across all western and eastern North Pacific colonies.

516-518: I would say that all results from Mantel’s should be presented together but testing for IBD between eastern and western Pacific colonies may be inappropriate (see above).

We agree, and have removed any testing of IBD across the full range.

547-548: See comments above. This is structure between two differentiated groups rather than IBD.

We agree.

Line 682: Analysis of highly differentiated loci in thick-billed murres (Tigano et al. 2017) has revealed a similar pattern, although not quantitively. The three Atlantic colonies that were more genetically similar at these outlier loci were also the colonies showing the highest overlap in wintering distributions. However, this pattern could have also been explained by colonization history following the retreat of the Pleistocene ice sheets or selection during the non-breeding time.

We added this result, and the Tigano et al. (2017) reference [lines 715-718].

Figure 1. Please show all colonies included in the study and highlight those that have also GLS data

We initially produced a figure showing all of the colonies, both in Japan and in North America. But this necessitated a figure of unacceptably large scale, creating issues around detail and visual separation of colonies. This was especially true of the western Pacific colonies - relative to the eastern Pacific, they are very close together. So, given that we did no geolocator tagging at the colonies in Japan for this study, that the focus of the work in Japan was to demonstrate the very strong east vs. west Pacific genetic split, we thought it better to go with a high quality map of the North American colonies rather than a low quality map showing all colonies. Locations of the sites in Japan are available in the Appendices, however.

Figure 2-4 Structure plot image quality is low and pop labels are hard to read. Make sure that images are high quality.

We made many changes to the figures, and believe we have addressed the comments of all three reviewers.

Fig. 4 What do colors represent? Add legend.

We added a sentence to the figure caption to explain the relevance of the different colours [lines 499-500].

Fig. 6 Add p-values on graphs and/or figure caption.

We made these additions.

Reviewer #2

My only real issue with the paper is that the abstract, intro, and discussion use language that makes it appear the authors found a causal relationship when they only tested a correlation. The study is framed as a test of the hypothesis that segregation in the wintering distributions drives genetic differentiation among breeding colonies. However, the paper does not test a causal mechanisms by which nonbreeding segregation would drive genetic differentiation, and as a paper the authors cite (Rayner et al 2011) points out, “Spatial segregation during the non-breeding season appears to provide an intrinsic barrier to gene flow among seabird populations that otherwise occupy nearby or overlapping regions during breeding, but how this is achieved remains unclear.” What the study does do is test if there is a correlation between genetic differentiation and degree of difference in space use during the nonbreeding season. I think the language can easily be revised to clarify that this distinction by removing words like “promotes” and “drives” that imply causation and replacing them with words like “correlates with” that make it clear what the authors are testing.

We completely agree, and addressed this issue in replies to the Editor’s comments.

The authors only set up one hypothesis test in the intro (winter segregation is correlated with genetic structure), but they test an additional hypothesis in the paper (distance between breeding colonies is correlated with genetic structure). I assume the authors test this to evaluate if average dispersal difference accounts for genetic structure, assuming all populations have equal mean dispersal differences (which is not a given). Since the authors present results of both tests and use them to say that overwintering overlap better accounts for genetic difference than mean isolating distance, both hypothesis tests should be set up in the intro.

We completely agree, and addressed this issue in comments to the Editor above

The authors also discuss two additional hypotheses in the discussion: genetic structure correlates with egg lay date or morphological differences like body size or mass (potentially indicating differences in foraging niche). These are also testable hypotheses, and it appears the authors have the data to test them instead of just discussing them. I think this could potentially strengthen the paper, though it’s not essential.

First, there is, as we discuss, a clear east-west difference in body size, but the larger birds from Asia do not interact or cohabitate with the smaller birds from North America - the two phenotypes are completely segregated throughout the year. As for structure being due to differences in timing of breeding, we did briefly test that hypothesis and concluded that it is unlikely to underlie the population genetic structuring in rhinoceros auklets, as it does in Cook’s petrel. We then provided the evidence we used to arrive at that conclusion - there is a latitudinal trend in timing of egg-laying (typically about one month later in the north than in the south, but with much variation), but the genetic variation does not match that latitudinal pattern [lines 698-705].

A limitation of the study is that the authors did not consider temporal segregation, only spatial segregation. Because of that, it is unclear to what degree these breeding colonies are potentially in contact during the nonbreeding season. The geolocator study cited from Japan (Takahashi et al. 2015) did find that space use varied with month, so I’m curious why the authors did not consider or account for that. I think the authors should address this in the discussion at a minimum.

We addressed this issue in comments to the Editor above.

Authors use inconsistent terminology for geolocators. I found the following terms used in the ms: geolocator, geologger, logger, and GLS device. Bird carrying geolocators were also referred to as loggered and geologgered. Best to pick one and use it consistently to avoid confusion.

We have addressed this issue – after initially introducing “light-level geolocator (GLS) tags”, we refer in the rest of the manuscript to “geolocator tags”, “tags” and “tagging”

325: When do the colonies initiate and end breeding? How does that relate to the “winter” dates selected?

We addressed this issue in replies to the Editor.

365, comma missing

We added the comma.

631-633: I found this sentence confusing. I recommend “in which the population showed the unusual behavior of remaining in close proximity to a breeding colony year-round and foraging…”

That amended that sentence [lines 665-667].

640-643: change word from “promotes” as you don’t have evidence of a mechanistic relationship.

We addressed this issue, as outlined in the replies to the Editor’s comments.

658: not “driven,” correlated with. Again, no causation has been shown.

We addressed this issue, as outlined in the replies to the Editor’s comments.

669-670: These are alternative hypotheses that are testable.

We addressed this issue as described above [lines 698-705].

Table 4: I think there’s an error in Middleton, which shows more track files recovered than geolocators recovered.

We corrected the error.

Fig 1: As this is the only full study area map, the Western Pacific sites should be included and differentiated with a different color or icon.

We addressed this issue in response to reviewer #1. The full map is in the Appendices.

Fig 4: Unclear what the colors represent. Please include a color legend. It would be logical for colors to be grouped by the five regions identified in figures 2 and 3, but there are more than five colors presented here. Please clarify.

We addressed this issue in response to reviewer #1.

Figure 5: The resolution of the file attached for figure 5 is low, so it’s hard to figure out what the gray lines are. The caption says they represent the geolocator tracks, but the look like dashed, horizontal lines, whereas geolocator tracks are typically presented as a single connected line. The gray lines here look more like shading, but it’s unclear how that shading was generated. Does this incorporate some metric of uncertainty around each point in the track or aggregate all tracks in some way? Please clarify.

Figure 5 is now provided at an increased resolution of 800 ppi. The TwilightFree algorithm provides the most likely daily point locations on a pre-determined grid of 1°×1° cell sizes, based on the each bird’s light record for that day. These underlying grid makes the points appear as lines, but they are point locations of all birds from each colony, with no uncertainty depicted. We provided those points to show the entire wintering distributions of all birds for each colony. The Utilization Distribution (UD) polygons are then derived from these points.

Submitted filename: Cover letter.doc

Click here for additional data file.

18 Sep 2020

Geolocator tagging links distributions in the non-breeding season to population genetic structure in a sentinel North Pacific seabird

PONE-D-20-09231R1

Dear Dr. Hipfner,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Janice L. Bossart

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

26 Oct 2020

PONE-D-20-09231R1

Geolocator tagging links distributions in the non-breeding season to population genetic structure in a sentinel North Pacific seabird

Dear Dr. Hipfner:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Janice L. Bossart

Academic Editor

PLOS ONE