Drivers of phenotypic divergence in a Mesoamerican highland bird

Genomic data

We extracted genomic DNA from 18 frozen tissue samples collected between 1991 and 2016 (P. b. bidentata: 2♂, P. b. citrea: 2 unknown sex, P. b. flammea: 1♀6♂, P. b. sanguinolenta: 4♀3♂; Table 1, Shown in a map in Fig. S1) from the MZFC tissue collection, as well as samples donated by other institutions, using the Qiagen DNeasy (Valencia, USA) kit protocol, the EPICENTRE MasterPure kit protocol as well as standard phenol-chloroform extraction. Samples deposited at the MZFC collection were obtained under a field collection permit provided by Instituto Nacional de Ecología, SEMARNAT, Mexico (FAUT-0169). We included two Piranga roseogularis and one Cardinalis cardinalis sample to act as outgroups in phylogenetic analyses. nextRad genotyping-by-sequencing libraries were obtained from our genomic DNA samples following the procedure described in Russello et al. (2015) at the SNPsaurus laboratory in the University of Oregon. Genomic DNA was first fragmented with Nextera reagent (Illumina, Inc.), which also ligates short adapter sequences to the ends of the resulting fragments. Since some of the samples had a high amount of degraded DNA, the Nextera reaction was scaled for fragmenting 20 ng of genomic DNA but 40 ng of genomic DNA was used for input to compensate, as well as to increase fragment sizes. Fragmented DNA was selectively amplified for 27 cycles at 74 degrees, with one of the primers matching the adapter and extending 10 nucleotides into the genomic DNA with the selective sequence GTGTAGAGCC. The nextRAD libraries were sequenced on a HiSeq 4000 platform with one lane of single-end 150 bp reads at an average depth of 54× (University of Oregon). Raw sequence reads are available at GenBank SRA (BioProject accession PRJNA757783).

We assembled the resulting reads into a SNP library using ipyrad (Eaton & Overcast, 2020) using the reference approach and mapping our reads to the Northern Cardinal (Cardinalis cardinalis) reference genome (Sin, Lu & Edwards, 2020, GenBank SRA BioProject accession PRJNA642398). We filtered reads for quality by trimming reads with a PHRED score of less than 43 from the 3′ end. Since many biological processes can affect the level of genetic differentiation between two samples, the clustering threshold at which two sequences are called as orthologous cannot be assumed by default, and must be optimized empirically for each study system (McCartney-Melstad, Gidiş & Shaffer, 2019). We selected a clustering threshold of 88% based on four metrics: (1) Pearson correlation coefficient between pairwise genetic distance and percentage of missing data, (2) cumulative variance explained by the first three Principal Components, (3) total SNPs recovered, and (4) total loci recovered (McCartney-Melstad, Gidiş & Shaffer, 2019). For the final alignment, we retained loci that were present in 75% of the samples. For analyses that require unlinked loci, we further filtered our dataset using the R packages SNPfiltR (DeRaad, 2021) and vcfR (Knaus & Grünwald, 2017), retaining only loci that were more than 1000bp away from one another.

Color data

We measured the color of seven feather patches (Crown, nape, tail throat, breast, upper belly and lower belly; shown in Fig. 1C) of 109 study skins (Table 1) held at the AMNH Ornithology Collection and Museo de Zoología de la Facultad de Ciencias (MZFC), using an OceanOptics USB 2000+ spectrophotometer paired with a pulsed xenon light source. Six of these patches are carotenoid-colored, and one (tail) is melanin-colored. The measuring probe was directed perpendicularly to the feather surface and ambient light was blocked by a drilled rubber stopper surrounding the tip of the probe. We obtained spectrographic measurements of feather reflectance in the avian visual range of 300-700 nm. Feather patch reflectance was measured relative to the reflectance of a WS-2 white standard (OceanOptics) and then processed in the R package pavo (Maia et al., 2019) to smooth the spectrographic curves and correct for electric noise. From the smoothed and noise-corrected spectrographic data of each feather patch we obtained three standard variables, hue (H3), chroma (S8) and absolute brightness (B1). We used H3 instead of the more commonly used H1 because the spectral curves of carotenoid coloration show a distinct “shoulder” shape, always having its peak (R_max) at 700nm, showing no variation and rendering the variable uninformative. We found that calculating the peak intensity at R_mid instead of R_max (H3) better reflects the observed variation in carotenoid coloration. In addition, we also obtained the red chroma (SR1) of each patch, calculated as the proportion of the area under the curve in the red region (between 580 nm and 700 nm), as a complementary measure of both hue and red spectral purity. We followed the same procedure for the tail patch, which shows melanin-based coloration, to make results comparable between feather patches. We modeled the visual response to the measured color under the noise-limited visual receptor model described by Vorobyev & Osorio (1998) as implemented in the R package pavo (Maia et al., 2019) to test if any observed differences in color translate into perceptual differences in the avian visual system. The visual model we used to simulate the quantum catch was that of the Blue Tit (Cyanistes caeruleus), as it is the phylogenetically closest to our study species from models available in the package. The source of illumination we chose was sunlight filtered through forest canopy, as it most accurately reflects the light conditions our study species is likely to be seen under in the wild.

