Strong phenotypic divergence in spite of low genetic structure in the endemic Mangrove Warbler subspecies (Setophaga petechia xanthotera) of Costa Rica

Abstract Despite the enormous advances in genetics, links between phenotypes and genotypes have been made for only a few nonmodel organisms. However, such links can be essential to understand mechanisms of ecological speciation. The Costa Rican endemic Mangrove Warbler subspecies provides an excellent subject to study differentiation with gene flow, as it is distributed along a strong precipitation gradient on the Pacific coast with no strong geographic barriers to isolate populations. Mangrove Warbler populations could be subject to divergent selection driven by precipitation, which influences soil salinity levels, which in turn influences forest structure and food resources. We used single nucleotide polymorphisms (SNPs) and morphological traits to examine the balance between neutral genetic and phenotypic divergence to determine whether selection has acted on traits and genes with functions related to specific environmental variables. We present evidence showing: (a) associations between environmental variables and SNPs, identifying candidate genes related to bill morphology (BMP) and osmoregulation, (b) absence of population genetic structure in neutrally evolving markers, (c) divergence in bill size across the precipitation gradient, and (d) strong phenotypic differentiation (P ST) which largely exceeds neutral genetic differentiation (F ST) in bill size. Our results indicate an important role for salinity, forest structure, and resource availability in maintaining phenotypic divergence of Mangrove Warblers through natural selection. Our findings add to the growing body of literature identifying the processes involved in phenotypic differentiation along environmental gradients in the face of gene flow.


| INTRODUC TI ON
Populations distributed along environmental gradients provide excellent systems to study the counteracting effects of natural selection and gene flow on population trait divergence at relatively small geographic scales (Bertrand et al., 2016;Doebeli & Dieckmann, 2003;Milá, Warren, Heeb, & Thébaud, 2010;Postma & Noordwijk, 2005;Richardson, Urban, Bolnick, & Skelly, 2014;Schluter, 2009;Seeholzer & Brumfield, 2018;Smith et al., 2005). If environmentally mediated divergent selection is strong enough to counteract the homogenizing effects of gene flow, trait variability could persist through time (Milá et al., 2010;Richardson et al., 2014;Schluter, 2009;Smith et al., 2005). This ecologically driven trait divergence could result in speciation if traits subject to natural selection secondarily influence sexual traits (Schluter, 2001). Since in this model, there is no physical separation between incipient species, there is still potential for interbreeding and gene flow during the speciation process (Bolnick & Fitzpatrick, 2007;Foote, 2018). For these reasons, it is essential to understand the conditions and processes under natural selection that can drive ecological divergence and reproductive isolation in continuously distributed populations, especially during early stages of formation of biological diversity (Bolnick & Fitzpatrick, 2007;Cheviron, Connaty, McClelland, & Storz, 2014).
The Pacific coast of Costa Rica presents a strong precipitation and salinity gradient in which yearly rainfall varies between F I G U R E 1 (a) Map showing the sampling locations of Mangrove Warbler populations and the precipitation gradient. The salinity gradient correlates with the precipitation gradient such that drier sites have higher salinity levels. Each individual sampling site is noted by different colors shown in the map. Population names correspond to 1 = Naranjo, 2 = Junquillal, 3 = Chira, 5 = Chomes, 6 = Tarcoles, 7 = Sierpe, 8 = Osa, 9 = Golfito. (b) Principal component analysis showing the morphological distribution of individuals along precipitation gradient in Costa Rica. Colors show average precipitation gradient and dashed lines with black numbers represent the average salinity levels 1,000 mm in the north and ~6,000 mm in the south creating a salinity gradient that has a strong influence on mangrove forest structure ( Figure 1a). Where rainfall is low, and salinity is high, canopy level rarely exceeds 20 m in height, while in areas with high rainfall and low salinity canopy can exceed 35 m (Jiménez, 1990). Several bird species are distributed along the gradient from which Mangrove Warbler (Setophaga petechia xanthotera) is restricted to the mangrove habitat. The environmental gradient of the Pacific coast of Costa Rica, then, has the potential to influence the divergence of traits involved in physiology, foraging, and possibly behavior of this insectivorous habitat specialist bird.
For example, insects accumulate excess salt in their exoskeleton in high salinity environments (Bradley, 2008). Salt regulation is particularly problematic for passerine species because they lack salt glands. Thus, differences in salinity and water availability should promote divergence of osmoregulatory genes or their expression, to help individuals deal with salt regulation and water loss at environmental extremes (Sabat, Maldonado, Rivera-Hutinel, & Farfan, 2004;Sugiura, Aste, Fujii, Shimada, & Saito, 2008). Additionally, differences in forest structure driven by the abiotic environment could influence size distribution of insect prey (Janzen & Schoener, 1968). Other characteristics of the forest such as understory density and overall forest interior structure are also affected by precipitation and salinity. Specifically, bill morphology might change along the gradient as a response to the change in overall resource size distribution (Grant & Grant, 2011), while wing, tail, and tarsus morphology could be influenced by the forest structure as these traits are directly related to flight performance and maneuverability (Milá et al., 2010;Pennycuick, 1968;Ricklefs, 2012;Thomas, 1996). Consequently, abiotic environment, resources, and forest structure can all have an impact on ecophysiological and morphological traits of Mangrove Warbler.
Prior to this study, we posited that Mangrove Warbler populations would have high levels of gene flow based on two facts: (a) There are no strong geographic barriers (e.g., mountains ranges, large river basins) to isolate the Mangrove Warbler populations along the gradient and (b) high levels of gene flow have been previously documented among populations of the Galapagos and Coco's island warbler (Setophaga petechia aeurolea) in spite of their strong geographic isolation (Chaves, Parker, & Smith, 2012). For these reasons, Mangrove Warblers on the Pacific coast of Costa Rica should provide an excellent system to study the effect of the environment on trait divergence with potential gene flow.
Studies that investigate the role of environmental gradients on adaptive variation in birds largely focus on divergence along elevational gradients (Cheviron & Brumfield, 2009;Milá, Wayne, Fitze, & Smith, 2009;Schluter, 2001;Seeholzer & Brumfield, 2018;Vines & Schluter, 2006). Few studies have attempted to determine the role of precipitation and salinity gradients in birds phenotypic variability (Bay et al., 2018), and even fewer studies have focused on understanding the genomic patterns behind trait variation along environmental gradients (Cheviron et al., 2014). In this study, we obtained samples from nine populations of Mangrove Warbler distributed along a steep precipitation gradient on the Pacific coast of Costa Rica to explore the role of the environmental gradient on phenotypic and genetic divergence. Our main objectives were to (a) identify genomic signals of selection, (b) determine levels of neutral genetic divergence among these populations (F ST ), (c) estimate levels of phenotypic divergence along the environmental gradient (phenotypic structure: P ST ), and (d) determine the balance between genetic and phenotypic divergence.

