Strong Philopatry, Isolation by Distance, and Local Habitat Have Promoted Genetic Structure in Heermann’s Gull connection between dispersal and genetic

: Philopatry can promote genetic differentiation among populations but remains un-described in many seabirds. Hence, we explored such associations in Heermann’s Gull. Philopatry was observed monitoring 998 gulls in Rasa Island, while genetic differences were examined in the Cardonosa, Rasa, and Isabel islands using the cytochrome b of 296 gulls. Adults returned repeatedly to its natal valley or to a very close distance from it under different modelled hypotheses. Likewise, the interaction between sex and distance indicated significant male-biased philopatry. Besides, low to high genetic differentiation was observed between the Rasa and Cardonosa islands ( Φ ST = 0–0.22) (both in the Midriff Islands Region), but higher genetic differentiation against Isabel Island ( Φ ST > 0.25) (in the Mexican Province region). Consistently, genetic structure among regions was observed using different approaches (AMOVA: Φ CT = 0.49; SAMOVA: F CT = 0.49; and BAPS: K = 2). Similarly, a pattern of isolation by distance (rM= 0.82, p = 0.03), agrees with lower estimates of scaled migration rates between regions than among islands of the same region. Overall, it is suggested that the genetic structure found in Heermann’s Gull has been promoted by physical and behavioral barriers.


Introduction
Philopatry, sensu stricto, is the return of individuals to their natal site for breeding while dispersal is the opposite, and, in its original sense, a high degree of philopatry implies minimal gene flow [1][2][3]. In turn, gene flow means individuals moving into and out of a population, influencing the partitioning of genetic variation among populations [4]. The pattern of movement of individuals can profoundly affect the genetic structure of populations; as generations go by, the fact that individuals are being born and breed in the same places should result in the formation of discrete populations, suggesting that the stronger the philopatry, the stronger the genetic structure [3,5]. If so, genetic signatures thereof can persist among separated populations within widely distributed and highly mobile species [6][7][8]. Nonetheless, not all the birds are philopatric, and, by extension, not all show reduced gene flow, as this phenomenon can vary significantly per species, and even within species [9,10]. To better describe the connection between dispersal and genetic  Table S1) were banded within their natal valleys (Orange continuous lines: ES, Valle del Estero; IT, Vallecito del Estero; CA, Valle de la Casita; GE, Gran Estación; VG, Valle de los Gallitos; TV, Tapete Verde; CR, Casitas Viejas Arriba; CB, Casitas Viejas Abajo; and VW, Valle de la W). Later, they were recaptured breeding as adults (n = 998)) see Table 1. Mitochondrial cytochrome b sequences (1033 bp) were obtained from 296 gulls in nests from twenty locations distributed on Cardonosa, Rasa and Isabel Islands (dark circles, numbered from 1 to 20) see Table 2. Map modified from Mellink [51], showing the classic regions Midriff Islands Region and Mexican Province [54], as well as the Large Marine Ecosystems [53]. Photo credit: Misael Mancilla. Table 1. Heermann's Gull individuals banded and recaptured on Rasa Island (n = 998). Individuals were banded at "Natal Valleys" and recaptured at either the same valley or at other "Recapture Valleys": IT, Vallecito del Estero; CA, Valle de la Casita; GE, Gran Estación Central; TV, Tapete Verde; CR, Casitas Viejas Arriba; CB, Casitas Viejas Abajo; and VW, Valle de la W. Valle Norte, Valle del Estero, and Valle de los Gallitos valleys were not monitored after banding.