Morphometric data

We examined 141 study skins of Piranga bidentata, held at the American Museum of Natural History (AMNH) and the Museo de Zoología de la Facultad de Ciencias (MZFC) zoological collections, covering all four-described subspecies and the majority of the species’ geographical distribution (Table 1). Since juvenile males have plumage that strongly resembles that of females, we removed from analysis any specimen with yellow plumage that lacked age data that identified it as an adult (e.g., skull ossification, presence of bursa or notes on gonad development). This also allows us to mitigate the error caused by allometric changes during different stages of individual development. We used a Mitutoyo digital caliper to take external morphological measurements of standard body parts. To reduce measurement error, we used the average of triplicate measurements of each body part for analysis. We discarded the measurements of tail length because we found it difficult to replicate and highly affected by terminal wear on feathers.

Climate data

We used the R package raster (Hijmans, 2019) to extract the 19 climate variables from the CHELSA 2.1 climatology data set (Karger et al., 2017; Karger et al., 2021) in a circular 1km buffer around the geographic coordinates of each voucher specimen. For population-level analyses, we calculated a convex hull containing all specimens belonging to that population and obtained the climate variables around the hull centroid. These climatic variables correspond roughly to temperature, precipitation and seasonality variables, and many are highly correlated. To reduce model overfitting and mitigate the risk of multicollinearity, we avoided using variables with more than |0.7| correlation in the same model as suggested by Dormann et al. (2013). To test whether there was a correlation between feather color and vegetation cover, we obtained an estimate the percentage tree cover from the Moderate Resolution Imaging Spectroradiometer (ModIS) satellite though their MOD44B.006 Terra Vegetation Continuous Fields Yearly 250 m dataset (DiMiceli et al., 2015), accessed through the Google Earth Engine portal (Gorelick et al., 2017) and operated by NASA. We calculated the median vegetation cover from 19 yearly images spanning the period from 2000 to 2018. Finally, we calculated pairwise geographic distance between each collection point and between population centroids using the great circle distance algorithm implemented in the R package geodist (Padgham & Sunmer, 2019).

Analysis

To first assess if genetic and color variation are phylogenetically structured (color groups corresponding to described subspecies), we constructed a phylogenetic tree with our SNPs data. We used the SVDQuartets (Chifman & Kubatko, 2014) algorithm implemented in PAUP* (Swofford, 2002) to obtain a coalescent tree topology. We generated 100 bootstrap replicates to assess branch support in the resulting consensus phylogenetic reconstruction, which we present in Fig. 2. We then assessed visually if members of the described subspecies clustered together into monophyletic groups. We also assessed these groups using the Bayesian clustering algorithm implemented in the software ADMIXTURE (Alexander, Novembre & Lange, 2009), using a K determined by minimizing cross-validation error. This allowed us to estimate the ancestry of individuals within our sample.

As the study species is obviously dichromatic, we performed separate analysis on colorimetric data by sex. We explored variation in color by performing one-way analysis of variance (ANOVA) on the calculated color variables (B1, S8, H3 and S1R) with the subspecies as a factor. The results of this analysis are in the Tables S1 and S2. To explore if intraspecific variation is consistent with Gloger’s Rule we constructed linear models with plumage brightness (B1) as a response variable and mean annual precipitation (bio12) and mean annual temperature (bio1) as predictors. Ventral coloration is under stronger sexual selection compared to dorsal coloration in many bird taxa (Shultz & Burns, 2017; Marcondes & Brumfield, 2019). Since ventral color is relatively uniform patterned in this species, we took the average of the four ventral patches as representative of overall carotenoid color in subsequent analyses.