| Mangrove Warbler sampling
Between 2013 and 2015, we visited nine localities distributed along the precipitation gradient on the Pacific coast of Costa Rica (Figure 1a).
At these localities, we captured 115 adult male and female Mangrove Warblers (S. p. xanthotera). For every individual captured, we measured body mass, bill length (from the nares to the tip) and height (top and bottom bill at the nostrils point), tarsus (from the inner bend of the tibio-tarsal articulation to the base of the toes), flattened wing cord length (bend of the wing to the tip of the longest primary feathers), and tail length (base of the tail to the tip of the longest feathers). These morphological traits were measured with a caliper (±0.005) except for body weight which was measured with a 100 g (±0.01) analog scale. These traits were chosen because they are expected to respond to ecological differences in habitat (Eroukhmanoff et al., 2013;Grant & Weiner, 1999;Maley, 2012;Ricklefs, 2012). We also obtained blood samples for genetic analyses, which were stored in lysis buffer (0.1 M Tris-HCl, 0.1 M EDTA, 0.01 M NaCl, 2% SDS).
We mapped all mangrove localities reported for Costa Rica using recent satellite images available at Centro Nacional de Alta Tecnología (CENAT). The nine populations chosen along the environmental gradient were selected based on their accessibility.
Although we considered including Atlantic populations in our sampling, we were not able to find any Mangrove Warbler individual at the Caribbean during the time of our study. Using digitized maps, we calculated the Euclidean and coastal distances among sampling populations. Euclidean distance refers to the straight line between sites (so might involve flying over water) while coastal distance is the distance of shoreline between two sites along the coast (assuming birds would not fly over water). Sampling sites were between 15 and 345 km apart using Euclidean distance and between 15 and 2,584 km apart using coastal distance.

| Environmental variables and habitat classification
To characterize the environment in each of the populations sampled, we used georeferenced environmental data sets from the Meteorological Institute of Costa Rica, WorldClim (Fick & Hijmans, 2017), and a database from the University of Costa Rica Marine Investigation Center (CIMAR). We gathered data for a set of eight environmental variables: mean annual temperature (BIO1), mean diurnal temperature range (BIO2), temperature seasonality (BIO4), annual precipitation (BIO12), precipitation of the wettest month (BIO13), precipitation of the driest month (BIO14), precipitation seasonality (BIO15), elevation (SRTM), and mean water salinity (ppt) (CIMAR). However, since precipitation variables (BIO12, BIO13, BIO14, BIO15) were highly correlated (R > 0.88 among the four variables), and we considered that precipitation was the principal factor affecting conditions along the gradient (since precipitation affects soils salinity levels which in turn influence the forest structure), consequently, we used only annual precipitation (BIO12) to describe the environmental gradient and for analysis of our morphological and genotypic data (Table 1).