Region Island Cluster Location Lat (N)/ Long (W) n (M, F) h S k Hd
Cardonosa Island is in the Midriff Islands Region of the Gulf of California (28°53′16″ N and 113°01′51″ W) [46,49]. It is an irregularly flattened island, 483 m long, 292 m wide, 18 m high, and ~0.14 km 2 in area [43]. The vegetation consists mainly of the coastal Saltbush (Atriplex barclayana), the Baja California cholla (Cylindropuntia alcahes), and the cardon (Pachycereus pringlei) [25,54,55]. Here, the Heermann's Gulls nest in two places: up in an eroded flat area covered by guano and down on a boulder and pebble beach ( Figure  S1a,b). We sampled the second location in 'Playa Rocosa', a rugged area with scarce vegetation ( Figure S1c). Rasa Island is located 4 nautical miles SE of Cardonosa, in the same biogeographic region (28°49′29″ N and 112°58′50″ W) [48][49][50]. It is a young flattened volcanic island (~10,000 years), 1170 m long, 703 m wide, 30 m high, and ~0.68 km 2 in area [43,54]. The relief consists of low flat areas covered by guano (valleys) and rocky slopes (hills), with relatively little change in the last fifty years [43,47,54] (Figure S2a-c). The vegetation consists of salt scrub in three esteros and desert scrub in the flats and hill areas, the last with a notable increase of chain-fruit cholla (Cylindropuntia fulgida) [47,54]. On Rasa Island, we sampled mainly in valley areas and surrounding rocky hills. Predation by the Yellowfooted Gull (Larus livens) ( Figure S2d) and Common Raven (Corvus corax) has been observed in Heermann's Gull chicks and eggs, and Peregrine Falcon (Falco peregrinus), also preys on adults and chicks of Heermann's Gulls [43,47,[55][56][57][58]. In normal years, initially, thousands of Heermann's Gulls arrive at the nesting area, and, a few days later the Royal Tern (Thalasseus maximus) and the Elegant Tern (Thalasseus elegans) arrive in large numbers, displacing many gulls from their valley nesting areas and establishing their nests in a frequently mixed colony [44] ( Figure S2a).
Isabel Island is located off the coast of Nayarit, in the Mexican Province (21°50′40″ N and 105°52′49″ W) [48,50]. This island is 1570 m long, 100-850 m wide, 30-85 m high, and ~1.06 km 2 in area [43]. Isabel Island has cliffs with steep slopes, rocky and sandy beaches, and vegetation patches with dense grass and trees areas of low deciduous tropical forest with strong seasonality [59]. Here, we found an area with approximately fifty scattered nests of Heermann's Gull near Brown Booby (Sula leucogaster) nests, over a mixed terrain of grass and solid tuff-breccia rock [59] (Figure S3a,b). Gull chicks were reared more successfully on nests among vegetation than in hard rock, and it is relatively common to find adult gulls stalking nests and consuming Blue-footed Booby (Sula nebouxii) eggs ( Figure  S3c), whereas the insular Milk-snake (Lampropeltis triangulum sinaloae) may also feed on gull chicks ( Figure S3d).