To further assess if genetic clusters match color groups, we performed linear discriminant analyses on ventral coloration variables on females (using plumage brightness, hue and chroma; B1, H3 and S8) and males (using plumage brightness, hue, overall chroma and red chroma; B1, H3, S8 and S1R) separately.

We followed a similar procedure for morphometric data, using one-way ANOVA to test for differences in each measured variable between subspecies (Table S3). Since one morphological variable was one order of magnitude larger than the rest, we rescaled all morphological data by log-converting it. To test if there were differences in body size among the described subspecies, we carried out an analysis of variance (ANOVA) with subspecies as a factor for each of the variables we measured. The results of these tests are shown in Table S3. To condense this morphological variation, we carried out a principal component analysis (PCA) on our morphological data matrix and obtained a set of principal components, the first of which explains 79.4% of the total variation, and is positively correlated almost entirely with wing chord. The second principal component explains 12.76% of variance and is positively correlated with the rest of the morphological variables. To test if this species follows the pattern known as Bergmann’s Rule, we constructed a linear regression model with the first PCA score value as a response variable and mean annual temperature (bio1) as the predictor.

To distinguish between the effect of genetic distance, geographic distance, and environmental variables, we used multiple regression on matrices (MRM) as implemented in the R package ecodist (Goslee & Urban, 2007), which tests for relationships between two or more distance matrices using matrix permutation. We also tested the relative effects of geographic distance and rainfall divergence on allele frequencies using a Bayesian MCMC framework implemented in the R package BEDASSLE (Bradburd, Ralph & Coop, 2013). We ran BEDASSLE for 1,000,000 generations using the beta binomial model, which allows populations to diverge from model expectations due to overdispersion. We used a matrix of Euclidean distances on rainfall, calculated from the CHELSA climatology data set (Karger et al., 2021). To test for the effects of selection, we used Bayescan 2.1 (Foll & Gaggiotti, 2008) to identify loci under selection in the color morphs. These tests allowed us to test the likelihood of the null model of isolation by distance (IBD), divergence due to drift, and local adaptation. Since we lack a one-to-one match between genomic sequences and color data, we carried out this analysis using the average spectrum per subspecies and an Fst distance matrix between the four subspecies.

Genomics

From our data set of 84,739 SNP loci we obtained a coalescent tree topology (shown in Fig. 2A) which grouped all of the Piranga bidentata samples into a single well-supported clade, supporting the monophyly of Piranga bidentata. The two western groups with orange male plumage (bidentata and flammea) cluster into one subclade with middling support. Samples from populations with red male plumage cluster together, although the topology doesn’t reflect the geographic structure of the samples, and we did not recover the subspecies P. b. sanguinolenta as monophyletic. Despite our limited sampling, we found some evidence for genetic structure in this species.

Our ADMIXTURE results (K = 2, Fig. 2B) show that most of our samples cluster unambiguously into two groups. These groups correspond to the orange color morph (subspecies bidentata and flammea) and the red color morph (subspecies citrea and sanguinolenta). Our two samples from the bidentata subspecies appear in the analysis as admixed between these two groups.

We performed two Bayescan runs of 100,000 generations each on a subset of 8,440 unlinked loci, one with prior neutral model odds of 100 and the other with odds of 1,000. In both runs we failed to detect any loci as statistically significant outliers in Fst, suggesting the markers in our data set are selectively neutral. However, this method risks low statistical power with small sample sizes, so these results should be interpreted carefully. A more extensive sampling is needed to fully assess selection in this taxon.

Our BEDASSLE run resulted in an aE/aD posterior distribution with a mean of 15,091 which suggests that the effect of rainfall on allele frequencies is relatively small compared to the effect of geographic distance (A difference of one mm of annual rainfall having a comparable effect as about 15km of lateral movement).

Color

The ANOVA of the color data (Table S1) showed a break in males’ carotenoid coloration for all four of the variables we computed (H3, S8, SR1 and B1), with males of the orange group (subspecies bidentata and flammea) having brighter, more orange and less spectrally saturated plumage color than males from the red group (subspecies sanguinolenta and citrea).

Our LDA result shows that for males (Fig. 3A), ventral color variables B1, H3, S8 and S1R distinguish two color groups congruent with the genetic clusters we detected. However, these predictors can’t distinguish between the described subspecies that form each of the groups. We found a similar pattern for females (Fig. 3B), although the separation is much less clear.