| Collection of genetic data and identification of SNPs
To estimate the variation in genetic structure and patterns of gene flow among populations of Mangrove Warbler distributed along the environmental gradient, we obtained SNPs using double-digest RADseq (ddRADseq). We first extracted DNA from the blood samples using the Qiagen PureGene DNA Isolation kit. DNA extractions were quantified using a NanoDrop ND8000 spectrophotometer and a Qubit 2.0 fluorometer with the DNA HS assay kit (Life Technologies) and checked for DNA quality on an agarose gel to select samples with appropriate DNA concentration (>20 ng/ µl) and quality. From the 115 individuals collected and for which we had morphological data, we only succeeded in extracting DNA, constructing libraries, and obtaining reliable genetic data from 68 of them (Table 1), following the method proposed by Peterson, Weber, Kay, Fisher, and Hoekstra (2012), using EcoRI and MspI.
Each individual RAD library was ligated to a unique molecular identifier using one of four types of DNA barcodes (either 8, 9, 10, 14 bp in length). Individuals were pooled together, fragments were subject to size selection to produce a mean fragment length of 250-440 bp, and then sequenced on a single Illumina NextSeq 500 run. The sequenced reads were quality-filtered, and individual barcode information was removed using the process_radtags program in STACKS (1.46) (Catchen, Amores, Hohenlohe, Cresko, & Postlethwait, 2011;Catchen, Hohenlohe, Bassham, Amores, & Cresko, 2013). Afterward, PCR clone sequences were eliminated with "clonefilter" in STACKS v1.30 .
We chose STACKS because it is a specialized software pipeline for building loci from short-read sequences with restriction enzyme-based data. We aligned the samples based on the reference genome available for Yellow Warbler (Setophaga petechia petechia) (Bay et al., 2018) using BWA-MEM V.0.7.9a (Li & Durbin, 2009).
This allowed for additional positioning information and facilitated detecting rare allele variants (Peterson et al., 2012) that are often removed from de novo assemblies, as they can be confounded with sequencing errors.
Subsequently, we used PSTACKS, CSTACKS, and SSTACKS (Catchen et al., , 2013Paris, Stevens, & Catchen, 2017) to identify SNPs. Parameters for these alignments included a terminal threshold of 500, a maximum number of mismatches allowed (M = 5), a minimum stack depth of three (m = 3) among reads with potentially variable sequences (Paris et al., 2017), and an indel penalty of 2 . In the population's module of STACKS and following consecutive filtering steps, we first retained RAD tags with a minimum stack depth (m) of 20 and a maximum stack depth of 100.
This step removed SNPs genotyped with too low coverage (m < 20) to be accurately called as well as SNPs genotyped with too high coverage (m > 100) that might reflect repetitive regions. Then, we retained SNPs genotyped in at least 80% of the individuals and 80% of the sampling locations and excluded markers showing heterozygosity >0.50 within samples (Hohenlohe, Amish, Catchen, Allendorf, & Luikart, 2011). We also removed markers out of Hardy-Weinberg equilibrium (p-value = .01) at more than 60% of the locations. SNPs with more than 30% missing data were also eliminated. Finally, we removed SNPs with very low frequency (MAF < 0.05), as these can create biases in quantifying genetic connectivity and should be removed when inferring demographic processes (Roesti, Hendry, Salzburger, & Berner, 2012). The number of SNPs kept after each filtering step can be found at the Appendix S1 in Table S1. Finally, we used the "Populations" module in STACKS to obtain individual genotypes per populations.
TA B L E 1 Coordinates, habitat, and a mean annual precipitation (BIO12), mean salinity levels (ppm), canopy height (m) of localities sampled and number of individuals by population, used in the genetic (SNPs) and morphological analyses