Analyses of Philopatry and Dispersal by Banding data
The Heermann's Gull nested synchronously in large numbers and was monitored on Rasa Island as part of other studies [44][45][46]60]. The banding was carried out in large flocks of fledglings moving together and restricted in chicken-wire enclosures in each natal valley ( Figure 1). Inside the enclosure, chicks were captured by hand by four banders, banded, and immediately released to the outside of the enclosure. Disturbance was avoided by waiting till chicks had reached fledgling stage, since fledglings are already moving freely from parents and territory at that age, returning to their territory only to be fed by parents. We recaptured all the individuals included in the analysis, and band identification was based on numerical characters, not colors. The recaptures were performed using walk-in traps placed on their nests [60]. Sex assignment of recaptured Heermann's Gull adults was carried out through mating behavior and somatic measurements (i.e., within members of each pair, males are heavier [60] and larger [61] than females).
To avoid duplicating records of individuals observed in more than one season, our calculations considered gulls that were observed breeding only for the first time they were recaptured. We organized a data matrix containing (a) the natal valley (where the banded bird was born), (b) the breeding valley (where the recovered adult was found nesting), (c) the distance between the natal and the breeding valleys (with reference to their location in Rasa Island; Figure 1 If strong natal philopatry occurs in the Heermann's Gull, we would expect nesting banded birds to be recaptured in, or close to, their natal valley, i.e., the recovery of nesting banded birds would be a negative function of distance to the natal valley [3,10]. Because different number of birds were banded in each valley, we analyzed the proportion of birds recaptured in each breeding valley that had been originally banded in each natal valley. The analysis was done through a Generalized Linear Model for binomial proportions, with logit link and Chi-squared deviance. Using a stepwise additive model, we tried as predictors of the proportion of birds breeding in a given valley: (a) the distance between natal and breeding valleys, (b) the fixed effect of the individual's sex, (c) the interaction between distance and sex, (d) the effect of re-sampling effort in the breeding valley, (e) the fixed of the natal and breeding valleys, and, lastly, (d) the factorial interaction between natal and breeding valleys.
Because at the time of banding the sex of chicks cannot be determined, we did not know the number of originally banded chicks in each sex in order to calculate sex-specific proportions of banded-to-recaptured. For this reason, we performed three different analyses under different assumptions to examine whether sexual bias is occurring (males prone to philopatry and females to dispersal). Model I: We first analyzed the proportion recaptured birds in each breeding valley in relation to the total number of recaptured birds from each natal valley. Model II: We assumed that the original sex ratio in the banded birds had been similar to the ratio found in the recaptured birds. Model III: Lastly, we also performed the analysis assuming that the original sex ratio in the banded birds had been 50%-50%. All analyses were done in the R package [62], using the function glm (Generalized Linear Model) with binomial error and logit link (see the logistic model script in Supplementary Material S3).

DNA Extraction, PCR Amplification, Sequencing, and Analyses of Genetic Diversity
DNA extraction, PCR amplification and sequencing was performed only with samples from Isabel Island (n = 10), which in turn were the most recently collected in 2017. All other sequences were those from a previous study by Ruiz et al. [49] and downloaded from GenBank (see below). Mitochondrial cytochrome b (Cyt-b) sequences were obtained from DNA isolated from ~60 µL of blood from the ulnar vein of Heermann's Gulls using the DNeasy ® Blood & Tissue kit (QIAGEN, Hilden, Germany). PCR reactions were performed in a final volume of 25 µL containing 50 ng of genomic DNA (1 µL), 10 mg mL-1 BSA (1 µL), DMSO 2% (1 µL), deionized water (10 µL), 2× Kapa Taq ReadyMix (Kapa Biosystems, R&D Cape Town, ZA) (10 µL), 10 mM primers L14841 (1 µL), and 10 mM primers H16065 (1µL) [49]. The PCR program and primers were previously described by Ruiz et al. [49], and the PCR products were outsourced for sequencing to Macrogen Inc. (Seoul, Korea). Then, the recovered sequences (1033 bp) were aligned and reviewed using Seaview v.4.5.4 [63]. Next, we identified the sequences as Cyt-b, by comparing them with previous GenBank sequences from Rasa and Cardonosa (accession numbers KT275317-KT275602; [49]). Finally, newly obtained sequences from Isabel Island were submitted to GenBank (accession numbers MT935687-MT935696). Jukes-Cantor 69 was selected as the best model fitting to molecular evolution by jModelTest v.2.1.7 [64]. Haplotype richness (h), number of segregating sites (S), number of nucleotide differences (k), haplotype diversity (Hd), and nucleotide diversity (π) were estimated using DnaSP v.6.12.01 [65]. Relationships among haplotypes from each island were inferred with a minimum spanning network by PopART v.1.7 [66]. To our knowledge, there is only another Cyt-b sequence for Larus heermanni from Grays Harbor, USA [67,68]. However, this was excluded from our analysis in the main text because it was shorter (290 bp) and the different diversity indices cannot be recovered from a single sequence. Instead, a phylogeny incorporating all the Cyt-b sequences reported to date is shown in the Supplementary Material (Table  S2 and Figure S4). Phylogenetic relationships were assessed using the automatic selection of substitution model [69], subtree pruning and regrafting (SPR) tree searching [70], and aLRT SH-like branch support [71], all the tools implemented in PhyML 3.0 [72].
On the other hand, to assess their significance, the pairwise ΦST values among the twenty locations were examined with an exact test of population differentiation, with 100,000 steps of Markov chains and 10,000 permutations, using Arlequin v.3.5.2.2 [83].

Genetic Structure
To assess population structure in the Heermann's Gull, three approaches were followed. In all of them, the twenty locations sampled were regarded as separated populations ( Figure 1). First, we assumed an a priori definition by analyzing molecular variance (AMOVA), which was used to test the statistical significance of the genetic partition within groupings inferred by geography and proximity (regions, islands, and clusters, as below in Table 2), with the Φ-statistics from haplotypic data using Arlequin v.3.5.2.2 [83]. Next, we used an a posteriori approach (SAMOVA) to search for discontinuity between the twenty locations; this analysis indirectly detects barriers and defines groups maximally differentiated [84]. With this second approach, the number of groups (K = n − 1) was identified using the highest FCT value. Significance tests were performed with 1000 permutations using SAMOVA v.2.0 [84]. Finally, a Bayesian analysis of population structure based on linked molecular information was carried out through clustering algorithms to assign individuals to clusters [73,85]. This third approach tests the population structure across populations employing a clustering model implemented in BAPS v.6 [86]. The clustering was implemented, running five independent replicates with linked loci, assuming a maximum number of groups (K = from 2 to 20) without specifying the sampling populations to geography.

Gene Flow
Migration rates between islands were estimated following a Bayesian approach implemented in Migrate v.4.4.3 [87]. The mutation scaled effective population size (Θ = xNeµ) and the estimate of migration per generation (M= m/µ) were the demographic parameters reported; both assume populations in migration-drift equilibrium [87]. This method (based on coalescence theory and a Monte Carlo Chains algorithm (MCMC)) considers the genealogical relationships of the samples to estimate demographic parameters, and provides reliable estimates of asymmetric gene flow [87]. The analysis was done using the Jukes-Cantor model, a constant mutation rate, and one long-chain of 10,000,000 steps (with three simultaneous replicates and 10% of the steps discarded as burn-in). To increase searching effectiveness, a static four-chain heating scheme was used with temperatures set to 1.0, 1.5, 3.0, and 1,000,000. Results were examined from independent runs, and the basic model of migration (migration free to vary between all the populations) was selected as the best model found having the highest probability, and each execution was accepted only when the estimates of effective sample sizes (ESS) were >1000. Pairwise comparisons among all locations could be problematic due to overparameterization and can inflate the effective population size estimates. Therefore, avoiding the violation of assumptions and to get more realistic gene flow estimates, we only made comparisons among islands considering each one as a population [87].

Correlation of Genetic Differentiation with Geographic Distances
To explore potential geographical isolation among populations, we used a Mantel test [88]. Under the assumption of isolation by distance (IBD), we expect genetic distance to increase with geographic distance. The Mantel test considers an island model, which assumes that the closer the individuals, the closer their relationship and, therefore, genetically more dissimilar from those further away [89]. Thus, the correlation between the matrix of linear geographic distances (natural logarithm of linear distances in kilometers) and the matrix of genetic distances (ΦST/1-ΦST [83,90]) was used according to Rousset [89]. Each location in the three islands was considered as a different population, and 1000 permutations were used to assess its significance.

Philopatry and Dispersal Behavior
A total of 998 Heermann's Gulls (714 males and 284 females) were recaptured in 2721 occasions (males = 1917 recaptures, and females = 804 recaptures) breeding from one to fourteen times as adults during their lifetime at seven monitored valleys in Rasa Island (Table 1). In total, 68.34% of the individuals were faithful to their natal valley for breeding (n = 682 gulls), often returning to the same nesting site or only a few meters away (namely, natal philopatry [3]). The remaining 31.66% (n = 316 gulls) dispersed to breed in other areas on Rasa Island and repeatedly returned to that same valley for breeding (namely, breeding philopatry [3]).
In all three models (Table 3), the best predictor of the location of the breeding site was the distance to the natal site with a negative coefficient, i.e., nesting birds tend to nest in the same valley where they were born, or in the close neighboring valleys (r 2 = 0.83 for Model I; 0.76 for Model II; and 0.72 for Model III; p < 0.0001 in all cases). Expectedly, the resampling effort had a small, but significant effect in all models (p ≤ 0.001 in all cases). Nonetheless, the interaction between sex and distance was also highly significant in all models, with a more pronounced negative slope for males than for females, indicating that males tend to be more philopatric. A significant interaction term between the natal and breeding valley suggests that some valleys tend to expel their offspring onto other valleys when they become adults, possibly due to differences in habitat quality. For example, Casitas Viejas Arriba is a significant (p < 0.001) expeller of offspring, many of which tend to choose Tapete Verde and Valle de la W as their preferred alternative for breeding as adults (see Table 1). Similarly, Valle de la W showed a very significant (p < 0.001) interaction term with itself, indicating that philopatry is higher here than in other natal sites, i.e., chicks raised in the valley tend to come back as adults to the same place avoiding other alternative sites. Finally, Model III, which assumed that banding at birth was done on a 50% ratio of males to females, indicated a very significant fixed effect of sex, suggesting that, if indeed sex ratio at banding was 50:50, then males are way more philopatric than females. Overall, the proportion of nesting adults that dispersed showed an inverse relationship with dispersal distance regardless of the model used, especially in males compared to females when assumed natal sex ratio is 1:1 (Model III) ( Figure 2).

Genetic Diversity
Forty-two Cytochrome b haplotypes were observed in 296 Heermann's Gulls from the Midriff Islands Region and the Mexican Province. Thirty-five haplotypes belonging to the Midriff Islands Region showed several segregating sites (S = 41), few nucleotide differences (k = 0.85), high haplotype (Hd = 0.61), and low nucleotide diversities (π = 0.00083). Conversely, only eight haplotypes were recovered from the Mexican Province, and showed fewer segregating sites (S = 12); nonetheless, there are higher values of nucleotide differences (k = 4.133), haplotype diversity (Hd = 0.93), and nucleotide diversity (π = 0.00401) ( Table 2). In Cardonosa Island, we found only three haplotypes, two segregating sites (S = 2), and the lowest values of nucleotide differences (k = 0.54), haplotype diversity (Hd = 0.51), and nucleotide diversity (π = 0.00053) ( Table 2). In Rasa Island, we found most of the haplotypes (h = 34), segregating sites (S = 40), and nucleotide differences (k = 0.88) ( Table 2). In Isabel Island, we recovered eight haplotypes and higher haplotype and nucleotide diversity values than in Cardonosa and Rasa Island. Despite sharing some haplotypes, we observed unique haplotypes on each island occurring in more than half of the locations ( Table 2).
The minimum spanning network indicates that only the haplotype H1 is shared among Cardonosa, Rasa, and Isabel Island (Figure 3), whereas the haplotype H2 was observed in both Cardonosa and Rasa Island, but not in Isabel Island. No other haplotype was observed in more than one island. Several low-frequency haplotypes with one to three mutational steps are derived from the two most common haplotypes (H1 and H2), and several low-frequency haplotypes with one to three mutational steps are derived. The occurrence of haplotype H1 in all three islands suggests that it is either very old and ancestral or a signal of genetic exchange even over long distances, or both. However, unique haplotypes suggest slight differentiation between islands, and, coupled with the resulting radial topology of relationships, there are signals of lineage sorting after a recent population expansion.

Genetic Differentiation
Comparisons between populations from different regions in all the tests revealed significant differentiation (Table 4). However, some tests revealed no significant differentiation regarding the lower geographic level. The Chi-square test (χ 2 ) showing large values suggested that the genetic differences among populations are better explained by selection or drift. Heterogeneity of haplotype frequencies among populations was detected at region, island, and cluster levels. Nonetheless, the Chi-square test is less robust if there are too many unique haplotypes or if the number of each haplotype in each location is small [74,76,77]. In contrast, the HST statistic showed low values with statistical significance (p < 0.01), indicating considerable diversity within each sample analyzed, and slight movement of haplotypes [76,77,86]. The HST statistic only showed significant genetic differentiation between regions. This could be explained by rare haplotypes grouped together and reflects the effect of population structure on the haplotype diversity [76]. Similarly, the values of KST* near to zero indicated no differentiation; however, they were statistically significant (p < 0.01), supporting genetic differentiation [91,92]. The sequencebased statistic KST indicated differentiation at all the geographic levels. The Z* statistic was also statistically significant (p < 0.01), so the null hypothesis (no differentiation) was rejected [92]. This statistic revealed significant differences between regions and islands; however, it detected no considerable differentiation among clusters and locations. The sequence-based statistic Snn revealed robust differences among regions and islands but a high genetic exchange among clusters and locations. Overall, the higher the geographic level, the more significant and stronger the differentiation observed in the Heermann's Gull subpopulations analyzed. Table 4. Genetic differentiation and gene flow estimates for the Heermann's Gull at different geographic levels (regions, islands, clusters, and locations). Genetic differentiation was estimated by haplotype (χ 2 and HST), and sequence-based statistics (KST*, Z*, and Snn), whereas the gene flow was calculated based on haplotype data (GST and Nm 1 ), nucleotide data based on diversity (NST and Nm 2 ), and nucleotide data based on global allele frequencies (FST and Nm 3 ). The GST, NST, and FST statistics showed the strongest genetic differentiation at the regions and islands levels; however, it was lower at the finest geographic levels (clusters and locations) ( Table 4). Conversely, the number of migrants (Nm) increased from higher to lower geographical level, a reflex of retention of homogeneity in the populations by lower exchange of individuals among regions and stronger exchange between locations. A widely accepted rule is that one migrant per generation will maintain a certain level of connectivity among subpopulations [93]. However, one to ten could be an appropriate number to minimize the loss of diversity within subpopulations [94,95]. In this way, despite the number of migrants per generation observed by nucleotide data (Nm 2 = 1.71 to 7.73, and Nm 3 = 1.71 to 7.75), it was not sufficient to avoid genetic differentiation among regions based on their nucleotide diversity (NST = 0.227) or random global departures of allele frequencies (FST = 0.226). In contrast, higher migration rates (Nm 1 = 17 to 145) showed lower genetic differentiation (GST < 0.028).
Paired comparisons of genetic differentiation between locations revealed significant differences in few instances, mostly between both Rasa and Cardonosa Islands and Isabel Island (Table 5). Global FST or ΦST values in seabirds are usually lower than the maximal theoretical (FST or ΦST = 1) [5,75], and we interpreted their values as low (0 to 0.05), moderate (0.05 to 0.15), high (0.15 to 0.25), and very high (>0.25) [96]. Thus, the lower the value, the higher the gene flow. Low or negative ΦST values were interpreted as low genetic differentiation and high genetic exchange between populations [83,90], i.e., the ΦST comparisons were very high (ranging 0.21 to 0.32) for the most distant population, Costa Fragata, versus all other locations, although not all were statistically significant. Comparisons between locations from Cardonosa and Rasa are predominantly undifferentiated, mainly with negative or zero ΦST values. However, low to moderate genetic differentiation within Rasa Island was observed between Valle de la W versus Vallecito del Estero, Valle Norte, Valle del Fin, Casitas Viejas Abajo, and Casitas Viejas Rocas, with significant ΦST values ranging from 0.09 to 0.22. Only Valle de la W from Rasa showed significant differentiation versus Cardonosa. Overall, more prominent differentiation between some locations than others may be due to differences in migration between breeding places, as revealed by banding on Rasa Island (see Table 1).

Genetic Structure.
The results of the AMOVA analysis are summarized in Table 6. Different patterns of structuring were found depending on groupings considered. For example, the inferred genetic sub-structuring among regions explained 48.57% of the variation, whereas comparisons among islands and clusters were 24.16% and 6.95%, respectively. However, the genetic differentiation was very high and significant among regions (ΦCT = 0.49, p = 0.05), although it was high, but not significant, between islands (ΦCT = 0.24, p = 0.10), and moderate, but not significant, among analysis clusters (ΦCT = 0.07, p = 0.12). Furthermore, from very high to moderate genetic differentiation was detected within populations (locations) and significant in all the groupings (regions, ΦST = 0.48; p = 0.0001; islands, ΦST = 0.23, p = 0.0001; and analysis clusters, ΦST = 0.05, p = 0.0001). Table 6. Analysis of molecular variance (AMOVA) of Cyt-b from the Heermann's Gull, including groupings (based on regions, islands, and clusters). For each case, it is showed the source of variation, degrees of freedom (d.f.), the sum of squares (S.S.), percentage of variation (%V), and Φ-statistics. Significant values are shown in bold (* p < 0.05, and *** p < 0.0001).

Grouping
Source of Variation d.f. The SAMOVA analysis revealed two maximally differentiated groups from each other when considering the twenty Hermann's Gull locations (Table S4). With this approach, we identified that location 1 (Costa Fragata, on Isabel Island), was grouped alone and against another group formed by all the other locations (K = 2, FCT = 0.486, p < 0.05). Similarly, the next maximally differentiated grouping was formed by the locations 7 (Valle de la W, Rasa Island), 20 (Costa Fragata, Isabel Island), and all other locations grouped (K = 3, FCT = 0.317, p < 0.05). All other K groupings showed lower, but significant, FCT values ranging from 0.105 to 0.177 (Table S4). Likewise, the Bayesian approach to estimate genetic structure implemented by BAPS also recovered two genetic groups (K = 2, log ML = −824.58, and KL divergence = 0.014). The larger one encompassed 283 individuals, mainly from gulls nesting on Cardonosa and Rasa Island (Midriff Islands Region), while the smaller one was formed by only 13 individuals, mostly from Isabel Island (Mexican Province region) (Figure 4a). Altogether, the three approaches are consistent in revealing strong genetic structure among Heermann's Gull breeding colonies nesting separately between the Midriff Islands Region and the Mexican Province biogeographic regions. Such genetic structure may result from differences at local habitat, disrupting the gene flow among Heermann's Gull subpopulations from both oceanic regimes.

Isolation by Distance
The correlation between pairwise genetic distances (ΦST/1-ΦST) and geographic distances was substantial and statistically significant (rM = 0.82, p = 0.03) (Figure 4c). Thus, an IBD model (isolation by distance), in addition to biogeographic barriers, could explain the mitochondrial genetic structure found in all the previous analyzes. Distance among locations ranged from a few meters to more than 1000 km (Table 5).

Discussion
This study provides the first evidence of population genetic structure among Heermann's Gulls breeding on islands of different biogeographic regions [48]. As in other studies of seabirds that breed in two or more ocean regimes [5,19,21,[26][27][28][29], our results support strong philopatry, local habitat, and distance as factors that reduce the gene flow among regions.

Philopatry and Dispersal Behavior
The results of our analysis confirm that Heermann's Gulls are strongly philopatric tending to return as nesting adults not only to the island where they were born but specifically to the island valley where they were hatched and reared, or at a very close distance to it ( Figure 2). Interestingly, the best predictor of the adult breeding ground was not the natal valley as a discrete factor but, rather, the quantitative distance to the place of birth (Table 3). A high proportion of the Heermann's Gull individuals were recaptured within 200 m of their hatching site, while a few moved greater distances elsewhere, as in the case of the Black-legged Kittiwake (Rissa tridactyla) [10]. This pattern may not be general in all species of birds, but there seems to be a strong tendency for philopatry among seabirds [5,75].
A second extremely consistent result, under any model hypothesis, is that males are more strongly philopatric than females (Table 3). Therefore, our study system is consistent with findings in other seabird species [1,10]. In seabirds, generally, the males will display higher philopatry than females, while the females will tend to nest farther away than males from their natal site [1,14]. Description of other patterns, however, have been found in land birds, such as the Brown Jay (Psilorhinus moio), in which philopatry is biased to females and dispersal to males [16]. Nonetheless, any interpretation should be taken with caution because any bias in the production of either sex or their survival (for which we have no control) may contribute to differences in dispersal among sexes. Such biases are not uncommon, as sex-biased production and mortality have been reported in other gulls [97][98][99]. For example, in the Black-headed Gull (Larus ridibundus), male-biased mortality has been reported: as males are larger, one possibility is that they are more prone to not fulfill their resource requirements [97]. Moreover, there is evidence that dimorphic birds overproduce the smaller sex under adverse conditions (usual females in Laridae) [61,98]. If so, females would be captured more than males because of this overproduction [98], which contrasts with our results.
We do not know the sex of the fledglings we banded initially. For this reason, we built three different models with varying fledgling sex ratio. Later, we recaptured more adult males than females, even though this is a slightly dimorphic species, with males slightly heavier and larger than females [51,61]. This can be unexpected because, as in the Black-headed Gull, males could have a higher mortality rate and, therefore, should be recaptured less [97]. However, selective factors in one species can be different from those of another, and effects may well vary [38]. Each species is in an environment with its own characteristics and natural selection pressures; in addition, each species may differ in the ways of reacting to these pressures, so we cannot expect to see the same results in two distinct species; on the contrary, what is interesting is to see the differences between them. A possible explanation for what we found is given by the Trivers-Willard hypothesis, which assumes differences in fitness returns with higher investment in the sex that gains more benefits [99]. In this way, an important part of the development of the natal population can be the more philopatric gender [100]. Then, as in the Banter See Common Stern (Sterna hirundo), male-biased natal recruitment in the Heermann's Gull suggests a greater influence of males in the natal population growth than their female equivalents [100]. The above suggests that more attention should be paid to gender bias in future studies of seabirds, as often only females are considered in population modeling and, in fact, males are left out.
Differences in the number of recaptured gulls with philopatry and dispersal behavior in each valley (Table 1) are probably related to the quality of local habitat and interspecific interactions indirectly affecting the survival [38,39,101]. If some local habitats are more suited than others (i.e., with higher quality), then some valleys will show a higher proportion of gulls recaptured with natal philopatry returning every year (philopatry benefits) [38,39]. Conversely, dispersers (i.e., gulls coming from any other nesting area) should be captured less frequently (selection against immigrants) [42]. For instance, a recent study found evidence of selection against immigrants, which were found to have a lower fitness than residents in wild seabird populations of three long-lived Procellariiformes species (Wandering Albatross Diomedea exulans, Southern Fulmar Fulmarus glacialoides, and Snow Petrel Pagodroma nivea) [42]. Other interactions may include the deterring by successive occupation of other birds arriving in large numbers or predation ( Figure S2), in addition to differences in egg sizes influenced by nest densities, habitat, and nest-site selection, as have been the cases for the Herring Gull (Larus argentatus) and the Common Gull (L. canus) [101,102]. Therefore, the nest-site selection must be a crucial decision because breeding synchronously and densely in certain areas might help to enhance offspring development and an effective defense against predation [30][31][32][33][34][35][36][37][38]44]. In the Heermann's Gull, as in many other species, producing offspring in the same environment experienced by the breeding adults at an early stage in life would be advantageous and favor the operation of group selection and local adaptation [26,[30][31][32]103].
In conclusion, both natal and breeding place appear to play an important role with possible benefits for philopatric individuals, but more research is needed to corroborate this. At the end, we do not rule out any potential bias due to different return rates by mortality, our ability to recapture the birds, or even some bias in the behavior of birds with greater territoriality that allows them to be captured more easily. The latter is matter of future research.