Under the noise-constrained visual receptor modeling (Vorobyev & Osorio, 1998), the average difference in color that we measured between the red and orange groups is higher than the perceptual threshold (i.e., larger than 1 JND, for just-noticeable difference Fig. 4A), which indicates that the colors are sufficiently different that birds are capable of discriminating between them. Notably, there was also a slight but noticeable and statistically significant difference between the West mainland and the Tres Marias orange birds, with the island birds having brighter orange plumage. The pairwise differences in ventral color between subspecies of the same color group (bidentata-flammea and sanguinolenta-citrea) are much smaller.

Females were generally consistent in their ventral carotenoid coloration, showing no significant variation in either hue or chroma among subspecies (Table S2). We found a small but statistically significant variation in plumage brightness among subspecies, which follows the same overall pattern as in the males. Although the pairwise differences in overall carotenoid color between subspecies are larger than one 1 JND, the threshold falls within the confidence interval (Fig. 4B).

For melanin coloration, we found a similar pattern and magnitude of divergence between subspecies in both males (Fig. 4C) and females (Fig. 4D).

For males, we found a significant negative relationship between plumage brightness (B1) and annual precipitation (bio12), for both carotenoid (ventral, F_1–73 = 26.6, adj. R² = 0.26, p > 0.001, Fig. 5) and melanin (tail, F_1–73 = 13.62, adj. R² = 0.15, p < 0.001) feather patches. We found a similar negative trend in female color, although the coefficients of determination were lower (ventral, F_1–48 = 12.62, adj.R² = 0.19, p < 0.001; tail, F_1–47 = 6.34, adj. R² = 0.10, p = 0.015). These results show that this species follows the ecogeographical pattern known as Gloger’s Rule, which predicts that among ectothermic species, individuals found in habitats with higher humidity will tend to have darker coloration.

We found a similar negative relationship between ventral plumage brightness and isothermality bio3 (males, F_1–73 = 27.42, adj. R² = 0.26, p > 0.001; females, F_1–48 = 8.76, adj. R² = 0.14, p = 0.005), as well as a statistically significant but weak negative relation between tail brightness and isothermality (F_1–73 = 7.22, adj. R² = 0.08, p = 0.009) in males.

To further test whether there was support for the hypothesis that Gloger’s Rule is at least partially due to selection of darker individuals in habitats with more dense vegetation, we tested for a relationship between plumage brightness and remotely sensed vegetation coverage calculated from ModIS measurements. We found a significant relationship between plumage brightness (F_1–60 = 16.55, adj. R² = 0.2, p < 0.001) and vegetation cover, but not for red chroma (F_1–60 = 1.37, adj. R² < 0.001, p = 0.246) and vegetation coverage, as calculated from ModIS measurements. The results of MRM analysis showed that for males there is statistically significant positive relationship between differences in breast feather color and geographic distance (F = 73.26, R² = 0.037, p = 0.001), and no significant relationship between breast color differences and genetic distance (F = 1.97, R² = 0.33, p = 0.38) or difference in vegetation cover (F = 9.26, R² = 0.004, p = 0.26). However, we found a significant positive relationship between differences in mean annual rainfall and differences in breast color (F = 439.229, R² = 0.19, p = 0.001). For females, we found no relationship between color distance and geographic distance (F = 16.064, R² = 0.015, p = 0.128) or differences in rainfall (F = 1.92, R² = 0.002, p = 0.544).

The MRM model with male breast color distance as response variable and Fst, geographic distance and rainfall distance matrices as joint predictors shows a not statistically significant but very large proportion of color variance is explained by these three factors (F = 26.39, R² = 0.975, p = 0.07). We obtained a similar result with female breast color distance, although the trend is much weaker (F = 0.50, R² = 0.43, p = 0.79)

Morphometrics

We found no significant sexual dimorphism in this group for five of the six morphological variables we measured. The only measurement for which sexual dimorphism was statistically significant was the wing chord. Of the six measurements, we found significant variation between subspecies in 5 of them (Table S3). The first two Principal Components resulting from our PCA explained 73.48% of the observed variation. Principal Component 1 was strongly positively correlated with wing chord, and Principal Component 2 was correlated with the other morphological variables we measured, and thus serves as a proxy for overall body size. We found that PC1 scores separate the subspecies bidentata from the rest, while PC2 separates the Mexican subspecies (bidentata and sanguinolenta) from the remaining two. However, all subspecies showed a high degree of overlap in morphometrics (See Supplemental Information). MRM results showed a small but significant relationship between our proxy for body size (PC2) and geographic distance (F = 202.57, R² = 0.02, p = 0.001) but not between PC2 and mean annual temperature (F = 0.3, R² = 0.0003, p = 0.84).