| Identification of outlier loci and correlation with the environment
Single nucleotide polymorphisms potentially under balancing and divergent selection should be removed when assessing genetic connectivity (gene flow) among populations (Beaumont & Nichols, 1996). We searched for loci with a level of population differentiation exceeding neutral expectations using two F ST based outlier analyses. First, we used the software OUTFLANK (Whitlock & Lotterhos, 2015), which calculates a likelihood based on a trimmed distribution of F ST values to infer the distribution of F ST for neutral markers.
Second, we detected outlier SNPs with BAYESCAN v. 2.1 (Foll & Gaggiotti, 2008) that estimates population-specific F ST coefficients using the Bayesian method and uses a cutoff based on the mode of the posterior distribution to detect SNPs under selection (Foll & Gaggiotti, 2008). SNPs with a posterior probability over 0.95 were considered as outliers, after running 100,000 iterations on all samples together (i.e., not pairwise, with remaining default parameters).
We specified a 'prior' odd of 10,000, which set the neutral model being 10,000 times more likely than the model with selection to minimize false positives (Whitlock & Lotterhos, 2015). Using the results of these two analyses, we divided our data set into two categories, neutral SNPs and SNPs under divergent selection. We used the neutral SNPs to calculate values of pairwise F ST .
To identify SNPs associated with environmental parameters, we used BayeScEnv (Villemereuil & Gaggiotti, 2015), an approach similar to BAYESCAN, but which aims to detect outlier loci associated with environmental parameters. BayeScEnv computes posterior probabilities of three models: a neutral model, a locus-specific model, and a local adaptation model linked to the environmental variable (Villemereuil & Gaggiotti, 2015). For this approach, we used all the SNPs 15,307 (neutral and under selection), and we used only mean annual precipitation (BIO12) as predictor variable. We ran BayeScEnv 10 times and averaged results over the 10 independent runs. We used the default parameters recommended for long runs to achieve convergence of MCMC (Foll & Gaggiotti, 2008;Villemereuil & Gaggiotti, 2015) and used a false discovery rate of 0.05 to reduce the number of false positives (Foll & Gaggiotti, 2008).
As an alternative approach to confirm if the outlier loci identified by BayeScEnv were consistent among methods, we used the latent factor mixed models (LFMM; Frichot, Schoville, Bouchard, & François, 2013), which measures the associations between genotype and phenotype or environmental variables while accounting for underlying population structure. We ran five separate MCMC runs with a latent factor of K = 1, based on preliminary structure and PCA, and using mean annual precipitation (BIO12). We used only average bill length for the phenotype-genotype association, as bill length and bill height were highly related (R = 0.8557, p = <2.2 e −16 ). p-Values from all five runs for the two independent genomewide association tests were combined and adjusted for multiple tests using a false discovery rate correction of 0.05. We used the default parameters recommended for long runs to achieve convergence of the MCMC (Frichot et al., 2013), as we did with BAYESCAN and BayeScEnv.
Finally, we also performed a redundancy analysis (RDA) using the entire set of loci to identify outlier loci and their correlation to environmental variables (Forester, Lasky, Wagner, & Urban, 2018).
Once obtained the first three axes from the RDA, we performed an outlier identification process by assuming that outlier loci were located above or below three SDs from the mean of the empirical distribution given by the scores of each axis. Finally, we correlated the observed allele frequency across populations with each of the four environmental variables to determine which was the environmental variable that explained the largest amount of variance in the structure of those outlier loci (Forester et al., 2018).