Gene Diversity, Genetic Structure, and Gene Flow
In Heermann's Gull, high haplotype diversity and low nucleotide diversity were observed in cytochrome b (Hd = 0.62 and π = 0.00098; Table 2). Both values of diversity indices are comparable, for instance, to mitochondrial ND2 and CR1 regions in the European Shag (Phalacrocorax aristotelis) (Hd = 0.95 and π = 0.006) [104], and to mitochondrial COI sequences in both the Black Storm-Petrel (Hd = 0.70 and π = 0.00291), and the Least-Storm Petrel (Hd = 0.95 and π = 0.00559) [25]. In seabirds with demographic expansions across a wide range, high haplotype diversity and low nucleotide diversities are probably the consequence of Pleistocene effects [49,105,106]. In our study system, unique haplotypes were found within each region, island, cluster, and even locations (Table 2, Figure 3); however, its distribution accounts for a small fraction of lineage sorting between the biogeographic regions Midriff Islands Region and Mexican Province. Overall, the diversity of cytochrome b haplotypes in Heermann's Gull agrees with demographic expansion [49], and it could be explained by regional genetic differences [106].
Given the large distances between Rasa and Cardonosa Island versus Isabel Island (>1000 km), it was natural to expect more robust differentiation if distance plays a major role as a factor of genetic structuring. Conversely, if the migratory behavior of Heermann's Gull individuals allows them to cover long distances during their lifespans (and there are no restrictions as to what places they can choose to nest and reproduce), then the genetic differentiation would be overturned by panmixia. Our results are consistent with the first scenario, as pairwise comparisons of ΦST showed very high values when comparing Heermann's Gull sequences from locations on Rasa and Cardonosa against Isabel Island (ΦST = 0.26 to 0.32; Table 5). However, although the ΦST values were usually low when locations within Rasa and Cardonosa are compared (ΦST = 0), there were instances of very high values (ΦST = 0.22), indicating that not only geographic distance but other factors (e.g., behavior) are also playing a role in genetic structuring at fine scales. Our results are comparable to pairwise ΦST differences among Atlantic Kittiwake (Rissa t. tridactyla) colonies (ΦST = 0.10 to 0.40) [107] and between the Red-billed Tropicbirds (Phaethon aethereus) populations on islands from the Mexican Ocean Pacific (ΦST = 0.0 to 0.20) [22].
On the other hand, very high genetic structure was observed between Heermann's Gull colonies from the Midriff Islands Region and the Mexican Province region. This was consistent through all three approaches implemented (AMOVA, ΦCT = 0.49, p < 0.05; SAM-OVA, FCT = 0.486, p < 0.05; and BAPS, log ML = -824.58, and KL divergence= 0.014). This pattern of population genetic structure may be the result of differences between ocean regimes and intrinsic ecological barriers. For example, this was the case between Atlantic and Pacific populations of the Black-legged Kittiwakes (Rissa tridactyla) (AMOVA, global ΦST = 0.53, p < 0.05) [107], and between the Atlantic and Indian populations of the Thinbilled Prion (Pachyptila belcheri) (AMOVA, FST = 0.35) [21]. However, this is not always the case. For instance, low genetic structuring among groups was found by Rocha-Olivares et al. [19], between colonies of Mexican Frigatebirds (Fregata magnificens) from the Mexican Ocean Pacific and the Mexican Caribbean (AMOVA, ΦCT = 0.015 p = 0.05).
Although, in the Heermann's Gull, the number of migrants per generation was more than one (Nm = 2 to 8), this was insufficient to create a single panmictic and homogenized population according to the genetic structure observed (Tables 6 and S3, Figure 4a). Moreover, the asymmetric gene flow observed between islands where Heermann's Gull breeds showed a southward decrease in the exchange of individuals ( Figure 4b). As in many other threatened seabirds, Heermann's Gull are usually found in fragmented and partially isolated populations, in which cryptic speciation may be found [6,[27][28][29]104]. Therefore, although the one migrant per generation rule simplifies assumptions about the connectivity and the loss of polymorphism among subpopulations, several factors may influence the genetic exchange between populations [5,[93][94][95]. In this last sense, our genetic analyses showed that physical distance is an essential factor shaping genetic structure among Heermann's Gull colonies breeding on islands of different ocean regimes (Figure 4c). Nonetheless, if the breeding success and reproductive decisions are occurring influenced by predation and multigenerational social interactions incentivized by inclusive fitness, then behavior may also play an important role in the formation of genetic structure in the Heermann's Gull as well as in other species [15,32,103]. Altogether, our results agree with the rationale of Friesen [5], as geographically separated populations (both by physical and non-physical barriers) tend to be less connected by gene flow, preventing dispersal and leading to very high genetic structure.