Plumage color

Plumage carotenoids are obtained exclusively from the diet, but the relationship between carotenoid ingestion and plumage color expression is not necessarily direct. The deposition of red and orange carotenoids (e.g., canthaxanthin and astaxanthin) requires metabolic conversion through ketolation of dietary yellow carotenoids (e.g., zeaxanthin). Although the metabolic pathway is not fully understood, evidence from several studies suggests that CYP2J19, a member of the P450 cytochrome gene family, is implicated in the ketolation reaction of dietary carotenoids (Mundy et al., 2016; Lopes et al., 2016; Emerling, 2018; Sin, Lu & Edwards, 2020). The enzymatic activity of CYP2J19 and other cytochrome proteins appears to be highly correlated with overall mitochondrial activity (Hill et al., 2019; Cantarero et al., 2020), which potentially maintains the honesty of modified carotenoid ornamentation as a signal of individual quality. While Koch et al. (2019) found that the evidence linking carotenoid deposition to several measures of individual performance is weak, other studies show evidence that carotenoid pigmentation is correlated with testosterone levels (Khalil et al., 2020).

While the exact chemical makeup of the feather carotenoids in this species has not yet been studied, closely related species of the genus Piranga derive their yellow feather coloration from dietary zeaxanthin and their red and orange coloration from a combination of several 4-keto-carotenoids such as canthaxanthin and astaxanthin (Hudon, 1991; Lopes et al., 2016). This seems to be largely conserved in the genus, with a few exceptions (Hudon, 1991). Since the difference in color that we observed appears to be continuously variable, it is likely that observed differences in color result from changes in the relative proportions of several carotenoids being deposited in feathers. This ratio is likely associated to some degree with the activity of CYP2J19-like gene products, which perform these ketolation reactions.

Another potential source of variation in coloration is the effect of differential sexual selection between populations. Our pairwise comparisons (Fig. 3) show similar patterns and magnitude of variation between subspecies for female carotenoid (3B), male melanin (3C) and female melanin (3D), but much larger differences in male ventral carotenoid colors (3A). As ventral coloration has been found to be under heightened sexual selection compared to other plumage patches (Shultz & Burns, 2017), this greater variation could be the result of stronger sexual selection for ketocarotenoid coloration in the red group (sanguinolenta and citrea). A series of studies on foraging behavior and characterization of dietary carotenoids available to this species would help illuminate other possible mechanisms modulating color expression.

Gloger’s rule

We found that in Piranga bidentata, darker and redder coloration in males is correlated with higher environmental humidity (measured as annual rainfall), compared to orange coloration. This is consistent with the ecogeographical pattern known as Gloger’s rule, which in its original formulation predicts that for endothermic organisms, darker coloration is found in hotter and more humid environments (Gloger, 1833). This pattern is widely observed in birds and mammals, and a study by Zink & Remsen (1986). found that a majority of North American bird species studied conform to it to some degree.