| Candidate gene identification
We used the outlier loci identified by all four methods (p < .08) and aligned the sequences to the reference genome of Zebra finch (Taeniopygia_guttata-3.2.4 reference Annotation Release 103 NCBI (https ://blast.ncbi.nlm.nih.gov/Blast.cgi), using the basic alignment search tool (Blast) from NCBI to align the sequences. We used "nr/ nt" database and with setting parameters max target sequences to 100, expect thresholds to 10, word size to 28 and max matches in a query to 0. We considered a locus homologous if the e-value returned was smaller than 1.0 e −10 . We determined whether there was any gene within 25 kb upstream or downstream of each candidate SNP to focus on genes likely to be within the same linkage group as our SNP (Bay et al., 2018). Then, we calculated the allele frequencies for each SNP under selection that was linked to a candidate gene previously associated with morphological, phenotypic, and metabolic functions related to environmental variables (Bay et al., 2018).
To calculate the allele frequencies, we only used the populations in which we had more than five individuals. If there was more than outlier SNP linked to the same candidate gene (e.g., BMP5), we calculated average allele frequencies across all SNPs.

| Population structure of neutral loci
We used ADMIXTURE to evaluate population structure with different numbers of hypothetical populations (k). We ran ADMIXTURE ver. 1.22 (Alexander, Novembre, & Lange, 2009) using 20,000 bootstrap replicates. We used k values between 1 and 9 with five iterations for each value; stabilization of parameters was checked for 100k length of burn-in and 100k MCMC simulations. To evaluate optimal partitioning in ADMIXTURE, cross-validation (CV) error values were computed for each k using a fivefold cross-validation procedure.
To determine how neutral genetic variation in our data set was distributed among populations, we conducted a principal component analysis (PCA) in the package ADEGENET 2.0.1 (Jombart & Ahmed, 2011), excluding the loci identified to be under selection by BAYESCAN and OUTFLANK (see above).
We identified SNPs that were under linkage disequilibrium (LD) by using the SNPRelate software assuming a threshold of 0.5 (Zheng et al., 2012(Zheng et al., , 2017. We found that 54% of our SNPs were in LD  (Goslee & Urban, 2007;Lichstein, 2007;Wang, 2013). Significance was assessed using 10,000 permutations of the distance matrices. F ST , coastal and precipitation distance matrices were scaled prior to analysis by subtracting the mean and dividing it by the standard deviation of the data set. Scaling of the predictor variables allows for direct comparison of the regression coefficients in order to understand the relative contribution of each independent variable over genetic distance.
Since Mangrove Warbler is a species restricted to Mangrove habitats, we believe that coastal distance between populations is a good proxy for dispersal distance and thus did not consider Euclidean distance in the Matrix Regression. Analyses in R were performed using version 3.5.2 (R Core Team, 2018).

| Phenotypic variability and traits under selection
To determine patterns of morphological variation along the gradient, we performed regressions for each trait against precipitation.
It is possible that changes in morphological traits are the result of allometry between body mass and other traits. To account for this, we also regressed the residuals of the relationship between body mass and each trait against precipitation. We also reduced the five morphological traits into two principal components (PC). Finally, we fitted a smooth surface using precipitation (BIO12) and salinity (ppm) as dependent variables and the morphological PC as coordinates for fitting. The regressions and principal component analysis were performed in R, and for fitting the smooth surface, we used ordisurf function in the VEGAN v.2.0-10 package (Oksanen et al., 2016).
To assess the level of phenotypic structure in our data, we compared neutral genetic differentiation (F ST ) to phenotypic differentiation (P ST ). To calculate P ST , we first estimated within-population and among-population variance using a linear mixed model with only intercept as a fixed effect and populations as random effects. We then used the within-population variance as the residual variance and between-population variance as the variance of the random effect. Confidence intervals of the within-and between-population variances were estimated using 1,000 parametric bootstrap replicates (Leinonen, McCairns, O'Hara, & Merilä, 2013). The variances and confidence intervals were estimated using the lmer and confint functions in the lme4 package (Bates, Maechler, Bolker, & Walker, 2015).
Using the estimated within-and between-population variances, we quantified the phenotypic divergence in a trait across populations using P ST for each morphological trait (except body mass; Brommer, 2011).
In this equation, 2 B and 2 W are the between-and within-population phenotypic variances, respectively, h 2 is the heritability of the trait under study, and the scalar c expresses the proportion of the between-population variance that is due to genetic effects across populations (Brommer, 2011). Under controlled conditions, phenotypic differences should be entirely due to additive genetic effects, so c/h 2 = 1 and P ST is equivalent to Q ST , and analogous to Q ST for a given quantitative trait (Wright, 1950). In wild populations, h 2 and c are usually difficult to estimate (Brommer, 2011) and nonadditive genetic effects such as selection can strongly influence the estimation of P ST (Brommer, 2011;Brommer, Hanski, Kekkonen, & Väisänen, 2014;Leinonen, Cano, Mäkinen, & Merilä, 2006;Pujol, Wilson, Ross, & Pannell, 2008). Consequently, we used a sequence of 100 values of c/h 2 between zero and two (Brommer, 2011). The objective of this approach is to estimate the value of c/h 2 for which P ST is larger than F ST . The smaller the critical value of c/h 2 for which P ST is larger than F ST , the more likely it is that selection influences morphological trait evolution.
The critical c/h 2 value thus reflects the robustness of the comparison between P ST and F ST (Brommer, 2011). As c/h 2 approaches 1, the morphological trait is assumed to evolve under neutral conditions. We repeated this procedure for each phenotypic trait measured in the field since all of them have been shown to be heritable in multiple bird species (Charmantier, Kruuk, Blondel, & Lambrechts, 2004;Teplitsky, Robinson, & Merilä, 2014). For interpretation, we use the value of c/h 2 for which P ST is larger than F ST computing P ST using the lower boundary of the confidence interval of 2 B and 2 W .

| Genomic signals of selection
After filtering, we obtained 15,307 SNPs (Table S1, Appendix S1). We removed outlier loci identified with BAYESCAN (20) and OUTFLANK (25) such that 15,262 putatively neutral loci remained to compute pairwise F ST . Using BayeScEnv and LFMM, we identified a total of 38 outlier SNPs associated with precipitation (BIO12) from the original loci set: 14 SNPs were found exclusively with the BayeScEnv approach, 11 SNPs were found exclusively with LFMM, and 12 SNPs were found with both methods (Figure 2, Table 2). In some cases, multiple SNPs mapped to the same locus. From the 38 SNPs, we identified 19 genes; these include functions such as supplying calcium to cardiac muscle (RYR2), neural regulation (NRG3), different cell processes (CCSER1, DDX10), GTPase activation (RIN3), protein kinase activation (LOC100229672), and transmembrane proteins (CCDC91). The strongest associations (lowest p value) between genotype and precipitation were upstream of genes with known function in avian morphology and osmoregulation (Figure 2).
For osmoregulation, the strongest associated genes were a sodium/ chloride exchanger (LOC100224232), potassium channel regulators (KCHN7, KCHN8), and aquaporin 1 (AQP1) (Figure 2, Table 2). In addition, we found strong associations with a candidate gene with known function in avian morphology, bone morphogenetic protein (BMP5), which plays a key role in bone and cartilage development ( Figure 2, Table 2).
We identified a total of 12 SNPs that were associated with bill size using the LFMM: two SNPs were found exclusively with this approach, while the other 10 SNPs were also identified in previous analyses with precipitation using either BayeScEnv or LFMM ( Figure 2, Table 2). From the 12 SNPs identified, seven genes had known functions including supplying calcium to cardiac muscle (RYR2), neural regulation (NRG3), cell processes (CCSER1), transmembrane proteins (CCDC91). Bone morphogenetic protein (BMP5) showed the strongest associations between genotype and phenotype, similar to what we found in our genotype-environment association (Figure 2, Table 2). We found that variation in allele frequencies followed environmental changes in precipitation (Figure 3).
Complementary, the first three axes of the redundancy analysis explained 84% of the variability in SNP loci. Although the correlation with environmental variables was not significant (p > .05), the first axis was positively related to Bio4 and the second and third axis were positively related with Bio1, Bio 12, and Bio 15 ( Figure S4). The outlier analysis using the RDA scores from the first three RDA iden-  Figure S4). In addition, we found candidate genes with known function in avian bone morphology (BMP1, COL21A1, LOC105759070, LOC100226932), which could play roles in the development of bone and cartilage (Table S2, Figure S4). TA B L E 2 Genes identified to be associated with precipitation, their scaffold, and position (in base pairs) according to the Yellow Warbler genome and the chromosome location according to the Zebra Finch genome

| Population structure
Pairwise F ST using putatively neutral loci ranged from 0.005 between two populations in the same habitat (drier end of the gradient) to 0.057 between two populations in different habitats (at the drier and wetter ends of the gradient), but none of the pairwise F ST values were significant (p > .05). The mean pairwise F ST was 0.015 ± 0.09 (Table 3). The F ST values were too low (<0.10) for reliable estimation of migration rates, which suggests ongoing gene flow among S. p. xanthotera populations (Meirmans, 2014).
Additionally, it did not observe significant population structure from ADMIXTURE analysis ( Figure S1, Appendix S1). The best fitting resolution according to the calculation of CV errors was k = 1 ( Figure S2, Appendix S1). Multiple matrix regression showed no influence of coastal or precipitation distance and genetic distances among populations suggesting little influence of isolation by distance and isolation by environment at least on neutral loci (R 2 = 0.06, F = 1.1, p = .36; coastal distance, β = 0.003, p = .5; precipitation distance, β = −0.006; p = .17).
We found high congruence between the lack of support of population structure from ADMIXTURE and the clustering pattern by PCA, which suggested no distinct population clustering (Figures S1 and S2). The populations of Mangrove Warblers did not cluster geographically or by environment according to their scores on the first two axes ( Figure S3, Appendix S1), which accounted for 4.2% and 3.8%, respectively, of the observed genetic variation. In addition, the PCA using only loci that are not in LD did not find distinct population clustering, with PC1 and PC2 accounting for 7.6% and 8.2% of genetic variation ( Figure S3). Also, the DAPC analysis did not identify any clusters within the data ( Figure S4, Appendix S1).

| Phenotypic trait differentiation and the role of environmental factors
We found that bill height, bill length, body mass, wing length, and tibia length were significantly related to mean annual precipitation (Figures   1b and 4; Table 4). Bill height and length increased with mean annual precipitation while wing length, tibia length, and body mass decreased F I G U R E 4 Relationship between morphological traits and mean annual precipitation (BIO12) in individuals of Mangrove Warbler. Solid black line represents the estimated model and grey polygon shows de 95% confidence interval of the regression. Panels with no regression lines indicate that the regression was not significant as precipitation increased (Figure 4). We found similar results when accounting for the allometric relationship between body mass and other phenotypic traits, so we present only the raw data. The principal component analysis showed that the first PC was related mostly to bill mor- was higher than F ST in two (bill height and bill length) of the six traits evaluated ( Figure 5, Table 5). The critical value of c/h 2 for bill length was 0.05, and for bill height, it was 0.13 ( Figure 5, Table 5).

| D ISCUSS I ON
Using an integrative approach, we found that Mangrove Warbler populations along the precipitation/salinity gradient of the Pacific coast of Costa Rica maintain significant phenotypic divergence despite the absence of genetic structure across most of the genome, suggesting high gene flow among populations. According to the four objectives of our study, our results show that (a) genes associated with bill growth and osmoregulatory pathways are associated with precipitation, (b) there is extremely low genetic structure in neutrally evolving loci, (c) morphological traits (bill size) change significantly along the gradient, and (d) bill phenotypic differentiation (P ST ) is substantially higher than genetic differentiation (F ST ). All these evidence points to the hypothesis that divergent natural selection at the gradient's extremes may be strong enough to counteract the homogenizing effect of gene flow potentially promoting initial steps of ecological speciation.

| Candidate genes at outlier loci
Combining all analysis used to identify outlier loci, we found 23 candidate genes related to osmoregulation processes. Some of these candidate genes are specialized to activate sodium-potassium channels which could help shed the excess inorganic ions and retain water (Maley, 2012). The candidate genes AQP1 and AQP4 code for aquaporin protein types 1 and 4. Aquaporins (AQPs) are a family of transport proteins that confer high-membrane water permeability in various tissues in animal, plants, and microorganisms. In chickens, for example, AQPs play a major role in regulating the total body water balance by concentrating or diluting uric acid (Sugiura et al., 2008).
In addition, the candidate gene NPNT is related to kidney development, which is the key organ associated with osmoregulation in passerine birds (Sabat et al., 2004). These outlier genes known to be involved in osmoregulation suggest that salt regulation and water availability may pose a challenge to Mangrove Warblers and potentially other passerine insectivorous species in the Pacific coast of Costa Rica. This physiological divergence in contrasting environments has also been found in a number of taxa, including killifish (Fuller, Mcghee, & Schrader, 2007;Whitehead & Crawford, 2006), sunflowers (Karrenberg, Edelist, Lexer, & Rieseberg, 2006), and other bird species (Maley, 2012).
In addition, the candidate genes COL21A1, LOC105759070, LOC100226932 are involved in cartilage development, which could be related to bill morphology. Comparisons of P ST /F ST indicated that variation found in bill size cannot be explained by genetic drift alone. These results along with the association between BMP and environment suggest that observed bill morphology variation in Mangrove Warbler individuals along the gradient may be the result of natural selection (Figure 3).
Although it is well known that bill size is highly heritable in birds (Eroukhmanoff et al., 2013;Grant & Weiner, 1999;Maley, 2012;Ricklefs, 2012), P ST was calculated only from phenotypic data and thus we cannot disentangle the contribution of plasticity and genetic variation in the observed trait variability. Since most documented phenotypic traits are affected by environmental conditions, at least some fraction of variation in bill size could be due to phenotypic plasticity. To measure the exact contribution of the additive genetic differentiation to bill size, it would be essential to calculate Q ST under common garden conditions (Pujol et al., 2008).
Such experiments, however, are unfeasible to perform in most birds and especially on an endemic endangered subspecies such as Mangrove Warbler.

| Phenotypic divergence in the presence of gene flow
The strong phenotypic divergence observed along the gradient contrasts with the lack of genetic structure. Similar results have been obtained in marine organisms where high dispersal distances lead to high levels of gene flow (Fuller et al., 2007;Whitehead & Crawford, 2006). Other studies in birds have also reported strong phenotypic differentiation with gene flow mainly along elevational gradients (Cheviron & Brumfield, 2009;Milá et al., 2009;Seeholzer & Brumfield, 2018 )-and within ( 2 W )-population variance components are followed in parentheses by the 95% confidence intervals from the posterior distribution of a Bayesian generalized linear mixed model. P ST values for each trait are computed at the null assumption of c/h 2 = 1 with 95% confidence intervals in parentheses. c/h 2 * is the critical value of c/h 2 at which the lower confidence interval for P ST exceeds the global F ST . Lower critical values represent more robust inferences of environmental selection. Wilson & Rannala, 2003;Migrate-n, Beerli & Palczewski, 2010).
Nonetheless, these methods are not accurate when divergence among populations is very recent and population structure is low (specifically F ST < 0.1; Faubet, Waples, & Gaggiotti, 2007;Samarasin, Shuter, Wright, & Rodd, 2017). Since Mangrove Warbler populations are hypothesized to have recently diverged from their ancestors (Chaves et al., 2012) and that we found low population structure, we believe that methods available for estimation of the number of migrants per generation in this case are not appropriate. It is known, however, that PCA of SNP data has a genealogical interpretation allowing us to make inferences about migration and connectivity among populations (McVean, 2009). Thus, from our results we can infer that migration among Mangrove Warbler populations is high because we did not observe defined clusters in our PCA based on putatively neutral loci.
It has been previously reported that bill height and length can respond to selection driven by food resource availability (Grant & Grant, 1993, 2011. Insect size distribution is positively influenced by precipitation (Janzen & Schoener, 1968) supporting our observations of smaller bill sizes in the drier habitat. The fact that bill morphology responds in the predicted direction and that we found associations between precipitation and genes in the BMP family, support the hypothesis that morphological divergence in Mangrove Warbler populations is maintained through natural selection. Alternatively, bill size has been associated with accessibility to prey, as relative longer bills length has been associated with foraging strategies that require access to prey hidden in deeper substrate (Wright & Steadman, 2012). It is possible that higher humidity in wetter forests provides habitat for a wider variety of epiphytes which in turn provide refuge for bird prey. Thus, longer bills at the wet end of the gradient can alternatively be explained by a wider range of depths and substrates at which warblers need to search for insects.
In a phylogenetic study of the Yellow Warbler complex (Chaves et al., 2012), the authors report no morphological divergence among the populations on the Galapagos islands, even along environmental gradients. These populations are as old as Mangrove Warbler populations suggesting that selection driven by differences in precipitation, salinity, and forest structure is strong. Our results are comparable to other studies that have shown that bill size responds to intense shortterm selection and, therefore, can evolve rapidly (Boag & Grant, 1981;Eroukhmanoff et al., 2013;Smith & Dhondt, 1980). Alternatively, there is growing evidence that avian bill morphology plays an important role in heat exchange and thermoregulation, even when this trait is strongly associated with diet and foraging niche (Grant & Grant, 1993, 2002Symonds & Tattersall, 2010;Tattersall, Andrade, & Abe, 2009). Some studies have found that bill size may be important for heat dissipation in high humidity habitats (Friedman, Harmáčková, Economo, & Remeš, 2017;Gardner et al., 2016). Having a longer bill could also help Mangrove Warbler individuals with thermoregulation, especially at the wetter end of the gradient as humidity increases with precipitation, supporting our observations of larger bill sizes in wetter habitats.
Future studies should focus on the influence of bill divergence on song characteristics in Mangrove Warblers. Natural selection on phenotypic traits can pleiotropically cause divergence in sexually selected traits (e.g., song) which may trigger ecological speciation (Caro, Caycedo-Rosales, Bowie, Slabbekoorn, & Cadena, 2013;Laiolo & Tella, 2006;Schluter, 2001). It is well known that bill morphology limits the pace and timing at which bird songs are delivered.
Both song traits are important in interspecific recognition and sexual selection (Podos & Nowicki, 2004;Seddon, 2005;Slabbekoorn & Smith, 2000). Preliminary data on Mangrove Warblers show substantial vocal variability among individuals at the extremes of the rainfall/salinity gradient, which points to the hypothesis that these populations might be undergoing initial stages of a speciation process ( Figure S5).
Based on our P ST /F ST comparisons, we cannot argue that selection influences tarsus length, wing length, and body size. However, we did find a trend in which tarsus and wing length decreased along with precipitation. Both traits have been previously reported to respond significantly to environmental changes in other bird species (Milá et al., 2010;Pennycuick, 1968;Thomas, 1996). Shorter wings could benefit individuals in the wetter end of the gradient, as shorter wings are more efficient for maneuvering in areas with denser forest (Pennycuick, 1968;Thomas, 1996). Furthermore, Mangrove Warbler individuals in drier habitats could benefit from longer legs in order to expand the diversity of perches and foraging methods when the food is scarce (Janzen & Schoener, 1968), as longer legs allow the use of a greater variety of perches (Wright & Steadman, 2012).

| CON CLUS ION
Our findings highlight the importance of understanding both phenotypic and genetic variation along with their links, when examining population differentiation processes in populations which are not geographic isolated. The Setophaga petechia species complex has abundant geographic intraspecific variation based on plumage color and pattern (Browning, 1994). At least nine subspecies are recognized in the aestiva group, eighteen in the petechia group, and sixteen in the erithachorides group. Such high variability might indicate that there is incipient diversification and our study suggests that such differentiation might be caused by environmental variability. Further studies in other groups within S. petechia can lead to a better understanding the early stages of the formation of biological diversity in a group in which numerous populations could potentially constitute incipient or full biological species.