Conclusions
In the Heermann's Gull, our banding data showed strong and male-biased natal philopatry in Rasa Island, while the genetic data revealed consistent genetic differentiation between colonies breeding in the Midriff Islands Region and the Mexican Province. Furthermore, despite their high dispersal potential, the Heermann's Gull displayed uneven gene flow and positive isolation by distance pattern between populations. The strong genetic structure observed among regions is supported by unique haplotypes in each subpopulation and may reflect philopatry as a crucial behavioral barrier reinforced by physical isolation (distance) and different local habitat conditions. Supplementary Materials: The following supporting information can be downloaded at: www.mdpi.com/article/10.3390/d14020108/s1; Figure S1.  Table S2. Letters at the tip of each label indicates the geographical origin: Grays Harbor (_Grays), Cardonosa (_Car), Rasa (_Ras), and Isabel (_Isa). The complete cytochrome b from Rissa tridactyla was used as the outgroup (GenBank: DQ385229). The HKY+G model was selected as the best substitution model. For each branch, the aLTRT support is shown. The most likely tree is shown, and their log likelihood tree was -2287.969; Table S1. Number of Heermann's Gull fledglings ringed in nine Rasa Island valleys between 1984 and 1993; Table S2. Comparison of sequences reported in the literature for the Heermann's Gull cytochrome b; Table S3. Accession numbers for each haplotype of the Heermann's Gull in the Gen Bank; Table S4. Fixation indexes for groups of populations maximally differentiated (K) and inferred by the Spatial Analysis of Molecular Variance (SAMOVA). The inferred groups based on the Heermann's Gull cytochrome b sequences from each location are described by the square brackets (Numbers between brackets correspond to the locations in Figure 1

Data Availability Statement:
The data presented in this study are available in Supplementary Materials. Also, can be found at the GenBank (https://www.ncbi.nlm.nih.gov/, accessed on 25 January 2022) or requested from the corresponding authors.