One of the most widely accepted hypotheses to explain Gloger’s rule is selection for darker and more cryptic plumage in environments with dense vegetation and thus lower environmental light levels (McNaught & Owens, 2002). Alternative explanations involving other, non-mutually exclusive processes have also been proposed. These include increased feather resistance to mechanical wear (Bonser, 1995), increased resistance to bacterial and fungal degradation of feathers (Burtt Jr & Ichida, 2004; Gunderson et al., 2008), and pleiotropic effects from selection on other traits which increase with humidity such as immune responsiveness and overall mitochondrial activity (Johnson & Hill 2013; Hill et al., 2019; Cantarero et al., 2020). In general, these proposed explanations suggest that Gloger’s rule results from an increase in pigment load in feathers, stemming from selection on other traits. We observed that plumage brightness of both melanin and carotenoid-colored patches is negatively correlated with rainfall (though not with a measure of vegetation cover) as predicted by Gloger’s Rule, similar to what has been observed in Australian birds (Delhey, 2018) and in other American birds (Ramírez-Barrera et al., 2019). However, this contrasts with the results of Ribot et al. (2019), who found a strong association between plumage brightness and remotely sensed vegetation cover. Taken together, these results fail to support the hypothesis that darker plumage in this species is linked to a darker environment caused by denser vegetation coverage. An alternative explanation that is compatible with our results would be that the mechanism in play is selection on stronger sexual signaling and not increased crypsis. Somewhat counterintuitively, additional pigment load increases spectral purity (chroma) but decreases overall brightness, since pigments function by selective absorption of incoming light. As ambient light decreases, more pigment is necessary to achieve the same level of background contrast, but it also decreases the overall amount of light being reflected. Our data does not support the hypothesis that differences in plumage color are caused primarily by fixed genetic factors. A mutation which results in reduced function of CYP2J19-like genes or other genes involved in the carotenoid metabolism pathway would be unlikely to become fixed outside of founder effect scenarios, as these mutations appear to occur in very low frequencies in the wild and would be rather strongly selected against. We found that the color of males from the Tres Marias islands is significantly different from mainland males, and that this difference is strong enough to be distinguishable by avian visual systems. The color of Tres Marias males is achromatically brighter and less saturated compared to their mainland counterparts. This is consistent with the overall pattern we found, as the climatic conditions on the islands are more xeric than on the mainland, and the dominant vegetation presents a more open canopy. However, many island species across many taxa present reduced sexual signaling when compared to their mainland sister groups, and it has been hypothesized that this is due to changes in the intensity of sexual selection compared to the mainland.

We found no significant relationship between plumage brightness or red chroma and vegetation cover, which also shows limited support for the hypothesis that increase of red coloration is due to an increase in dietary carotenoid availability linked to the environment’s primary productivity as was the case in the work of Jacot et al. (2010). The strong positive correlation between feather color and rainfall suggests that the mechanism modulating carotenoid deposition is related to humidity, such as selection acting on immune response or mitochondrial activity. A direct measurement of mitochondrial activity or quantitative transcription of genes associated with immune response would allow us to directly test this hypothesis.

The lack of relationship between color variation and geographical distance indicates that the mechanism underlying this variation is not isolation by distance. This is in contrast to the findings of Núñez Zapata et al. (2018) and by Ramírez-Barrera et al. (2019), who used a similar approach and found that for Turdus assimilis and Habia rubica respectively, Isolation by Distance is the model that best explains color variation. This, coupled with the lack of genetic structure that correlates to geography suggests that there may be more gene flow than previously thought and/or that carotenoid-based coloration in this species is fairly plastic.

Coloration is a complex and likely highly modular phenotypic trait. While carotenoid-based coloration is undoubtedly under some degree of fixed genetic regulation, both experimental (e.g., McGraw & Hill, 2000; Jacot et al., 2010) and comparative studies such as Ramírez-Barrera et al. (2019) and this one show that it is relatively plastic and subject to a large degree of environmental modulation. Our results show that a single gross-scale environmental measurement like mean annual rainfall can have a disproportionate effect on the coloration of certain taxa. However, the exact mechanisms that mediate the observed color differences have yet to be tested, and experimental studies that explicitly test these hypotheses are needed.

Body size

Our results on body size are not conclusive; we found a statistically significant relationship between body size and geographic distance, suggesting that the mechanism driving body size variation is isolation by distance. This is similar to what was found in Cranioleuca antisiensis, where the mechanism underlying a dramatic change in body size between populations was found to be isolation by distance (Seeholzer & Brumfield, 2018). Similar results were also found by Ramírez-Barrera et al. (2019) in Habia rubica. These studies make apparent that variation in body size in these neotropical birds is mostly driven by geographic isolation. However, although the relationship we found was statistically significant, the effect size was very small and we would be cautious about extracting strong conclusions from it. A genomic scan of alleles linked to variation in body size would further clarify the role of genetic variation on phenotypic variation this species.

Taken together, our results show that coloration and body size evolution are seemingly decoupled in this species. In conclusion, we found that climate is an important driver of phenotypic divergence in this species. An experimental study manipulating availability of dietary carotenoids would help illuminate the proximate mechanisms underlying this divergence. Although our phylogenetic analysis shows that the two samples from Central America (ssp. citrea) separate into a well-supported clade, a more extensive sampling is required to clarify their taxonomic status.

Finally, although we find evidence of genetic structure in this taxon, we suggest that a broader sampling and more in-depth genomic analysis is necessary to fully assess the relationship between genomes, environment and phenotypes.