S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e L I F E S C I E N C E S copyright © 2024 the Authors, some rights Sequential introgression of a carotenoid processing reserved; exclusive licensee American gene underlies sexual ornament diversity in a Association for the Advancement of genus of manakins Science. No claim to original U.S. Haw Chuan Lim1,2*†, Kevin F. P. Bennett3,4*†‡, Nicholas M. Justyn5§, Matthew J. Powers5¶, Government Works. Kira M. Long6#, Sarah E. Kingston7, Willow R. Lindsay8, James B. Pease9**, Matthew J. Fuxjager10 distributed under a , creative commons Peri E. Bolton11,4**, Christopher N. Balakrishnan11,12, Lainy B. Day8, Thomas J. Parsons4, Attribution Jeffrey D. Brawn13, Geoffrey E. Hill5, Michael J. Braun3,4* Noncommercial License 4.0 (cc BY- Nc). In a hybrid zone between two tropical lekking birds, yellow male plumage of one species has introgressed asym- metrically replacing white plumage of another via sexual selection. Here, we present a detailed analysis of the plumage trait to uncover its physical and genetic bases and trace its evolutionary history. We determine that the carotenoid lutein underlies the yellow phenotype and describe microstructural feather features likely to en- hance color appearance. These same features reduce predicted water shedding capacity of feathers, a potential liability in the tropics. Through genome- scale DNA sequencing of hybrids and each species in the genus, we identify BCO2 as the major gene responsible for the color polymorphism. The BCO2 gene tree and genome- wide allele frequency patterns suggest that carotenoid-p igmented collars initially arose in a third species and reached the hybrid zone through historical gene flow. Complex interplay between sexual selection and hybridization has thus shaped phenotypes of these species, where conspicuous sexual traits are key to male reproductive success. INTRODUCTION implications for the trait’s evolvability. While genes underlying most Sexual selection has produced many remarkable physical and be- display traits remain unknown (4), rapid advances in genomics are havioral traits in animals (1). More broadly, sexual selection shapes creating unique opportunities to study sexual selection and reveal biodiversity by promoting speciation, driving rapid diversification of how it can work in concert—or at times, in conflict—with natural extravagant phenotypes, and accelerating molecular evolution (2, 3). selection to produce the diversity we see in nature (5). Gaining insight into the evolution of sexually selected traits includ- Parallel advances in speciation research are revealing the preva- ing ornamentation and other reproductive signals will be facilitated lence of hybridization and introgression in shaping the evolution by understanding their genetic architecture. Whether a trait depends of organisms (6). Introgression can enlarge the pool of genetic on few or many genes and the degree to which the underlying genes variation available to a species and assemble existing variants in are involved in other traits or metabolic pathways have important unique combinations, which in turn become the raw materials that fuel further evolution and diversification (7). In addition, genetic 1 elements underlying particular advantageous traits can sometimes department of Biology, George Mason University, Fairfax, vA 22030, USA. 2National Zoo and conservation Biology institute, Smithsonian institution, Washington, dc introgress across species boundaries despite processes that purge 20013, USA. 3department of Biology and Biological Sciences Graduate Program, foreign alleles in the genome (8). Although these topics have re- University of Maryland, college Park, Md 20742, USA. 4department of vertebrate ceived increased attention in the genomics era, important work Zoology, National Museum of Natural history, Smithsonian institution, Washington, dc 20013, USA. 5department of Biological Sciences, Auburn University, Auburn, AL remains to establish the frequency of introgressive hybridization 36830, USA. 6Program in ecology evolution and conservation Biology, University of and its role in driving phenotypic changes, adaptive radiation, diver- illinois Urbana-c hampaign, Urbana, iL 61801, USA. 7Sea education Association, Woods sification, and even reticulate evolution (9). Some of the perti- hole, MA 02543, USA. 8department of Biology and interdisciplinary Neuroscience Minor, University of Mississippi, University, MS 38677, USA. 9department of Biology, nent challenges include identifying introgressed traits that are Wake Forest University, Winston-S alem, Nc 27109, USA. 10department of ecology under positive selection, measuring their fitness effects, determin- evolution and Organismal Biology, Brown University, Providence, Ri 02912, USA. ing their genetic architecture, and deciphering the interplay be- 11department of Biology, east carolina University, Greenville, Nc 27858, USA. 12divi- tween positive and negative selection, which can limit the spread sion of environmental Biology, National Science Foundation, Alexandria, vA 22314, USA. 13department of Natural Resources and environmental Sciences, University of or persistence of introgressed elements in new environments or illinois, Urbana, iL 61801, USA. genomic backgrounds. *corresponding author. email: hlim22@ gmu. edu (h.c.L.); bennettkf@ si. edu (K.F.P.B.); To date, most research on introgressive hybridization has investi- braunm@ si. edu (M.J.B.) †these authors contributed equally to this work. gated traits whose introgression was driven by natural selection [e.g., ‡Present address: department of Biology, Pennsylvania State University, State college, (10–12)], with fewer studies focusing on traits involved in mate prefer- PA 16801, USA. ence or sexual selection [but see (13–16)]. One possible explanation §Present address: department of Biosciences, Swansea University, Swansea, SA2 for the apparent paucity of systems demonstrating such introgression 8PP, UK. ¶Present address: department of integrative Biology, Oregon State University, corvallis, is that sexual selection tends to be associated with accelerated pheno- OR, 97333, USA. typic evolution and increased discrimination among diverging species #Present address: center for conservation Genomics, Smithsonian National Zoo and on the basis of reproductive traits (17, 18). Cases where introgression conservation Biology institute, Washington, dc, 20008, USA. **Present address: department of evolution, ecology, and Organismal Biology, the of sexually selected traits has been identified deserve careful study, Ohio State University, columbus, Oh 43210, USA. both because their prevalence may be underappreciated and because Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 1 of 20 S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e they provide special opportunities to explore mechanisms of sexual Previous studies demonstrated a notable asymmetric intro- selection, such as perceptual bias (19) or honest signaling (20). gression of male secondary sexual traits between two Manacus The manakins (family Pipridae) are a group of birds in which species—M. candei (white-c ollared manakin) and M. vitellinus males are under strong sexual selection. Manakins exhibit lek mating (golden-c ollared manakin)—in northwestern Panama (13, 23, 29). systems, including colorful male plumage and elaborate mating In this region, the yellow collar plumage of male golden-c ollared displays (21). They also display interesting hybridization dynamics: manakins has introgressed beyond the genomic center of a hybrid There are at least nine intergeneric crosses known in manakins, sup- zone between the species, producing hybrid populations that look porting early observations of increased rates of intergeneric hybrid- like golden- collared manakins, but have genomic backgrounds of ization in avian clades with high sexual dimorphism (22). Manakins white- collared manakins (Fig. 1C). The spread of yellow collars into also display interspecific hybrid zones (13, 23–25) and one instance M. candei has not extended across the entirety of its range but has of apparent hybrid speciation (24), which has rarely been document- stabilized at the Río Changuinola, a river roughly 50 km from the ed in birds. Among manakins, the species in the genus Manacus per- center of the hybrid zone (23, 29). Other secondary sexual traits of form the most acrobatic courtship display (see movie S1) (26). The male golden-c ollared manakins, such as their olive- green belly four species and several subspecies of Manacus replace each other plumage, narrow hindneck collar, aggressive behavior on leks, and geographically in humid lowland forests throughout the Neotropics, at least one dance maneuver, have also introgressed asymmetrically and hybrid zones occur in three areas where their ranges abut but not to the extent of golden collars (13, 23, 29–31). Introgression (Fig. 1A). Males differ conspicuously in the colors of their collar and of these traits appears to be driven by female mate choice (28) and belly plumage, which include bright whites, vibrant yellows, saturated male- male competition (30), but the failure of these traits to sweep oranges, olive-g reens, and grays (Fig. 1B), and plumage coloration to fixation in M. candei populations suggests that countervailing has been shown to correlate with male mating success (27, 28). forces are at play (32). A B M. manacus (cis-Andean) M. manacus (trans-Andean) 100 100 M. vitellinus 100 M. aurantiacus 100 M. candei C Plumagetransition 2 4 3 D manacus candei Hybrid vitellinus aurant. 5 6 0.6 8 Genomic 0.4 hybrid zone 9 12 ~280 km 0.2 0 25 50 0.0 10 oll ar elly lar elly lar elly yB l B l lla r ell llar ll y Kilometers C Co Co B Co B Co B e Fig. 1. Overview of Manacus system. (A) distribution of major Manacus lineages. cross-h atched regions are hybrid zones between M. vitellinus and M. manacus (119). (B) the species tree of Manacus with maximum likelihood bootstrap support values shown (120). (C) hybrid zone between M. vitellinus (orange region) and M. candei (blue region), with plumage introgression shown in orange and blue stripes. Sampling sites for populations used in whole-g enome resequencing (circles), RAdseq (triangles), or both methods (squares). Site 12 is located at Gamboa in central Panama. (D) Lutein concentrations of Manacus collar and belly feathers. circles indicate replicate indi- viduals, and diamonds indicate means. Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 2 of 20 inola Río Chan gu Lutein (μg mg-1) S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e Because of the unusual introgression dynamics present in this whether lutein- pigmented Manacus belly feathers appeared olive- system and its occurrence in a group marked by intense sexual green or yellow is a function of melanin deposition in the barbules. selection, we sought to understand the evolutionary processes that The variable degree of melanosome deposition in hybrids reflects the gave rise to the diversity of plumage color in Manacus manakins. In fact that the geographic extent of olive- green belly color introgres- this study, we first characterized the biochemical and microstruc- sion trailed that of yellow collar coloration at the time the specimens tural contributors to color variation in male collar and belly plum- were collected in 1991 (23); however, belly color has introgressed age and then used whole- genome resequencing data to identify a further in the past 30 years (29). major gene involved in the yellow/white collar polymorphism and trace its evolutionary history in the genus. Feather microstructure Aspects of feather microstructure differed by species, plumage patch, and color. Carotenoid-p igmented collar feathers of M. vitellinus and RESULTS aurantiacus had significantly wider barbs than white collar feathers Feather pigmentation of M. manacus and candei (Table 2, Fig. 2A, and fig. S3). Hybrids’ Using quantitative high- performance liquid chromatography (HPLC) collar feathers had barbs of intermediate width, but they were sig- with carotenoid standards, we found only a single carotenoid pig- nificantly different only from aurantiacus. Barbule morphology ment—lutein—in male collar and belly feathers of Manacus species or also varied among the taxa. M. manacus and candei collar feathers hybrids (fig. S1). The concentration of lutein largely determined the had barbules typical of avian contour feathers, while vitellinus and coloration of collar plumage (Figs. 1D and 2A and table S2A). White aurantiacus collar feathers lacked barbules on roughly the distal collar feathers of M. manacus and candei had no detectable lutein, third of the feather (Fig. 2A). Hybrids had barbules that appeared while yellow collar feathers of M. vitellinus and birds in the plumage smaller than those of M. manacus or candei but not as greatly re- introgression region (hereafter, “hybrids”) had intermediate lutein duced as those of vitellinus or aurantiacus (Fig. 2A). levels. The orange collar feathers of Manacus aurantiacus had lutein To explore the effect of barb width on collar color, we modeled concentrations that were significantly higher than each of the other correlations between barb width and colorimetric variables from taxa (Table 1). In belly feathers, M. manacus had no lutein. Among reflectance spectra. Because barb width was positively correlated the taxa with lutein-p igmented bellies, the only significant difference with lutein concentration (P = 0.0071), we included lutein concen- was between M. aurantiacus, which had the highest concentration, tration as a covariate in each model. Even controlling for lutein, barb and M. vitellinus, which had the lowest (Table 1 and table S2A). width was significantly correlated with feather brightness and hue in We modeled the effect of lutein concentrations on colorimetric collar feathers, with wider barbs showing decreased brightness and variables calculated from reflectance spectra using pavo (Fig. 2, B more red- shifted hues (Table 2). and D) (33). For collar feathers, these analyses showed that lutein To further explore the effect of feather microstructure on color, concentration was significantly and positively correlated with yel- we compared yellow belly feathers of M. candei, which have fully low saturation, carotenoid (i.e., all yellow- to- red wavelengths) chro- formed barbules (Fig. 2C), to yellow collar feathers of vitellinus and ma, and hue (Fig. 2, E and F, and Table 1). Lutein concentration hybrids, which have reduced or missing barbules (Fig. 2A). There did not significantly predict collar feather brightness and was nega- was no significant difference in lutein concentration among these tively correlated with ultraviolet (UV) wavelength saturation (fig. S2 three (Table 1 and table S3A), but both M. vitellinus and hybrid col- and Table 1). These relationships were unchanged when belly lar feathers had significantly higher yellow saturation, higher carot- feathers (which contained visible melanosomes) were included in enoid chroma, and lower UV saturation—essentially, purer and more the analysis along with a melanization term (table S2D). saturated yellow coloration—than candei belly feathers (table S3). Using light microscopy, we qualitatively characterized the visible Microstructural changes can affect the physical properties of melanization patterns of Manacus collar and belly feathers. Typical feathers in addition to their color. A property of avian contour feathers feathers are bipinnate structures composed of a central rachis ex- that is especially important in the humid tropical forests where these tending from the skin, barbs extending from the rachis, and barbules manakins live is their ability to shed water. To investigate potential extending from the barbs (Fig. 2A). We found no visible evidence fitness consequences of the observed differences in feather micro- of melanin in collar feathers of any Manacus species or hybrids structure, we estimated water repellency of Manacus collar feathers (Fig. 2A). Yellow belly feathers of M. candei also contained no visible using a model from textiles science that has frequently been used in melanin (Fig. 2C). Belly feathers of trans-A ndean M. manacus (ssp. ornithology (34, 35). When applied to feathers, this model estimates abditivus), which appear gray to the naked eye, had white barbs and the beading of water on feather vanes using a simple structural barbules, with black stripes running the full length of the barbules index dependent only on measurements of barb width and barb (Fig. 2C). We presume these stripes to be clusters of melanosomes. spacing. Our water repellency estimates indicated that M. vitellinus Belly feathers of M. vitellinus, aurantiacus, and candei x vitellinus and M. aurantiacus, the two species with carotenoid- pigmented hybrids contained varying amounts of lutein in barbs and melano- collar feathers and accompanying microstructural modifications, somes in barbules, producing dark olive-g reen shades in vitellinus, had the poorest ability to shed water when compared to the species dull yellow to olive in aurantiacus, and yellow to dark olive- green in with white collars (Fig. 2G and fig. S3C), while hybrids showed hybrids (Fig. 2C). Belly feathers of M. vitellinus generally had heavily intermediate water shedding ability. Differences in barb width ap- melanized barbules. M. aurantiacus had some melanized barbules in peared to drive these differences in estimated water repellency the center of belly feathers, but the barbules became small, lost mela- among Manacus collar feathers (fig. S3A) because barb spacing nization, and were usually orange toward the distal end. In hybrids, was similar among taxa, with a statistically significant difference the presence of melanosomes in barbules was sporadic, making belly only between white-c ollared M. candei and M. manacus (fig. S3B feathers olive-g reen in some individuals and yellow in others. Thus, and table S4). Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 3 of 20 S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e M. manacus M. candei M. candeix vitellinus M. vitellinus M. aurantiacus A Rachis Barb Barbules 10× 10× 10× 10× 10× B 60 60 60 60 60 40 40 40 40 40 20 20 20 20 20 0 0 0 0 0 350 450 550 650 350 450 550 650 350 450 550 650 350 450 550 650 350 450 550 650 Wavelength (nm) C 40× 40× 40× 40× 40× D 60 60 60 60 60 40 40 40 40 40 20 20 20 20 20 0 0 0 0 0 350 450 550 650 350 450 550 650 350 450 550 650 350 450 550 650 350 450 550 650 Wavelength (nm) E F G Hybrid 540 0.85 M. vitellinus M. aurantiacus 530 Carotenoids 0.80 520 Hybrids * 0.75 β = 0.14, P = 0.045 510 β = 42.27, P = 0.003 No carotenoids 0.70 500 0.0 0.2 0.4 0.6 0.0 0.2 0.4 0.6 0 50 100 Lutein (μg mg- 1) Lutein (μg mg- 1) Water repellency (R + D)/R Fig. 2. Feather color, microstructure, and water repellency differences among Manacus species. (A and C) Microscope images of collar (A) and belly (c) feathers of Manacus species showing variation in color and microstructural morphology of barbs and barbules. M. vitellinus, M. aurantiacus, and, to a lesser extent, hybrids display reduced number and size of barbules toward the distal ends of the collar feather barbs. (B and D) Reflectance spectra of collar (B) and belly (d) plumage of Manacus spe- cies. (E and F) the modeled effect of lutein concentration in collar feathers on yellow- red wavelength chroma (e) and hue (F). AU, arbitrary units. (G) Predicted water repel- lency of Manacus collar feathers. “carotenoids” refers to M. aurantiacus and M. vitellinus, hybrids refers to M. vitellinus x M. candei hybrids, and “no carotenoids” refers to M. candei and M. manacus. the asterisk indicates a statistically significant (α < 0.05) difference between groups. Genetic basis of plumage color genomic similarity (figs. S4 and S5) and ongoing gene flow between To investigate the genetic basis of traits that have introgressed asym- these cross-r iver focal populations (tables S7 to S9 and figs. S6 to metrically across the vitellinus-c andei hybrid zone, we resequenced S8), the color difference between them, and the relatively deep whole genomes of two parental candei and four parental vitellinus divergence of the parental species (~0.97 million years; table S9) (sites 2 and 12, Fig. 1C), as well as 10 males from each of two focal combine to create ideal conditions to isolate the genetic basis of populations separated by less than 8 km across the Río Changuinola collar coloration. Several other M. vitellinus traits (belly color, collar (sites 3 and 4, Fig. 1C). Genomes of these focal populations are pre- width, and aggression at leks) have introgressed to varying degrees dominantly derived from candei (29, 36); but while males from site (13, 23, 29, 30) but were not fixed at site 4 on the east bank of the Río 3 have white collars, males from site 4 have yellow collars (23). The Changuinola when the samples sequenced here were collected. As a Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 4 of 20 Carotenoid chroma (AU) Hue (AU) S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e Table 1. Feather lutein variation and effects on color in Manacus. Results from statistical models that contrast lutein concentration in collar feathers among species, contrast lutein concentration in belly feathers among species, and estimate the relationship between lutein concentration and reflectance measurements in pigmented collar feathers. Lutein concentration data for each species are in table S2A. Bolded values correspond to contrasts or variables where differences or effects are statistically significant at P = 0.05. Lutein concentration - collar C ontrast Estimate (β) SE df t.ratio P value Lower CI Upper CI M . vitellinus −0.0281 0.0491 25 −0.572 0.9780 −0.172 0.116 - hybrid M. vitellinus - −0.2773 0.0491 25 −5.651 0.0001 −0.421 −0.133 M. aurantiacus M. vitellinus - 0.3648 0.0491 25 7.434 <0.0001 0.221 0.509 M. manacus M. vitellinus - 0.3648 0.0491 25 7.434 <0.0001 0.221 0.509 M. candei h ybrid - −0.2493 0.0491 25 −5.080 0.0003 −0.393 −0.105 M. aurantiacus hybrid - 0.3928 0.0491 25 8.005 <0.0001 0.249 0.537 M. manacus hybrid - M. candei 0.3928 0.0491 25 8.005 <0.0001 0.249 0.537 M. aurantiacus - 0.6421 0.0491 25 13.085 <0.0001 0.498 0.786 M. manacus M. aurantiacus - 0.6421 0.0491 25 13.085 <0.0001 0498 0.786 M. candei M. manacus - 0.0000 0.0491 25 0.000 1.000 −0.144 0.144 M. candei Lutein concentration - belly Contrast Estimate (β) SE df t.ratio P value Lower CI Upper CI M . vitellinus −0.0755 0.0354 25 −2137 0.2365 −0.1794 0.0283 - hybrid M. vitellinus - −0.1149 0.0354 25 −3.251 0.0248 −0.2188 −0.0111 M. aurantiacus M. vitellinus - 0.2261 0.0354 25 6.394 <0.0001 0.1222 0.3299 M. manacus M . vitellinus - −0.0934 0.0354 25 −2.641 0.0929 −0.1972 0.0105 M. candei hybrid - −0.0394 0.0354 25 −1.114 0.7975 −0.1432 0.0644 M. aurantiacus h ybrid - 0.3016 0.0354 25 8.531 <0.0001 0.1978 0.4054 M. manacus hybrid - M. candei −0.0178 0.0354 25 −0.504 0.9862 −0.1217 0.0860 M. aurantiacus - 0.3410 0.0354 25 9.645 <0.0001 0.2372 0.4448 M. manacus M. aurantiacus - 0.0216 0.0354 25 0.610 0.9722 −0.0823 0.1254 M. candei M. manacus - −0.3194 0.0354 25 −9.035 <0.0001 −0.4232 −0.2156 M. candei Effect of lutein concentration on pigmented collar feather colorimetric variables D ependent Estimate (β) SE P value Lower CI Upper CI variable Y ellow saturation* 0.07579 0.02559 0.00919 0.02154 0.13 Uv saturation† −0.09862 0.03237 0.00769 −0.1672 −0.03 c arotenoid 0.14216 0.06521 0.0446 0.003906 0.2804 chroma‡ M ean brightness§ −2.794 6.877 0.69 −17.3735 11.785 Mid hue¶ 42.273 12.102 0.00301 16.618 67.9271 *Pavo variable S1Yellow. †Pavo variable S1Uv. ‡Pavo variable S9carotchroma. §Pavo variable b2meanbright. ¶Pavo variable h3huermid. Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 5 of 20 S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e result, we expected the genomic signal of collar coloration to be frequency at site 4, where all males have yellow collars. There were strongest when comparing our focal sites. 192 such dSNPs across the entire 1.1-G b genome, 186 of which were We mapped the resequencing data to a reference genome for clustered in a few regions on four putative chromosomes (2, 20, 24, M. vitellinus, aligned reference scaffolds to chromosomes from an- and Z; Fig. 3A). No other chromosome had more than one dSNP other manakin assembly (Chiroxiphia lanceolata), and identified (data S1). We further filtered dSNPs to consider only those fixed at single- nucleotide polymorphisms (SNPs) that fit the allele frequen- either site 3 or 4. Of 53 such “fixed” dSNPs, 46 of them were in a cy pattern expected for genomic regions responsible for introgress- haplotype block of 72 kb on chromosome 24 (Fig. 3, A and I). This ing vitellinus traits. Specifically, we filtered the data for strongly block contains three genes, including Beta-C arotene Oxygenase 2 differentiated SNPs (“dSNPs”) that were missing no more than 50% (BCO2), which encodes an enzyme that cleaves yellow carotenoids, of genotypes in parental candei and vitellinus populations (sites 2 including lutein, into colorless products (Fig. 3, I and J) (37). The and 12, Fig. 1C), were fixed for different alleles between parental two other genes in this block (Fig. 3J and fig. S28) were TEX12, populations, had ≥80% M. candei allele frequency at site 3, where which is involved in spermatogenesis (38), and PTS, which is in- all males have white collars, and had ≥80% M. vitellinus allele volved in synthesis of an amino acid hydroxylase cofactor (39). A 25 dSNPs River “fixed” dSNPs F 1.00 Chr. 24 Chr. 2 0 B 0.5 F 0.75ST 0.50 0 C 3e-3 Da 0.25 0 0 D 3 RND Site 2 3 4 5 6 8 9 10 G 1500 0 E 0.8 fd 1000 500 0 Chr. 1 1A 2 3 4 4A 5 6 7 8 10 12 14 19 24 Z 0 I Homozygous candei Homozygous vitellinus Parental Site 2 3 4 5 6 8 9 10 Heterozygous Missing data candei H Belly 20 60 Collar Site 3 40 15 Site 4 20 Beard 0 10 Parental 0 40 80 vitellinus 1.4 1.5 1.6 1.7 1.8 1.9 Mb km from site 2 (M. candei) J PTS BCO2 TEX12 Genes dSNPs “fixed” dSNPs 1.62 1.63 1.64 1.65 Mb Fig. 3. Genomic evidence implicating BCO2 in the collar color difference between M. vitellinus and M. candei. (A) density of SNPs showing strong association with collar color in sites 3 and 4 (dSNPs) in 10- kb windows. Fixed dSNPs, overlaid in front, are a more strongly differentiated subset of dSNPs. (B to E) differentiation and intro- gression metrics averaged in 10-k b windows across the genome. Autosomal and Z chromosome 99.9% outliers are indicated in red. Scans of 25- and 50-k b windows produced similar patterns (fig. S8). (B) and (c) differentiation metrics FSt and Da, respectively, for sites 3 and 4. (d) and (e) introgression metrics relative node depth (RNd) and fd, respectively. (F to H) Genetic and phenotypic clines across the hybrid zone from sites 2 to 10. (F) clines for 4750 SNPs from RAdseq data. clines with centers dis- placed by more than 50 km from the main cluster of clines are highlighted in pink (chromosome 24), teal (chromosome 2), or black (other chromosomes). (G) histogram of all SNP cline center locations. (h) Phenotypic clines in collar color (% reflectance at 490 nm), belly color (% reflectance at 665 nm), and beard length. (I) individual geno- types for all SNPs fixed for different alleles between parental M. vitellinus (site 12) and M. candei (site 2) in the region of chromosome 24 around carotenoid metabolism gene BCO2. the bold black lines and pink shading indicate the inferred boundaries of the core block of dNA that has introgressed from M. vitellinus as far as site 4. (J) Loca- tions of dSNPs (black tick marks) and fixed dSNPs (orange tick marks) around the BCO2 gene. exons are shown as blue boxes. Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 6 of 20 # of cline centers Allele frequency Reflectance (%) Beard length (mm) S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e Table 2. Collar feather barb width variation and effects on color in Manacus. Results from statistical models that contrast barb width among species’ collar feathers and estimate the relationship between barb width and reflectance measurements in pigmented collars, controlling for lutein concentration. Bolded values correspond to contrasts or variables where differences or effects are statistically significant at P = 0.05. B arb width C ontrast Estimate (β) SE df t.ratio P value Lower CI Upper CI M. manacus - −1.97085 2.974 25 −0.6627 0.962556 −10.7045 6.762774 M. candei M. manacus −6.1235 2.974 25 −2.0590 0.268657 −14.8571 2.610125 - hybrid M. manacus - −14.6635 2.971 25 −4.9347 <0.0001 −23.3918 −5.93523 M. vitellinus M. manacus - −25.8807 2.974 25 −8.7021 <0.0001 −34.6143 −17.1471 M. aurantiacus M. candei - hybrid −4.15265 2.974 25 −1.3963 0.635563 −12.8863 4.580973 M. candei - −12.6926 2.971 25 −4.2715 0.002105 −21.4209 −3.96438 M. vitellinus M. candei - −23.9098 2.974 25 −8.0394 <0.0001 −32.6434 −15.1762 M. aurantiacus hybrid - −8.54 2.971 25 −2.8740 0.0573 −17.2683 0.188264 M. vitellinus h ybrid - −19.7572 2.974 25 −6.6431 <0.0001 −28.4908 −11.0236 M. aurantiacus M. vitellinus - −11.2172 2.971 25 −3.7749 0.007185 −19.9454 −2.48892 M. aurantiacus Effect of barb width on pigmented collar feather colorimetric variables Dependent Estimate (β) SE P value Lower CI Upper CI variable Yellow saturation* 1.91 × 10−4 5.41 × 10−4 0.7296 −9.63 × 10−4 1.34 × 10−3 Uv saturation† −6.39 × 10−4 6.59 × 10−4 0.3479 −2.04 × 10−3 7.66 × 10−4 carotenoid −1.18 × 10−3 1.38 × 10−3 0.4084 −4.12 × 10−3 1.77 × 10−2 chroma‡ Mean brightness§ −0.2986 0.1272 0.0331 −0.5698 −0.0274 M id hue¶ 0.5448 0.2135 0.0221 0.0898 0.9998 *Pavo variable S1Yellow. †Pavo variable S1Uv. ‡Pavo variable S9carotchroma. §Pavo variable b2meanbright. ¶Pavo variable h3huermid. To test the geographic concordance of SNP markers with intro- concordant signals among them increase the likelihood of detecting gressed plumage traits, we examined a separate restriction site– true positives (40, 41). We therefore used FST and Da (42) to gauge associated DNA sequencing (RADseq) dataset from 152 individuals levels of genetic differentiation between focal sites 3 and 4, and across the entire hybrid zone (Fig. 1C and table S14). From this we used fd (43) and relative node depth (RND) (44) to infer intro- genome-w ide dataset, we retained 4750 high-q uality SNPs show- gression from parental M. vitellinus (site 12) into the yellow- collared ing clear clinal variation in allele frequency within the hybrid zone hybrid population (site 4). We treated as outliers all 10-k b windows transect. Cline centers for the vast majority of these SNPs fell near with extreme values for at least one statistic exceeding the 99.9th site 9 (Fig. 3F), marking the major genomic transition from M. vitellinus percentile (Fig. 3, B to E). These analyses revealed several genomic to M. candei, in agreement with previous analyses (23, 36). How- neighborhoods with high concentrations of strongly supported win- ever, cline centers for 34 SNPs were located more than 50 km to dows, each containing outliers for multiple statistics (fig. S9). These the northwest (Fig. 3G), falling near the Río Changuinola and neighborhoods included regions in chromosomes 2 (690 kb) and 24 the limit of color introgression (Fig. 3H). Of these 34 SNPs with (140 kb around BCO2) and the Z sex chromosome (590 kb) (Fig. 3, displaced cline centers, 10 are clustered in a genomic span of B to E; figs. S10 to S33; and table S10; see Supplementary Text for 768 kb on chromosome 24 surrounding the BCO2 gene (Fig. 3F details). In addition to carotenoid metabolism, genes adjacent to or and data S2). within strongly supported genomic windows span a wide range of To further localize DNA blocks that have introgressed, we con- functions related to angiogenesis, feather/follicle development, ducted sliding window scans of the resequencing data, calculating solute transport, and spermatogenesis (table S11). Gene ontology– metrics of differentiation and introgression in non-o verlapping based analyses revealed pathways or molecular functions that were 10- kb genomic windows (Fig. 3, B to E). The use of multiple enriched with genes found within or around outlier genomic win- summary statistics leverages different aspects of the data, and dows (tables S12 and S13). Among these pathways were several with Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 7 of 20 S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e plausible connections to traits of interest, but none survived multi- and 384 kb around two other genes also on chromosome 24 but ple testing correction. See Supplementary Text for further discussion. with no known relationship to coloration (Fig. 4B). Topologies of BCO2 is the strongest candidate for a gene of major effect driving the control gene trees (Fig. 4, D and E) matched the strongly sup- introgression of collar color in this system. The genomic region ported species tree (Fig. 1B), with M. aurantiacus sister to candei around BCO2 contained multiple windowed outliers, displaced SNP and vitellinus sister to trans- Andean manacus. In the BCO2 gene clines, and numerous SNPs fixed for the vitellinus allele in site 4 tree, however, vitellinus was sister to aurantiacus, and the vitellinus- and almost fixed (except for two heterozygous individuals) for the aurantiacus clade was sister to candei, both with strong branch sup- candei allele in site 3 (Fig. 3I). BCO2 has been implicated in carotenoid port (Fig. 4C). The fact that vitellinus moves from its position on color polymorphisms in a variety of organisms [e.g., (45–47)], in- the species tree to become sister to aurantiacus in the gene tree cluding cases of interspecific introgression (48, 49). Because BCO2 rather than the reverse indicates that the BCO2 sequence now cleaves yellow carotenoids into colorless products when active (37), present in vitellinus originated in the aurantiacus genomic back- we expect the pigmented collar phenotypes of Manacus to be associ- ground (Fig. 4A). ated with reduced BCO2 activity in relevant tissues during feather Because discordance between gene trees and species trees can growth, either due to reduced gene expression or inactivity of the arise from several processes, we also looked for genome- wide evi- enzyme (48,  50–52). The amino acid sequences predicted by the dence of introgression between vitellinus and aurantiacus using al- white-c ollared and yellow-c ollared haplotypes are identical, and lele frequency patterns. We calculated fd in 10- kb windows using nucleotide substitutions were concentrated in introns and noncod- trans- Andean M. manacus, vitellinus, aurantiacus, and cis-A ndean ing regions 5′ or 3′ of the gene (Fig. 3J). Thus, mutations that are manacus as test taxa and identified outliers exceeding the 99.9th functionally relevant to the yellow phenotype likely cause reduced percentile. A window 10 kb upstream of BCO2 had the fifth highest gene expression in M. vitellinus either constitutively or in key tissues fd value of ~78,000 genome- wide windows (Fig. 4, F and G). at crucial times for feather growth. Available RNA sequencing data To localize the signal of introgression from M. aurantiacus to showed that BCO2 expression levels were higher in brain and lower vitellinus at a finer scale, we repeated the fd calculation in 2- kb win- in muscle tissues of M. vitellinus when compared to those of related dows across chromosome 24 (Fig. 4I). We also estimated local gene species, suggesting that the gene is under tissue-s pecific regulation trees based on 100- bp non- overlapping windows and calculated the (fig. S34). frequency of trees in 10- and 2-k b windows showing the topology Our suite of analyses identified several other genes with plausi- expected under gene flow, with M. aurantiacus sister to vitellinus ble connections to introgressing traits across the hybrid zone. Of (Fig. 4, H and J). These two metrics—fd and frequency of “gene flow” 192 dSNPs in the genome (data S1), 4 fell in a region on chromo- topologies—identified different stretches of DNA in the vicinity of some 20 between 25 and 114 kb upstream of Agouti Signaling Pro- BCO2 as containing the strongest signals, with the highest fd peak tein (ASIP), which encodes a polypeptide involved in regulating occurring just upstream of the gene and the highest proportion of melanogenesis and commonly underlying melanin- based color gene flow trees occurring just downstream of the gene (Fig. 4, J and variation [e.g., (53,  54)]. Sliding windows (table S10) and dSNPs I). Together, these results all support directional gene flow from (data S1) both implicated Retinol Dehydrogenase 10 (RDH10) on M. aurantiacus into vitellinus at BCO2, including a considerable in- chromosome 2. This gene has not previously been associated with trogression signal within the same block that has introgressed across color variation, but it encodes an enzyme involved in the linked the M. vitellinus- candei hybrid zone. metabolism of carotenoids and retinoids (37). A third region of in- terest, identified by dSNPs, sliding windows, and shifted RADseq clines, fell on chromosome 2 between the genes R-s pondin 2 DISCUSSION (RSPO2) and Eukaryotic Translation Initiation Factor 3 subunit E Feather color and microstructure (EIF3E) (table S10 and data S1 and S2). RSPO2 has the most obvi- Through HPLC, spectrophotometry, and light microscopy, we ous connection to plumage: It is a regulator of Wnt (55), which has found that color production in belly and collar feather patches been implicated in other instances of plumage patterning (56), and among the major lineages of Manacus can be explained largely by RSPO2 controls major variation in domestic dog hair growth (57), the concentration of the carotenoid lutein, the presence or absence which could indicate a connection to follicle development. For of melanin, and feather microstructure, including barb morpholo- more discussion of introgressing genes, see Supplementary Text. gy and the presence or absence of barbules. Differences in lutein concentrations of pigmented collar feathers History of BCO2 gene flow in Manacus affected saturation, chroma, and hue, as predicted by simulations of The two Manacus species with carotenoid- pigmented collars carotenoid absorbance (60). Orange feathers of M. aurantiacus con- (M. vitellinus and aurantiacus) are not sister taxa (Fig. 1B), but tained more lutein than yellow feathers of M. vitellinus and vitellinus both use lutein to produce their orange/yellow coloration (Fig. 2 x candei hybrids. Yellow bellies of M. candei and some hybrids con- and fig. S1). We wondered whether pigmented collars had a single tained only lutein, while olive-g reen bellies of hybrids, vitellinus, origin in this genus or might have arisen independently. The two and aurantiacus contained both lutein and melanized feather bar- species are not now in contact, but their ranges come within ~25 km bules. Gray bellies of M. manacus contained only melanin. Adding in central Panama (58,  59), and hybridization could have oc- or reducing carotenoids in a plumage patch may be one of the sim- curred through long- distance dispersal or during past periods of plest ways to create differences in color (60). On the basis of our contact. To test this question, we sequenced whole genomes of ad- light microscopy observations, we believe that adding melanized ditional manakins representing the remaining major lineages of the barbules could be another “simple” mechanism by which plumage genus (table S15). We constructed phylogenetic trees from a 48-k b color can be changed, although melanism is not associated with di- region centered around BCO2 and from “control” regions of 306 chromatism as often as deposition of carotenoids (61). Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 8 of 20 S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e A M. manacus M. manacus M. manacus (cis-Andean) (cis-Andean) (cis-Andean) M. manacus M. manacus M. manacus (trans-Andean) (trans-Andean) (trans-Andean) M. vitellinus M. vitellinus M. vitellinus M. aurantiacus M. aurantiacus M. aurantiacus M. candei M. candei M. candei No gene flow aurantiacus vitellinus vitellinus aurantiacus C 79 M. manacus F 0.8 87 (cis-Andean) 100 M. manacus (trans-Andean) 64 fd 0.4 60 32 BCO2 78 M. vitellinus63 0.0 99 86 Chr. 1 1A 2 3 4 4A 5 6 7 8 10 12 15 20 28 Z M. aurantiacus 100 Position in genomeB 100 G 0.8 100 M. candei 10-kb windows Substitutions per site: 0.001 D fd 0.4 M. manacus 1.6 Mb 95 (cis-Andean) 0.0 60 H 100 M. manacus2.8 Mb 0.06(trans-Andean) 80 100 91 0.04TBCEL 100 M. vitellinus 0.02 100 100 90 0.0 M. aurantiacus 0 Mb 2 Mb 4 Mb 6 Mb 5.4 Mb 100 Position on chr. 24 100 I 100 M. candei E 0.8 BCO2 2-kb windows 88 M. manacus 98 (cis-Andean) fd 0.4 47 M. manacus 62 (trans-Andean) J 0 DCPS 68 7886 0.2 M. vitellinus candeiM. vitellinus introgression block 98 0.1 69 100 M. aurantiacus 100 0 100 1.2 Mb 1.4 Mb 1.6 Mb 1.8 Mb 2.0 Mb 100 M. candei Position on chr. 24 Fig. 4. Evidence for historical BCO2 introgression between M. aurantiacus and M. vitellinus. (A) Gene trees expected under three different potential introgression scenarios. Lack of introgression results in a topology similar to the species tree. directional introgression results in a discordant topology where the gene flow recipient shifts to the branch of the gene flow donor. (B) Locations of BCO2 and two control genes (TBCEL and DCPS) on chromosome 24. (C to E) Gene trees generated using three regions on chromosome 24 for all major lineages in Manacus. Branch labels are maximum likelihood bootstrap support. (c) A 48- kb region including the gene BCO2, (d) a 384-k b region including the gene TBCEL, and (e) a 306- kb region including the gene DCPS. (F and G) introgression metric fd calculated between M. aurantiacus and M. vitellinus in 10-k b windows with 99.9% outliers in red. (F) the entire genome and (G) just chromosome 24. (H) Proportions of 100-b p trees in 10- kb windows in which M. vitellinus is sister to aurantiacus, which is the expected topology in a gene flow scenario. trees included M. vitellinus, M. aurantiacus, M. candei, and trans- Andean M. manacus. (I) fd and (J) gene flow topologies from (h) calculated in 2- kb windows in the region surrounding BCO2. the pink shaded regions correspond to the vitellinus- to- candei introgression block denoted by bold black lines in Fig. 3i. The impact of feather microstructure on carotenoid coloration is feathers develop without barbules, or if the barbules are gradually understudied, but the common co- occurrence of feather carot- lost throughout the season, as both processes have been observed enoids with wide barbs lacking barbules has been known for de- (63, 64). Lutein deposition in yellow feathers may inhibit full forma- cades (62). In accordance with these previous observations, we tion of hierarchical feather structure (63), but our results suggest found wider barbs and reduced barbules in Manacus collar feathers that it is not the sole factor causing structural differences in Manacus containing lutein and a significant positive relationship between feathers. Collar feathers of vitellinus had lutein concentrations simi- barb width and lutein concentration. It is unclear whether these lar to hybrid collars and candei bellies, but both of the latter retained Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 9 of 20 Chromosome 24 Freq(M.aur + M.vit) Freq(M.aur + M.vit) S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e more barbules than vitellinus collars. We also found that feather mi- well- supported BCO2 gene tree suggests that the allele for pig- crostructure, including barb width and, potentially, barbule mor- mented collars arose in aurantiacus and spread to vitellinus. While phology, affect color by shifting peak reflectance toward redder hues incomplete lineage sorting can also produce gene trees that do not and increasing color saturation and purity. These results suggest that match the species tree, we consider it unlikely to be the cause of the altering feather microstructure may be an important component of BCO2 gene tree topology. Incomplete lineage sorting is most likely birds’ color signaling and that this process is likely refined by selec- when internodes are short (68), but this is not the case for Manacus tion. Barb shape, including greater width and flatness, has previ- (69). In the study of Leite et  al. (69), all marker types (ultracon- ously been shown to increase color saturation in other birds (65), served elements and exons) and analytical methods (concatenation and a lack of barbules has been hypothesized to reduce diffuse and gene tree reconciliation) strongly supported the sister relation- scattering of light (63). Further analysis of morphological data in ship between candei and aurantiacus. M. vitellinus and aurantiacus manakins and other species, including more quantitative barbule are not now in contact, but they are separated by ~25 km in central measurement, will improve tests of this hypothesis. Populations that Panama (58, 59) and could have hybridized in the past. While we have received introgression resulting in de novo deposition of carot- cannot identify the phenotypic outcome of this introgression with enoids, as in the manakin hybrid zone, are particularly fruitful the same confidence as the vitellinus- candei introgression, a com- systems in which to conduct such tests. bination of evidence strongly supports the notion that BCO2 in- trogression from aurantiacus to vitellinus had an important impact Genetic basis of plumage color on vitellinus coloration. This evidence includes the connection we Three complementary genomic approaches, dSNPs, RADseq clines, discovered between BCO2 and vitellinus coloration, the known im- and sliding window differentiation and introgression metrics, all iden- portance of BCO2 to plumage color in other species (46,  47), the tified alleles around and within the gene BCO2 as the major genetic known importance of collar color for sexual selection in Manacus difference between cross- river manakin populations that differ in collar (27, 28), and the fact that the signal of introgression from aurantiacus color (Fig. 3). These analyses also confirmed that, in yellow- collared to vitellinus largely overlaps with the vitellinus-t o-c andei introgres- hybrid populations, the block of DNA containing BCO2 is derived sion block. The most likely scenario, therefore, is one in which the from M. vitellinus and reached admixed populations via introgres- BCO2 allele for pigmented collars arose once in M. aurantiacus, was sion. A growing body of evidence has pointed to BCO2 as a common subsequently transferred through hybridization to vitellinus, spread target of selection controlling simple carotenoid- based color differ- throughout its range, and then partially spread into a third species, ences in the integument of many vertebrates, including salmon flesh M. candei, in western Panama. (51), mammal fat (45), lizard skin (50), bird bare parts (48, 52, 66), Three- species introgression has rarely been reported in nature and bird feathers (46, 47). BCO2 alleles also appear to have spread (70–77). Even rarer are instances in which such introgression is between colorful species several times in parulid warblers (49). clearly connected to positive selection for the movement of traits When BCO2 is active, it cleaves carotenoids into colorless prod- between species, as occurred with the spread of seasonal camouflage ucts (37). As a result, changes to BCO2 that result in greater carot- pelage among rabbits (76) and, most likely, colorful plumage in enoid content typically involve loss- of-f unction mutations (45, 66) or Setophaga warblers (49). Multispecies introgression may be com- reduced expression (47, 48, 50). We found no differences between the mon early in rapid radiations like those of African cichlids (77), BCO2 coding sequences of yellow- and white- collared males, so loss Heliconius butterflies (78), and Andean siskins (79). Whether situa- of function at the protein level is unlikely. Nevertheless, the signal of tions like the present one involving sexually selected traits are com- introgression was localized in and around the gene, thus pointing to mon awaits further study, and the continued advance of genomic a mutation or mutations in cis-r egulatory elements in an intron or technologies and methods will facilitate that process. flanking sequence as the most likely genetic cause of the yellow phe- notype. Although BCO2 expression was lower in M. vitellinus muscle Sexual selection: What traits and which mechanisms? than several relatives, its expression in the cerebellum appeared to be The asymmetric introgression of adult male plumage traits in Manacus higher (fig. S34). Therefore, we favor a hypothesis involving precisely has been attributed to female choice and/or male- male competition localized and/or timed reduction in gene expression or enzyme ac- [reviewed in (32)]. In a study on mixed leks with both yellow- and tivity during feather development. Such tissue-s pecific gene expres- white- collared males, Stein and Uy (28) found that yellow-c ollared sion regulates other sexual traits in manakins (67), and skin- specific males received more visits from females and copulated more often control of BCO2 expression underlies the yellow skin of domestic than white- collared males, implicating female choice in the plumage chickens (48). Gene expression analysis in growing feathers, liver, introgression. Supporting the male- male competition hypothesis, kidney, gut, and other tissues involved in carotenoid processing and McDonald et al. (30) showed that yellow-c ollared manakin males, deposition is needed to address this hypothesis. whether pure M. vitellinus or hybrids, were more physically aggres- Several other genes were identified by one or a combination of sive than white- collared birds, which likely confers an advantage analyses that point to other instances of introgression across the during dominance interactions. However, since yellow collars and Manacus hybrid zone, including ASIP, RDH10, and RSPO2. Wheth- olive- green bellies were correlated in their study populations, it is er any of these genes are causally related to introgressing plumage or plausible that belly color was the plumage trait directly associated behavioral traits deserves further study. with increased aggression in that study. Stein and Uy (28), in a simi- lar study on mixed- color leks, did not detect a difference in aggres- Historical gene flow at BCO2 sion between yellow- and white- collared males, but this study also We also identified historical introgression at BCO2 between M. aurantiacus did not control for belly color. We have shown here that the darker and M. vitellinus. Allele frequency patterns showed a strong introgression olive- green belly colors of vitellinus and some hybrids result from signal at BCO2 between the two species, and the topology of the the deposition of melanosomes in feather barbules. In a number of Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 10 of 20 S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e vertebrate systems, melanocortins have been known to produce (34). Water repellency of hybrid collars spanned nearly the entire pleiotropic effects on both coloration and behavior, with darker col- range of values for white and pigmented collars (Fig. 2G), suggesting oration often associated with greater aggression, especially in birds that individuals at the lower end of this range may incur a fitness cost (80). Of note, Long et al. (29) recently showed that olive green belly (86) in the hybrid zone in western Panama, an area where rainfall is color, while initially trailing yellow collars in the Manacus system, extremely high (87). has continued to spread into admixed populations and has reached Collar width is another trait that has introgressed asymmetri- the Río Changuinola within the past 30 years. Thus, yellow collars, cally along with carotenoid pigmentation in this hybrid zone (13). olive bellies, and enhanced aggressiveness may all be targets of The pigmented collar on the dorsal surface (hindneck) of M. vitellinus sexual selection in this system, and their relative importance re- is much narrower than the white dorsal collar of M. candei (Fig. 5). mains to be determined. We speculate that this narrow hindneck collar trait may compensate In combination with previous results, evidence presented here for the reduced water repellency of lutein- pigmented collar feathers, that colorful collars likely originated in M. aurantiacus and spread limiting potential negative impacts of their reduced water shedding to vitellinus before reaching the vitellinus x candei hybrid zone and capacity. Previous work demonstrated that the narrow dorsal collar spreading into candei populations supports the hypothesis that col- of vitellinus has partially introgressed into hybrid populations but orful collars are generally attractive within Manacus. Positive selec- trails behind the limit of yellow collar introgression (13). Hybrids at tion on carotenoid- based plumage colors is implicated in other bird site 4 on the east bank of the Río Changuinola have bright yellow systems, such as the multiple cross- species introgressions of BCO2 collars like M. vitellinus, but their dorsal collar width is very broad among Setophaga warblers (49), the asymmetrical introgression of like M. candei (Fig. 5). It is possible that hybrids with this combina- red plumage in fairy-w rens (14), and the attractiveness of brighter tion of traits may experience negative fitness effects due to reduced red plumage in house finches (20). Selection for more intense colors water shedding capacity of their broad yellow collars. Natural selec- has been interpreted as potentially resulting from honest signaling tion may act to ameliorate this cost by reducing the width of the for higher male quality (81) or perceptual bias for particular signals dorsal collar, driving the observed introgression of narrow collars (19). These explanations for yellow plumage introgression would from vitellinus into candei. Future studies could test this idea by ex- predict that it should continue to spread whenever possible. amining the physiological effects of intense rainfall on hybrid mana- In Manacus, therefore, the fact that plumage introgression from kins with different combinations of these traits. M. vitellinus into candei appears to have been stalled for more than We note that our water repellency estimates are based on a sim- 30 years at the Río Changuinola (29) and has not swept through ple model originally designed to estimate the water shedding capac- trans- Andean M. manacus populations despite two apparent ity of porous textiles (34, 35). A recent review details complexities of contact zones in Colombia (Fig. 1A) (82) presents a conundrum. feather structure that are not captured by this model and empha- Several explanations for stalled introgression have been proposed, sizes the likely importance of barbule microstructure in water repel- including the river acting as a barrier to gene flow, positive frequency- lency (88). While barbule microstructure is not considered in our dependent selection for yellow, variable light environments, or a model estimates, it seems likely that the absence of barbules on the combination of these factors (23,  28,  83). To date, none of these distal third of pigmented Manacus collar feathers would only in- proposals is supported by convincing evidence (32). A recent mod- crease their wettability, reinforcing the conclusions reached here. el of mate choice focused on how females acquire mating prefer- In sum, our results show that a case of asymmetric sexual trait ences predicts that, under certain conditions, female preference introgression in the manakin genus Manacus likely results from a will vary over time and space (84). This possibility also deserves regulatory change in the carotenoid processing gene BCO2 and that consideration in the Manacus system. Evidence from the present this introgression was preceded by historical gene flow at the study linking BCO2 to collar color variation in the hybrid zone same gene from a third species in the genus. Such three-s pecies should facilitate further detailed investigation of the evolutionary introgression events may be more common than now appreciated mechanisms acting on color in Manacus. early in radiations, but reported cases involving key traits, like coloration in a lekking bird, have so far been rare. We show that the Do colorful collars come at a cost? We hypothesize that the spread of carotenoid- pigmented collars, while clearly advantageous for sexual selection in some cases, may M. vitellinus Hybrid M. candei come at a fitness cost for natural selection. We have shown here that microstructural modifications accompany lutein pigmentation in Manacus collar feathers and are likely to affect their physical perfor- mance in shedding water (see Materials and Methods for further discussion of the water repellency model). Water shedding is a criti- cal function of feathers in the tropical rainforests where manakins live. The potential negative fitness effects of reduced water shed- ding capacity are numerous and varied (85). Our estimates of water repellency predict that Manacus species with pigmented collar feathers have significantly reduced water shedding performance compared to Fig. 5. Collar width differences between M. vitellinus, M. candei, and their hy-brids. Specimens with representative dorsal collars from the hybrid zone region species with white collars. Lower predicted water repellency is due to were photographed: M. vitellinus are from site 10, hybrids are from site 4, and the wider barbs of pigmented feathers, which reduce pore size and M. candei are from site 2. Yellow feathers have reduced water shedding capacity increase the ratio of solid surface to air in contact with water drop- relative to white feathers, so we propose that the narrow collar of M. vitellinus may lets, reducing net surface tension and enhancing wettability of feathers compensate by reducing the available wettable surface. Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 11 of 20 S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e microstructure of carotenoid-p igmented feathers has likely been re- (hybrids, M. vitellinus, and M. aurantiacus). Next, we investigated fined by selection to enhance coloration, and we propose both a po- the effect of increasing lutein concentration on colorimetric vari- tential fitness detriment that results from those modifications and a ables obtained during reflectance spectroscopy of collar feathers potential compensatory mechanism. These findings reveal the intri- only (to remove the effect of barbule morphology present in belly cate interplay of sexual selection, hybridization, and introgression feathers) using a mixed effects model that incorporated specimen involved in shaping the ornamentation and behavior of these birds identity to account for non- independence of data within species/ and hint at additional levels of complexity yet to be explored. hybrid form. These variables were saturation in the yellow wave- lengths of light (S1yellow), UV reflectance saturation (S1UV), ca- rotenoid chroma (S9), mean sample brightness (B2), and sample MATERIALS AND METHODS midwavelength hue (H3). Next, we compared feather brightness (as Feather color and pigment analysis measured by reflectance spectroscopy) of collar feathers of each of We examined the collar and belly feathers of four species in the genus the five taxa. Next, we modeled the relationship between barb Manacus (M. manacus, M. candei, M. vitellinus, and M. aurantiacus), width and lutein concentration and between barb width and the as well as from M. candei x M. vitellinus hybrids. Five feathers were colorimetric variables described above. Last, we compared the con- plucked from each of two plumage patches, the belly and collar, of six centration of lutein and the colorimetric variables above between museum specimens of each taxon from the US National Museum of yellow belly feathers of M. candei and yellow collar feathers of Natural History (NMNH). Specimen information is available in table M. vitellinus and hybrids. S1. We used a combination of light microscopy, reflectance spectros- From each model analysis, we report the estimated effect size (β: copy, and HPLC to characterize the pigments in feathers and the the difference between groups or the slope of the relationship be- mechanisms of color production. We sought to determine which ca- tween two variables), the confidence interval around the estimated rotenoids were present and their concentrations and to visually de- effect size, and the associated P value (Table 1 and tables S2 and S3). termine the presence and location of melanin in feather morphology. All supplemental feather figures and graphs were built in R using the Each feather was examined under reflected light before and after ggpubr v.0.2.1 (91) and cowplot v.1.0.0 (92) packages. carotenoid extraction. The presence or absence of barbules along the barbs and the presence of pigmentation (carotenoids and/or melani- Generation of genome resequencing data and zation) in the barbs and barbules were noted. Pictures were taken reference- based scaffolding using an AmScope Microscope Digital Camera at a magnification of Samples used for genome resequencing were collected in northwest ×10 and ×40. or central Panama. All birds were collected in 1980–1990s and pre- When performing UV-v isible spectroscopy, five feathers were served as voucher specimens at the US Museum of Natural History overlapped and taped down flat, ventral side up, onto a black card- as round or flat study skins or skeletal specimens. All samples were board background. This was done to simulate how the feathers may previously studied by Brumfield et al. (23), so we follow their popu- have been placed on the bird during life and to eliminate any back- lation/site designations for consistency. We sequenced 12 M. candei, ground noise when sampling the feathers. Three measurements 4 M. vitellinus, and 10 hybrids from a population with a predomi- were taken in a dark room at a 90° angle of incidence from each set nantly candei genomic background but some introgression from of feathers using an Ocean Optics USB4000 spectrophotometer and vitellinus (table S5). Two M. candei individuals represented paren- an Ocean Optics PX-2 pulsed xenon light source relative to a white tals from site 2 near the border with Costa Rica and 10 were from reflectance standard (Labsphere Inc., North Sutton, NH, USA). site 3 on the west bank of the Río Changuinola (Fig. 1C). The 10 Spectra were measured using OceanView v1.67 (Ocean Insight, FL, hybrids were from site 4 on the east bank of the Río Changuinola USA). Spectral measurements were taken from the left, center, and (Fig. 1C). The four parental M. vitellinus were from site 12, ~250 km right of the distal end of the overlapped group of feathers. Capturing east of the Rio Changuinola, far from the plumage or genomic con- multiple measurements from each plumage patch helps account for tact zones (Fig. 1C). Collar and belly plumage color of the birds any variation in color within or between feathers. The resulting from sites 3 and 4 was confirmed by examining the museum speci- spectral data were analyzed and graphed using the pavo package v. mens. DNA, extracted through a phenol-c hloroform and ethanol 2.4.0 (33) in R v. 4.3.3 (89). precipitation method, was shipped to the McDonnell Genome Insti- For carotenoid analysis, we sampled two feathers in each feather tute, Washington University for library preparation and DNA se- patch from each species and hybrids. All feathers were cut approxi- quencing across three lanes of the Illumina HiSeq X platform (2 × mately one-t hird of the way up from the base, removing the bottom 150 bp paired- end reads). Raw sequencing reads are available from portion of the feather that was heavily melanized and saving the top National Center for Biotechnology Information (NCBI) at BioPro- two- thirds of the feather. Carotenoids were extracted from cut ject PRJNA1129471. feather sections and analyzed using HPLC. See Supplementary Text The average number of raw FASTQ reads per individual was for detailed methods. 97.6 million (range: 64.4 to 123.6 million). We used Cutadapt im- We performed statistical analyses of feather data in R using lin- plemented in Trim Galore v0.6.4 (https://github.com/FelixKrueger/ ear models corrected for multiple comparisons using the emmeans TrimGalore) to conduct quality and adapter trimming (default package v.1.4.5 (90) and linear mixed effects models. First, we com- parameters). The average read length and number of reads per sam- pared the concentration of lutein (square root transformed to ple after quality trimming were 145.4 bp and 96.3 million (range: achieve normality) in the belly or collar feathers among the five 63.7 to 121.6 million), respectively. We used bowtie2 v2.3.5 (93) to species and hybrid included in this study. Next, we compared the map trimmed and filtered reads to the annotated reference genome concentration of lutein in the belly versus the collar feathers of the of M. vitellinus (NCBI RefSeq accession GCF_001715985.3). At three taxa which had yellow pigmentation in both feather patches the time of read mapping, this reference genome was the only one Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 12 of 20 S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e available for the genus. Samtools v1.9 was then used to sort and Pure and Applied Chemistry ambiguity codes. Before NeighborNet index the bam files (94). This reference genome (scaffold N50 = network construction, we used the HKY85 substitution model, a 17.9 Mbp, contig N50 = 290.6 kbp) was generated with a combina- model appropriate for closely related species, to calculate genetic tion of Illumina HiSeq short reads and Pacific Biosystem long reads distances among individuals. To simplify the network, we filtered and assembled with MaSuRCA v3.1.1 (95). edges by setting maximum dimension to 4 (i.e., removing splits that Across all samples, the average overall botwtie2 alignment rate caused boxes in the network to have dimensions > 4). was 94.6% (range: 91.5 to 96.5%). Next, we followed recommended To conduct principal components analysis, we first used ANGSD best practices and used GATK v3.8.1 to conduct variant calling (96). v0.918 to generate filtered genotype posterior probabilities (com- Briefly, we used the HaplotypeCaller tool to conduct calling of SNPs mand: - doMaf 1 - doMajorMinor 1 - SNP_pval 1e- 6 - uniqueOnly and insertions/deletions (indels) simultaneously via local de novo 1 - remove_bads 1 - only_propoer_pairs 1 - C 50 – baq 1 -m inMapQ assembly of haplotypes. This generated per-s ample intermediate 20 - minQ 13 - minMaf 0.05 - minind 21 - GL 1 - doGeno 32 -d oPost GVCF files, which were then combined with the GenotypeGVCFs 1) (101). Only SNPs with data in at least 21 of 26 individuals and a tool to produce the final raw VCF file. We then selected for SNPs minimum allele frequency of 5% were kept. This created a genotype and used the GATK tool VariantFiltration to conduct both SNP- posterior probabilities output containing 13.1 million SNPs. Follow- and genotype-l evel filtering. We used the expressions “QualByDepth ing this, we used ngsCovar v1 to randomly subsample 1 million < 2.0 || FisherStrand > 60.0 || StrandOddsRatio > 3.0 || RMSMap- SNPs to reduce linkage disequilibrium and to calculate the covari- pingQuality < 40.0 || ReadPosRankSum < −3.0 || ReadPosRank- ance matrix among individuals (102). Eigenvector decomposition Sum > 3.0 || BaseQRankSum < −3.0 || BaseQRankSum > 3.0 || was then conducted in the R statistical package to generate the prin- MQRankSum < −3.0 || MQRankSum > 3.0” and “DP < 5.0 || GQ < cipal components (89). 13.0” to remove low-q uality or spurious SNPs and genotypes, re- spectively, during the filtering process. Last, we retained SNPs that Modeling population divergence and gene flow contained genotype information in at least 21 individuals (of the Using the genome resequencing data and an allele frequency spec- total 26), producing a filtered VCF file containing 19.8 million SNPs. trum approach, we estimated population divergence histories be- At the individual level, the proportion of SNPs with missing geno- tween M. vitellinus (N = 4), M. candei (N = 10, site 3 only), and type information ranged from 2.9 to 22.5% (table S3). hybrids (N = 10, site 4) and between the latter two populations only Because a downstream analysis (ABBA- BABA analysis) required (table S3). For the three-p opulation analysis, we removed any SNPs the use of an outgroup taxon to determine derived- ancestral allelic that have missing data across the three populations, creating a full states, we downloaded Illumina HiSeq 2000 FASTQ reads from a data matrix. Following this, we kept only SNPs that were at least whole-g enome sequencing project of Pipra filicauda (NCBI Run Ac- 1000 bp from each other to reduce linkage disequilibrium among cession: SRR6885520; 105.5 × 1012 bp, 349.2 × 106 150 bp paired- them. This resulted in a folded three- dimensional (3D) site frequen- end reads). Following the above workflow, we mapped P. filicauda cy spectrum that contained 608,905 segregating sites. Next, we reads to the M. vitellinus reference genome and generated a separate modified an optimization routine (moments_Run_Optimizations. VCF file with outgroup data. This VCF contained 44.5 million SNPs py) in moments_pipeline (github.com/dportik/moments_pipeline) after SNP-l evel, genotype- level, and missingness filtering. to evaluate the fit of five candidate population divergence models To place M. vitellinus scaffolds in their approximate chromosomal (table S6 and fig. S6) to the data. The population branching history positions, we conducted reference- guided chromosome-l evel of all three-p opulation models was (M. vitellinus, (M. candei, hy- scaffolding, using the genome of C. lanceolata (lance-t ailed mana- brids)), and the tested models varied with respective to whether mi- kin) as the reference (NCBI RefSeq accession GCF_009829145.1). gration was allowed between pairs of populations. The lower and This high- quality genome assembly consisted of full chromosomes, upper bounds of population size (0, 40), time of split (0.0001, 10), including 33 autosomal macro- and microchromosomes, and two and migration rate (0, 50) parameters were kept the same when test- sex chromosomes (W and Z). We used RaGOO v1 (97) with default ing all models. These parameter bounds have been verified to be parameters to conduct scaffolding. As expected, given the high level appropriate through multiple preliminary moments runs. of synteny among bird species, the scaffolding process resulted in After the preliminary runs were completed, for each demographic high confidence in terms of placing and orienting each M. vitellinus model being tested, we first conducted four to five independent scaffold against C. lanceolata chromosomes (average positional moments_pipeline optimization runs, with each run containing four confidence score = 77.8%, SD = 27.6%; average orientation confi- sequential rounds and each round containing 10 (round 1) or 20 dence score = 98.2%, SD = 7.3%). (rounds 2 to 4) independent replicate moments runs. While starting parameters of the first- round replicates were random, parameters Principal components and phylogenetic network analysis from the best scoring replicate after each round was completed were To reveal relative divergence among populations and evolutionary re- used to initiate perturbed starting parameters for replicates of the lationships among individuals, we conducted phylogenetic network subsequent round. The levels of perturbation of starting values (i.e., analysis using the NeighborNet algorithm implemented in SplitsTree “fold” value) for the four sequential rounds were 3, 2, 1, and 1, re- v5 (98, 99). We took the above filtered VCF file containing 19.8 mil- spectively. lion SNPs and removed SNPs with minor allele frequency of less than When parameter optimization of each model was completed, we 5%. This was followed by thinning of SNPs so that no two sites were inspected log likelihood values across runs to ensure convergence. within 10,000 bp of each other. This produced a final VCF file con- After all candidate models were evaluated with moments_pipeline, taining 103,885 SNPs. Using PGDSpider v2.1.1.5 (100), we converted we then used the Akaike information criterion to select the best de- the reduced VCF file into a FASTA file (one sequence per individual) mographic model. Using this best demographic model, we conducted while collapsing heterozygous genotypes into International Union of one final, long moments run (maximum number of iterations = Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 13 of 20 S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e 100), using demographic parameters from the best round 4 moments_ Because the plumage differences between M. candei and hybrids pipeline replicate as starting values. To estimate parameter uncer- are likely due to introgression of these traits from M. vitellinus into tainties, we first generated 100 bootstrapped folded spectra and used the hybrid population (13, 23), we calculated introgression statistics the Godambe Information matrix, an appropriate approach for com- to identify genomic regions in hybrids that originated in M. vitellinus. posite likelihoods and SNPs that are not fully unlinked, to calculate We predicted that regions underlying the plumage color differences 95% confidence intervals (103). would show signatures of both differentiation and introgression. We The same optimization and modeling approaches were applied to calculated fd, which is a modified version of the ABBA- BABA D sta- a 2D site frequency spectrum between M. candei from site 3 and tistic intended to be applied to population samples in small genomic hybrids from site 4 (number of segregating sites = 399,572). We windows using the ABBABABAwindows.py script (github.com/ tested four candidate models that varied with respect to whether or simonhmartin/genomics_general) (43). fd is an estimator of the pro- when migration occurred between the two populations (table S6). portion of a genomic region that is shared due to introgression be- To convert substitution rate–scaled parameter estimates to absolute tween populations P2 and P3 in a set of populations with the units, we used a genome-w ide substitution rate of 3.3 × 10−9 substi- relationship (((P1, P2), P3), O). fd is a more sensitive and less biased tutions/year and a generation time of 2.5 years. statistic than D when applied to small genomic regions (43), and unlike FST and Da, fd is not sensitive to local recombination rate Genomic window-b ased divergence and variation since it depends on the ratio of various allelic states across introgression metrics a tree rather than measures of genetic diversity (104). In our calcula- Because of our sampling scheme along the Río Changuinola—very tions, we used P1 = M. candei (site 3, N = 10), P2 = hybrids (site 4, similar populations genetically but with a major color difference— N = 10), and P3 = M. vitellinus (site 12, N = 4), and the outgroup we expected genomic regions underlying color to show strong sig- was P. filicauda. We omitted windows with fewer than 100 biallelic nals of both differentiation across the river and introgression from SNPs used to calculate fd, as recommended by the author of ABBA- M. vitellinus. To locate these regions, we calculated differentiation BABAwindows.py. and introgression statistics in sliding windows along genomic scaf- In addition to fd, we estimated the RND of the split between folds. Because any single population genetic statistic may be affected site 3 M. candei and site 4 hybrids to the split between hybrids and by confounding factors and because we expected a region underly- M. vitellinus (site 12) by calculating sites 3 to 4 DXY ÷ sites 4 to 12 ing color to be robust to such issues, we used multiple analysis DXY (44). Assuming a consistent substitution rate across lineages, methods to assess both differentiation and introgression (40,  41). values greater than one represent genomic windows where the di- We calculated these statistics in non- overlapping windows of 10, 25, vergence between hybrids and M. vitellinus is more recent than and 50 kb (fig. S35). Results were similar between window sizes, and the hybrid/M. candei divergence. Because the genomes of hybrids the smallest size allowed us to narrow down regions of interest more are predominantly derived from M. candei (23,  29,  36), these precisely. For that reason, we used 10-k b windows in downstream windows should be very rare in the genome and indicate intro- analyses and present those results in Fig. 3. We omitted all scaffolds gression from vitellinus. smaller than 10 kb, which accounted for 0.25% of the total genome A total of 266 autosomal (of 91,913) and 26 Z chromosome (of length. We identified focal regions potentially underlying the plum- 7079) 10-k b windows were considered 99.9th percentile outliers for age color difference by selecting windows that fall above the 99.9th at least one of the four summary statistics. A substantial proportion percentile of the differentiation and introgression statistics (de- of these windows were flagged as outliers by two or more summary scribed below). The avian Z chromosome has a smaller effective statistics (autosomes: 32.3%; Z chromosome: 23.1%) (fig. S9). Among population size than autosomes (macro- and microchromosomes), the windows flagged by two summary statistics (n = 75), the major- and as a result, differentiation and introgression statistics behaved ity were flagged by either the two introgression statistics, fd and RND differently on Z than on autosomes (fig. S36). Therefore, we calcu- (autosomes: 37 of 75; Z chromosome: 4 of 6), or the two differentia- lated Z chromosome 99.9th percentile thresholds separately from tion statistics, FST and Da (autosomes: 36 of 75; Z chromosome: 1 of autosome thresholds. 6). The fact that outlier windows for introgression statistics did not We predicted that genomic regions underlying the plumage color often overlap with differentiation statistics outliers suggests that high differences between M. vitellinus and M. candei would show substan- differentiation (between hybrids and M. candei in site 3) of some tial genetic differentiation between birds from sites 3 and 4 on either windows might have arisen from processes other than introgression side of the Río Changuinola (Fig. 1C and table S3). We used the pop- (from M. vitellinus into hybrids), including local selection or drift in genWindows.py script (github.com/simonhmartin/genomics_general) one of the two M. candei populations after their lineages diverged or to calculate the nucleotide diversity, π, of each population, the fixa- increased FST or Da due to low π and/or recombination rates (104). tion index, FST, and sequence divergence, DXY (42), between the Nevertheless, a number of autosomal outlier windows were supported two populations. FST is the proportion of total genetic variation by three or all four summary statistics, indicating that these were strong accounted for by differences between populations, while DXY is the candidates for windows containing loci that underlie the plumage average of the pairwise sequence divergence between population introgression (fig. S9). samples. The method we used calculates variant site measures of DXY To facilitate the inspection of gene content of genomic regions and π, so we calculated invariant sites in each window and adjusted that received strong support from our window- based analyses, we DXY and π values accordingly (see get.popgen.data.R in Supplemen- selected focal windows based on the following approach. We consid- tary Code for details). Because DXY is heavily influenced by nucleo- ered an autosomal window (window size = 10 kb) a focal window if tide diversity, π, (Pearson’s r = 0.984; fig. S37A), we instead used Da it was marked as a 99.9th percentile outlier by at least three of the (42) for downstream analyses, which we calculated by subtracting four summary statistics (num_outliers ≥ 3). In addition, we used the the mean π of the two populations from DXY (fig. S37B). rollmean function of the R package zoo (105) to calculate centered Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 14 of 20 S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e moving average of num_outliers; for each window, we calculated the Cline analysis average of its num_outliers and those of the two windows that were We sought to increase the geographic resolution of our introgres- immediately up or downstream (MAV3). Next, an autosomal win- sion results by fitting geographic clines to phenotypic traits SNPs dow was considered a focal window if its num_outliers ≥ 3 or MAV3 genotyped using RADseq from additional sampling sites in the ≥ 1. The use of moving averages as a criterion for identifying focal Manacus distribution. We sequenced 152 individuals from eight windows allowed us to find genomic regions of interest where outlier sites (Fig. 1C, sites 2 to 10) to capture variation in hybrid as well as windows for the various introgression/differentiation statistics were parental populations, following the sampling transect from (23). not exactly coincident but could be found in adjacent windows. Be- The locations and sample sizes of each site are available in table S14. cause no Z chromosome windows had num_outliers ≥ 3 (fig. S9B), For phenotypic clines, we measured three male traits (beard length, the following criteria were used to identify focal windows: num_ collar color, and belly color) using the whole museum skins of male outliers ≥ 2 or MAV3 ≥ 1. Manacus manakins that were measured by (23) and deposited at Using these criteria, we found a total of 74 and 8 focal windows NMNH. Detailed methods are described in (29). For the genetic in autosomes and the Z chromosome, respectively. We further clines, we used DNA samples archived at NMNH from individuals organized the focal windows into 24 genomic neighborhoods captured by Parsons et al. (13) and Brumfield et al. (23) and then where focal windows were contiguous or no more than 1000 kb sequenced in (29). Detailed methods on DNA extraction, sequenc- (100- window length) from each other (table S10 and figs. S10 to ing, SNP calling, filtering, including cline width credibility interval S33). These genomic neighborhoods were found only in six auto- filters, and genetic and phenotypic cline fitting are described in (29) somes and the Z chromosome. They each spanned 10 to 960 kb, following the methods for the “historical” dataset. Our final SNP containing zero to eight ab initio annotated protein coding genes. cline dataset consisted of 4750 geographic clines with precise cline Three such genomic neighborhoods stood out as being highly centers. Raw reads are available at NCBI BioProject PRJNA893627. prominent in that they each contained one (genomic neighbor- hood 11) or two (genomic neighborhoods 9 and 19) 10-k b win- BCO2 expression dows that were outliers for all four summary statistics (table S10). We explored Manacus lineage–specific gene expression of BCO2 A list of genes, and their functions, found in strongly supported using gene expression data from scapulohumeralis and pectoralis windows can be found in table S11. muscles in six species of manakins and a suboscine outgroup from (67). Muscle reads were mapped to the P. filicauda genome Pathway and gene ontology term (GCF_003945595.1) using the splice aware- aligner STAR v 2.7 overrepresentation analysis (108) and counted against gene features using featureCounts v2.0.1 To test whether the outlier-m arked regions of the genome were en- (109). Data on BCO2 expression in cerebellum came from three riched with genes from defined molecular pathways or functions, manakin species and a suboscine outgroup (table S16). Cerebel- we used the Protein Analysis Through Evolutionary Relationships lums were collected from breeding males of each species, and RNA (PANTHER) classification system (106). We filtered the M. vitellinus was extracted using a Qiagen RNeasy kit and sequenced in 2 × 100 bp reference genome for isomers annotated under multiple protein IDs reads on an Illumina NextSeq 500. Reads are available at NCBI and consolidated to a single isomer per gene. The filtered M. vitellinus BioProject PRJNA321179. Cerebellum reads counts were estimated reference genome amino acid sequences for all proteins were que- using kallisto, using 30 sample bootstraps and transcripts available ried against the PANTHER hidden Markov model library using the for M. vitellinus (GCF_001715985.3), and converted to counts us- hmmscan option to return the best hit (107). Once functional clas- ing R package tximport. We compared lineage- specific gene expres- sifications were obtained for as many of the proteins in the genome sion in these muscle and brain tissues using PhyDGET. This method as possible, genes (again, filtered down to single isomer representa- assigns Bayes factors (>|1.5|) to infer significant changes in gene tion) that appeared within a 10-k b outlier window supported by any expression along specific branches. of the four differentiation or introgression statistics or 20 kb up- or downstream of the outlier window were defined as the outlier subset Historical BCO2 gene flow in Manacus from the full genome complement of proteins. We chose an addi- We prepared whole-g enome sequencing libraries for nine addi- tional area that was +20 kb of the outlier window because it could tional individuals to encompass the major lineages within Manacus, contain a distal cis- regulatory element (e.g., enhancer and insulator) including three M. aurantiacus, three M. manacus from east of the with its associated gene. Andes (cis-A ndean), and three M. manacus from west of the Andes We used the PANTHER Gene List Analysis tool to perform a (trans-A ndean) (table S15). The libraries were sequenced by Novogene PANTHER pathways statistical overrepresentation test; in short, Corp. on one Illumina HiSeq lane, producing 217 million total the functional classifications for the entire genome complement reads averaging 23.8 million reads per sample. Raw sequencing served as a basis for expected composition of an outlier subset. reads are available from NCBI at BioProject PRJNA1129471. Reads Without pathway enrichment, one would expect to observe the were trimmed using Cutadapt implemented in Trim Galore and same proportion of genes from a given pathway in a subset as are aligned to the M. vitellinus reference genome using bwa mem present in the entire genome complement. Statistical deviation (110). We marked duplicates using picard MarkDuplicates. We fol- from expected values in our observed outlier gene set was assessed lowed the GATK best practices workflow for variant calling, in- with a Fisher’s exact test including a false discover rate correction cluding producing and filtering a VCF file, and then used the for multiple tests. We also used the PANTHER gene list analysis filtered VCF to apply a recalibration table to the original bam files tool to assess whether any GO- Slim (Gene Ontology Slim) term using GATK v3.8.1.0 (96). The terms we used for removing low- categories (molecular function) were overrepresented in the same quality or spurious SNPs and genotypes were “QD < 2.0 || FS > manner as the pathways. 60.0 || SOR > 3.0 || MQ < 40.0 || ReadPosRankSum < −3.0 || Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 15 of 20 S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e ReadPosRankSum > 3.0 || BaseQRankSum < −3.0 || BaseQRank- For a finer-s cale view of introgression evidence in the BCO2 re- Sum > 3.0 || MQRankSum < −3.0 || MQRankSum > 3.0” and “DP < gion, we calculated fd on chromosome 24 using 2- kb windows and a 4.0 || DP > 45.0 || GQ < 3.0,” followed by a filter to remove indels and 10- site minimum cutoff. Using the all-s pecies VCF file described sites missing in more than five of nine individuals in VCFtools above, we also made trees in 100- bp windows across the chromo- v0.1.16 (111). some. For this analysis, we subset the main VCF into a window and We selected two control genes on chromosome 24 similar in then obtained the consensus sequence for each of four species size to BCO2, more than 1 Mb away, and with no known connec- (M. vitellinus, trans- Andean M. manacus, parental M. candei, and tion to coloration, DCPS and TBCEL. Our rationale was to deter- M. aurantiacus) using bcftools v1.18 (114), aligned the sequences mine whether any introgression signal we detected at BCO2 was with Mafft v7.520 (112), and built a tree using fasttree v2.1.11 (115). independent of its genomic context. We identified genomic re- Next, we identified all trees in which M. aurantiacus and M. vitellinus gions to use for building gene trees by anchoring the ends of the were sister to one another and calculated the frequency of such trees regions of interest in exons of flanking genes to aid sequence in consecutive bins of 20 and 100 windows across the chromosome, alignment. We used a much shorter sequence for the BCO2 gene corresponding to 2- and 10- kb bin sizes, respectively. tree (48 kb, from 26 kb upstream to 5 kb downstream) than the control trees (306 and 384 kb, respectively) because of the likeli- Feather water shedding in Manacus hood that recombination has reduced signals of shared ancestry The waterproofing of a bird’s feathers is determined by two forces: that may have been present if introgression occurred at BCO2. water repellency and water penetration resistance (34). Water repel- Using bam files from each individual, including two M. candei lency is described as the ability for a feather to cause water to bead and four M. vitellinus whole-g enome resequencing samples used in and shed from its surface, and water penetration resistance can be the above analyses, we produced a single consensus sequence per defined as the force required to push water through the barbs of individual using Samtools consensus (94) and then aligned the se- overlapping feathers (34). These forces are in partial physical oppo- quences using Mafft (112). To facilitate tree rooting, we included the sition to each other, as they are both determined by the diameter of corresponding genomic regions from the reference genomes of a feather’s barbs and the spacing between adjacent barbs (34). P. filicauda (accession ASM394559v2) and C. lanceolata (Accession Mathematical estimations of water repellency and water penetra- bChiLan1.pri) in the alignments. tion resistance have been developed on the basis of the similarity of We used IQ-T ree v2 (113) to produce maximum likelihood gene morphology of feather structure to adjacent and overlapping cylin- trees from each alignment with 10,000 bootstrap replicates, and in drical surfaces (34). These estimates have been applied to multiple each case, we allowed IQ-T ree to select the best model of sequence bird species that inhabit aquatic and arboreal habitats (35, 116). Ter- evolution using Bayesian information criterion. The resulting trees restrial and arboreal birds are believed to benefit the most from were rooted at the midpoints. maximizing water repellency, while aquatic birds benefit most from Next, to place the evidence of gene flow at BCO2 in the context penetration resistance (86, 116). Here, we use these estimates to ex- of genome-w ide signals of introgression between M. aurantiacus amine the collar feathers of the four species in the genus Manacus and M. vitellinus, we calculated ABBA- BABA D statistics. To pro- (M. manacus, M. candei, M. vitellinus, and M. aurantiacus) and the duce a VCF file for this analysis containing all species, we con- M. candei x M. vitellinus hybrids that attained vitellinus plumage ducted an additional round of SNP calling using bam files from the coloration through genetic introgression. We compared feathers initial M. vitellinus, candei, and hybrid sequencing as well as bam from each species and the hybrids to look for differences in feather files from the subsequent M. manacus and M. aurantiacus sequenc- structure that could lead to fitness-r elated consequences in terms of ing as inputs for GATK HaplotypeCaller and CombineGVCFs and waterproofing ability. We used light microscopy, photographical then filtered using the following terms: “QD < 2.0 || FS > 60.0 || measurement, and wettability modeling to estimate the ability of SOR > 3.0 || MQ < 40.0 || ReadPosRankSum < −3.0 || ReadPos- feathers to repel water and resist water penetration. RankSum > 3.0 || BaseQRankSum < −3.0 || BaseQRankSum > 3.0 For microstructural analyses, we used the same collar feathers || MQRankSum < −3.0 || MQRankSum > 3.0” and “DP < 4.0 || GQ < from the color analyses described above. From the collar feathers, we 3.0.” Last, we filtered out indels and any sites missing in more than sampled the most intact feather from each of the six individuals per 10 of 15 individuals. taxon. The feathers chosen were not used in any color analysis be- ABBA-B ABA statistics require a four- species tree, and the cal- yond light microscopy—they were not subjected to HPLC or spec- culation of the fd statistic makes the assumption that, in the tree troscopy before our measurement for this analysis. Each feather was structure (((P1,P2),P3,)P4), gene flow is directional from P3 to P2 examined under reflected light on a glass calibration slide with an (43). Gene flow in the opposite direction will still result in elevated etched scale accurate to 0.1 mm. Feathers were placed on this scale, fd, but the signal will be reduced. We used the genomics_general and pictures were taken using an AmScope Microscope Digital Cam- package (https://github.com/simonhmartin/genomics_general) to era at a magnification of ×40. Photographs were imported into FIJI calculate fd across the genome in 10- kb windows (43) with trans- software, and morphological measurements were recorded using the Andean M. manacus as P1, M. vitellinus as P2, M. aurantiacus as calibration slide scale as reference. We recorded feather measure- P3, and cis-A ndean M. manacus as P4. We changed values of fd to ments in the middle of the feather to avoid any fraying of the feather zero in windows where D was negative, where fd was negative, or edges influencing the data according to the recommendation in (34). where fd was greater than one, as recommended by the author. Re- We recorded two main metrics in accordance with the methods of moving windows with fewer than 100 sites used to calculate fd is (34, 35): feather barb width and the spacing between feather barbs. recommended by the author to reduce stochastic variation in fd (fig. We also recorded barb length from the base of the barb to the distal S38A), but this cutoff was not feasible for our dataset (fig. S38B). tip. From each feather, we recorded these measurements using five Instead, we used a 50- site cutoff. barbs and, therefore, retained feather ID in all statistical models to Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 16 of 20 S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e control for pseudoreplication and non-i ndependence of barb mea- 6. S. A. taylor, e. L. Larson, insights from genomes into the evolutionary importance and surements within each feather. prevalence of hybridization in nature. Nat. Ecol. Evol. 3, 170–177 (2019). We calculated water repellency using Eq. 1 below, as de- 7. d. A. Marques, J. i. Meier, O. Seehausen, A combinatorial view on speciation and adaptive radiation. Trends Ecol. Evol. 34, 531–544 (2019). scribed in (34) 8. B. M. Moran, c. Payne, Q. Langdon, d. L. Powell, Y. Brandvain, M. Schumer, the genomic consequences of hybridization. eLife 10, e69016 (2021). (R +D)∕R (1) 9. d. Berner, W. Salzburger, the genomics of organismal diversification illuminated by where R is equal to the radius of each barb and D is equal to half the adaptive radiations. Trends Genet. 31, 491–499 (2015). 10. K. J. Liu, e. Steinberg, A. Yozzo, Y. Song, M. h. Kohn, L. Nakhleh, interspecific introgressive distance between adjacent barbs. Lower values of the repellency esti- origin of genomic diversity in the house mouse. Proc. Natl. Acad. Sci. U.S.A. 112, 196–201 mate indicate reduced water repellency and higher values represent (2015). better water repellency. To estimate water penetration resistance, we 11. A. Suarez-G onzalez, c. A. hefer, c. christe, O. corea, c. Lexer, Q. c. B. cronk, c. J. douglas, used Eq. 2 below, also as described in (34) Genomic and functional approaches reveal a case of adaptive introgression from Populus balsamifera (balsam poplar) in P. trichocarpa (black cottonwood). Mol. Ecol. 25, γ 2427–2442 (2016). P = (2) 12. M. R. Jones, L. S. Mills, P. c. Alves, c. M. callahan, J. M. Alves, d. J. R. Lafferty, F. M. Jiggins, R� J. d. Jensen, J. Melo- Ferreira, J. M. Good, Adaptive introgression underlies polymorphic where P is the pressure required to force water between feather seasonal camouflage in snowshoe hares. Science 360, 1355–1358 (2018). barbs in grams per square centimeter, γ is the surface tension of wa- 13. t. J. Parsons, S. L. Olson, M. J. Braun, Unidirectional spread of secondary sexual plumage traits across an avian hybrid zone. Science 260, 1643–1646 (1993). ter, and R′ is the principal radius of curvature of a feather barb, cal- 14. d. t. Baldassarre, t. A. White, J. Karubian, M. S. Webster, Genomic and morphological culated using Eq. 3 below analysis of a semipermeable avian hybrid zone suggests asymmetrical introgression of a sexual signal. Evolution 68, 2644–2657 (2014). √ { }2 R + D 15. G. M. While, S. Michaelides, R. J. P. heathcote, h. e. A. MacGregor, N. Zajac, J. Beninde, � = θ + − 2R cos sin θ (3) P. carazo, G. Pérez i de Lanuza, R. Sacchi, M. A. L. Zuffi, t. horváthová, B. Fresnillo, R U. Schulte, M. veith, A. hochkirch, t. Uller, Sexual selection drives asymmetric where θ is the contact of angle of water with the feather surface. As introgression in wall lizards. Ecol. Lett. 18, 1366–1375 (2015). in Rijke et al. (34, 35), we set γ to be equal to the surface tension of 16. S. e. Lipshutz, J. i. Meier, G. e. derryberry, M. J. Miller, O. Seehausen, e. P. derryberry, water at 20°C (72.86 dynes/cm) and θ to be equal to 90°. differential introgression of a female competitive trait in a hybrid zone between sex-r ole We calculated differences among taxa in barb width, barb spac- reversed species. Evolution 73, 188–201 (2019). 17. O. Seehausen, Y. terai, i. S. Magalhaes, K. L. carleton, h. d. J. Mrosso, R. Miyagi, ing, barb length, water repellency [(R + D)/R], and water penetra- i. van der Sluijs, M. v. Schneider, M. e. Maan, h. tachida, h. imai, N. Okada, Speciation tion resistance using linear mixed effects models in R using the lme4 through sensory drive in cichlid fish. Nature 455, 620–626 (2008). package (117) and corrected for multiple comparisons using the 18. N. Seddon, c. A. Botero, J. A. tobias, P. O. dunn, h. e. A. MacGregor, d. R. Rubenstein, emmeans package v.1.8.5 (90). In this model, we used taxon as a J. A. c. Uy, J. t. Weir, L. A. Whittingham, R. J. Safran, Sexual selection accelerates fixed effect and feather ID as a random effect to control for the fact signal evolution during speciation in birds. Proc. R. Soc. B Biol. Sci. 280, 20131065 (2013). that five barb measurements came from each feather. We plotted our 19. M. J. Ryan, M. e. cummings, Perceptual biases and mate choice. Annu. Rev. Ecol. Evol. Syst. results using the cowplot v.1.1.1 (92), ggpubr v.0.6.0 (91), and patch- 44, 437–459 (2013). work v1.1.2 (118) packages. From each model analysis, we report 20. G. e. hill, Plumage coloration is a sexually selected indicator of male quality. Nature 350, the estimated effect size (β: the difference between groups), the SE, 337–339 (1991). 21. J. d. Ligon, The Evolution of Avian Breeding Systems (Oxford Univ. Press, 1999). and the confidence interval around the estimated effect size and the 22. M. Â. Marini, S. J. hackett, A multifaceted approach to the characterization of an associated P value in table S4. intergeneric hybrid manakin (Pipridae) from Brazil. Auk 119, 1114–1120 (2002). 23. R. t. Brumfield, R. W. Jernigan, d. B. Mcdonald, M. J. Braun, evolutionary implications of divergent clines in an avian (Manacus: Aves) hybrid zone. Evolution 55, 2070–2087 Supplementary Materials (2001). 24. A. O. Barrera-G uzmán, A. Aleixo, M. d. Shawkey, J. t. Weir, hybrid speciation leads to The PDF file includes: novel male secondary sexual ornamentation of an Amazonian bird. Proc. Natl. Acad. Sci. Supplementary text U.S.A. 115, e218–e225 (2018). Figs. S1 to S40 25. A. e. Moncrieff, B. c. Faircloth, R. t. Brumfield, Systematics of Lepidothrix manakins (Aves: tables S1 to S16 Passeriformes: Pipridae) using RAdcap markers. Mol. Phylogenet. Evol. 173, 107525 Legend for movie S1 (2022). Legends for data S1 and S2 26. W. R. Lindsay, J. t. houck, c. e. Giuliano, L. B. day, Acrobatic courtship display coevolves References with brain size in manakins (Pipridae). Brain Behav. Evol. 85, 29–36 (2015). 27. A. c. Stein, J. A. c. Uy, Plumage brightness predicts male mating success in the lekking Other Supplementary Material for this manuscript includes the following: golden-c ollared manakin, Manacus vitellinus. Behav. Ecol. 17, 41–47 (2006). Movie S1 28. A. c. Stein, J. A. c. Uy, Unidirectional introgression of a sexually selected trait across an data S1 and S2 avian hybrid zone: A role for female choice? Evolution 60, 1476–1485 (2006). 29. K. M. Long, A. G. Rivera- colón, K. F. P. Bennett, J. M. catchen, M. J. Braun, J. d. Brawn, REFERENCES AND NOTES Ongoing introgression of a secondary sexual plumage trait in a stable avian hybrid zone. 1. M. J. Ryan, A Taste for the Beautiful: The Evolution of Attraction (Princeton Univ. Press, Evolution 78, 1539–1553 (2024). 2018). 30. d. B. Mcdonald, R. P. clay, R. t. Brumfield, M. J. Braun, Sexual selection on plumage and 2. M. Andersson, Sexual Selection (Princeton Univ. Press, 1994). behavior in an avian hybrid zone: experimental tests of male- male interactions. Evolution 3. W. J. Swanson, v. d. vacquier, the rapid evolution of reproductive proteins. Nat. Rev. 55, 1443–1451 (2001). Genet. 3, 137–144 (2002). 31. J. Barske, M. J. Fuxjager, c. ciofi, c. Natali, B. A. Schlinger, t. Billo, L. Fusani, Beyond 4. G. S. Wilkinson, F. Breden, J. e. Mank, M. G. Ritchie, A. d. higginson, J. Radwan, J. Jaquiery, plumage: Acrobatic courtship displays show intermediate patterns in manakin hybrids. W. Salzburger, e. Arriero, S. M. Barribeau, P. c. Phillips, S. c. P. Renn, L. Rowe, the locus of Anim. Behav. 198, 195–205 (2023). sexual selection: Moving sexual selection studies into the post-g enomics era. J. Evol. Biol. 32. K. F. P. Bennett, h. c. Lim, M. J. Braun, Sexual selection and introgression in avian hybrid 28, 739–755 (2015). zones: Spotlight on Manacus. Integr. Comp. Biol. 61, 1291–1309 (2021). 5. L. Rowe, h. d. Rundle, the alignment of natural and sexual selection. Annu. Rev. Ecol. Evol. 33. R. Maia, h. Gruson, J. A. endler, t. e. White, pavo 2: New tools for the spectral and spatial Syst. 52, 499–517 (2021). analysis of colour in r. Methods Ecol. Evol. 10, 1097–1107 (2019). Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 17 of 20 S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e 34. A. M. Rijke, Wettability and phylogenetic development of feather structure in water birds. 59. N. d. Sly, Golden- collared Manakin (Manacus vitellinus), version 2.0, Birds of the World J. Exp. Biol. 52, 469–479 (1970). (2023). https://birdsoftheworld.org/bow/species/gocman1/cur/introduction. 35. A. M. Rijke, W. A. Jesser, S. W. evans, h. Bouwman, Water repellency and feather structure 60. S. Andersson, M. Prager, “Quantifying colors” in Bird Coloration, G. hill, K. McGraw, eds. of the Blue Swallow Hirundo atrocaerulea. Ostrich 71, 143–145 (2000). (harvard Univ. Press, 2006), pp. 81. 36. t. L. Parchman, Z. Gompert, M. J. Braun, R. t. Brumfield, d. B. Mcdonald, J. A. c. Uy, 61. d. A. Gray, carotenoids and Sexual dichromatism in North American Passerine Birds. Am. G. Zhang, e. d. Jarvis, B. A. Schlinger, c. A. Buerkle, the genomic consequences of Nat. 148, 453–480 (1996). adaptive divergence and reproductive isolation between species of manakins. Mol. Ecol. 62. S. L. Olson, Specializations of some carotenoid- bearing feathers. Condor 72, 424–430 22, 3304–3317 (2013). (1970). 37. M. A. K. Widjaja-A dhi, M. Golczak, the molecular aspects of absorption and metabolism 63. A. L. Potticary, e. S. Morrison, A. v. Badyaev, turning induced plasticity into refined of carotenoids and retinoids in vertebrates. Biochim. Biophys. Acta Mol. Cell Biol. Lipids adaptations during range expansion. Nat. Commun. 11, 3254 (2020). 1865, 158571 (2020). 64. d. M. troy, A. h. Brush, Pigments and feather structure of the redpolls, Carduelis flammea 38. h. Bellil, F. Ghieh, e. hermel, B. Mandon- Pepin, F. vialard, human testis- expressed (teX) and C. hornemanni. Condor 85, 443–446 (1983). genes: A review focused on spermatogenesis and male fertility. Basic Clin. Androl. 31, 9 65. d. e. Mccoy, A. J. Shultz, c. vidoudez, e. van der heide, J. e. dall, S. A. trauger, d. haig, (2021). Microstructures amplify carotenoid plumage signals in tanagers. Sci. Rep. 11, 8582 39. B. thöny, G. Auerbach, N. Blau, tetrahydrobiopterin biosynthesis, regeneration and (2021). functions. Biochem J 347, 1–16 (2000). 66. M. A. Gazda, M. B. toomey, P. M. Araújo, R. J. Lopes, S. Afonso, c. A. Myers, K. Serres, 40. t. e. cruickshank, M. W. hahn, Reanalysis suggests that genomic islands of speciation P. d. Kiser, G. e. hill, J. c. corbo, M. carneiro, Genetic basis of de novo appearance of are due to reduced diversity, not reduced gene flow. Mol. Ecol. 23, 3133–3157 carotenoid ornamentation in bare parts of canaries. Mol. Biol. Evol. 37, 1317–1328 (2020). (2014). 67. J. B. Pease, R. J. driver, d. A. de la cerda, L. B. day, W. R. Lindsay, B. A. Schlinger, 41. R. Burri, interpreting differentiation landscapes in the light of long- term linked selection. e. R. Schuppe, c. N. Balakrishnan, M. J. Fuxjager, Layered evolution of gene expression in Evol. Lett. 1, 118–131 (2017). “superfast” muscles for courtship. Proc. Natl. Acad. Sci. U.S.A. 119, e2119671119 (2022). 42. M. Nei, Molecular Evolutionary Genetics (columbia Univ. Press, 1987). 68. W. P. Maddison, Gene trees in Species trees. Syst. Biol. 46, 523–536 (1997). 43. S. h. Martin, J. W. davey, c. d. Jiggins, evaluating the use of ABBA–BABA statistics to 69. R. N. Leite, R. t. Kimball, e. L. Braun, e. P. derryberry, P. A. hosner, G. e. derryberry, locate introgressed loci. Mol. Biol. Evol. 32, 244–257 (2015). M. Anciães, J. S. McKay, A. Aleixo, c. c. Ribas, R. t. Brumfield, J. cracraft, Phylogenomics of 44. J. L. Feder, X. Xie, J. Rull, S. velez, A. Forbes, B. Leung, h. dambroski, K. e. Filchak, M. Aluja, manakins (Aves: Pipridae) using alternative locus filtering strategies based on Mayr, dobzhansky, and Bush and the complexities of sympatric speciation in Rhagoletis. informativeness. Mol. Phylogenet. Evol. 155, 107013 (2021). Proc. Natl. Acad. Sci. U.S.A. 102, 6573–6580 (2005). 70. d. B. Mcdonald, t. L. Parchman, M. R. Bower, W. A. hubert, F. J. Rahel, An introduced and a 45. J. Strychalski, P. Brym, U. czarnik, A. Gugołek, A novel AAt-d eletion mutation in the native vertebrate hybridize to form a genetic bridge to a second native species. Proc. coding sequence of the BCO2 gene in yellow- fat rabbits. J. Appl. Genetics 56, 535–537 Natl. Acad. Sci. U.S.A. 105, 10837–10842 (2008). (2015). 71. P. J. Wilson, L. Y. Rutledge, t. J. Wheeldon, B. R. Patterson, B. N. White, Y- chromosome 46. d. P. L. toews, S. A. taylor, R. vallender, A. Brelsford, B. G. Butcher, P. W. Messer, i. J. Lovette, evidence supports widespread signatures of three- species Canis hybridization in eastern Plumage genes and little else distinguish the genomes of hybridizing warblers. Curr. Biol. North America. Ecol. Evol. 2, 2325–2332 (2012). 26, 2313–2318 (2016). 72. M. L. haines, J. Melville, J. Sumner, N. clemann, d. G. chapple, d. Stuart- Fox, Geographic 47. M. A. Gazda, P. M. Araújo, R. J. Lopes, M. B. toomey, P. Andrade, S. Afonso, c. Marques, variation in hybridization and ecological differentiation between three syntopic, L. Nunes, P. Pereira, S. trigo, G. e. hill, J. c. corbo, M. carneiro, A genetic mechanism for morphologically similar species of montane lizards. Mol. Ecol. 25, 2887–2903 (2016). sexual dichromatism in birds. Science 368, 1270–1274 (2020). 73. Q. K. Langdon, d. Peris, e. P. Baker, d. A. Opulente, h.- v. Nguyen, U. Bond, P. Gonçalves, 48. J. eriksson, G. Larson, U. Gunnarsson, B. Bed’hom, M. tixier- Boichard, L. Strömstedt, J. P. Sampaio, d. Libkind, c. t. hittinger, Fermentation innovation through complex d. Wright, A. Jungerius, A. vereijken, e. Randi, P. Jensen, L. Andersson, identification of the hybridization of wild and domesticated yeasts. Nat. Ecol. Evol. 3, 1576–1586 (2019). yellow skin gene reveals a hybrid origin of the domestic chicken. PLOS Genet. 4, 74. P. R. Grant, B. R. Grant, triad hybridization via a conduit species. Proc. Natl. Acad. Sci. U.S.A. e1000010 (2008). 117, 7888–7896 (2020). 49. M. d. Baiz, A. W. Wood, A. Brelsford, i. J. Lovette, d. P. L. toews, Pigmentation genes show 75. L. Natola, S. S. Seneviratne, d. irwin, Population genomics of an emergent tri- species evidence of repeated divergence and multiple bouts of introgression in Setophaga hybrid zone. Mol. Ecol. 31, 5356–5367 (2022). warblers. Curr. Biol. 31, 643–649.e3 (2021). 76. M. S. Ferreira, t. J. thurman, M. R. Jones, L. Farelo, A. v. Kumar, S. M. e. Mortimer, 50. P. Andrade, c. Pinho, G. Pérez i de Lanuza, S. Afonso, J. Brejcha, c.- J. Rubin, O. Wallerman, J. R. demboski, L. S. Mills, P. c. Alves, J. Melo- Ferreira, J. M. Good, the evolution of P. Pereira, S. J. Sabatino, A. Bellati, d. Pellitteri- Rosa, Z. Bosakova, i. Bunikis, M. A. carretero, white-t ailed jackrabbit camouflage in response to past and future seasonal climates. N. Feiner, P. Marsik, F. Paupério, d. Salvi, L. Soler, G. M. While, t. Uller, e. Font, L. Andersson, Science 379, 1238–1242 (2023). M. carneiro, Regulatory changes in pterin and carotenoid genes underlie balanced color 77. J. i. Meier, M. d. McGee, d. A. Marques, S. Mwaiko, M. Kishe, S. Wandera, d. Neumann, polymorphisms in the wall lizard. Proc. Natl. Acad. Sci. U.S.A. 116, 5633–5642 (2019). h. Mrosso, L. J. chapman, c. A. chapman, L. Kaufman, A. taabu-M unyaho, c. e. Wagner, 51. S. J. Lehnert, K. A. christensen, W. e. vandersteen, d. Sakhrani, t. e. Pitcher, J. W. heath, R. Bruggmann, L. excoffier, O. Seehausen, cycles of fusion and fission enabled rapid B. F. Koop, d. d. heath, R. h. devlin, carotenoid pigmentation in salmon: variation in parallel adaptive radiations in African cichlids. Science 381, eade2833 (2023). expression at BCO2- l locus controls a key fitness trait affecting red coloration. Proc. R. Soc. 78. N. B. edelman, P. B. Frandsen, M. Miyagi, B. clavijo, J. davey, R. B. dikow, B Biol. Sci. 286, 20191588 (2019). G. García-A ccinelli, S. M. van Belleghem, N. Patterson, d. e. Neafsey, R. challis, S. Kumar, 52. e. d. enbody, c. G. Sprehn, A. Abzhanov, h. Bi, M. P. dobreva, O. G. Osborne, c.- J. Rubin, G. R. P. Moreira, c. Salazar, M. chouteau, B. A. counterman, R. Papa, M. Blaxter, R. d. Reed, P. R. Grant, B. R. Grant, L. Andersson, A multispecies BCO2 beak color polymorphism in K. K. dasmahapatra, M. Kronforst, M. Joron, c. d. Jiggins, W. O. McMillan, F. di Palma, the darwin’s finch radiation. Curr. Biol. 31, 5597–5604.e7 (2021). A. J. Blumberg, J. Wakeley, d. Jaffe, J. Mallet, Genomic architecture and introgression 53. c. R. Linnen, e. P. Kingsley, J. d. Jensen, h. e. hoekstra, On the origin and spread of an shape a butterfly radiation. Science 366, 594–599 (2019). adaptive allele in deer mice. Science 325, 1095–1098 (2009). 79. e. J. Beckman, P. M. Benham, Z. A. cheviron, c. c. Witt, detecting introgression despite 54. M. Abolins- Abols, e. Kornobis, P. Ribeca, K. Wakamatsu, M. P. Peterson, e. d. Ketterson, phylogenetic uncertainty: the case of the South American siskins. Mol. Ecol. 27, B. Milá, differential gene regulation underlies variation in melanic plumage coloration in 4350–4367 (2018). the dark-e yed junco (Junco hyemalis). Mol. Ecol. 27, 4501–4515 (2018). 80. A.-L . ducrest, L. Keller, A. Roulin, Pleiotropy in the melanocortin system, coloration and 55. Y.- R. Jin, J. K. Yoon, the R- spondin family of proteins: emerging regulators of WNt behavioural syndromes. Trends Ecol. Evol. 23, 502–510 (2008). signaling. Int. J. Biochem. Cell Biol. 44, 2278–2287 (2012). 81. R. J. Weaver, e. S. A. Santos, A. M. tucker, A. e. Wilson, G. e. hill, carotenoid metabolism 56. A. i. vickrey, R. Bruders, Z. Kronenberg, e. Mackey, R. J. Bohlender, e. t. Maclary, R. Maynez, strengthens the link between feather coloration and individual quality. Nat. Commun. 9, e. J. Osborne, K. P. Johnson, c. d. huff, M. Yandell, M. d. Shapiro, introgression of 73 (2018). regulatory alleles and a missense coding mutation drive plumage pattern diversity in the 82. R. t. Brumfield, M. J. Braun, Phylogenetic relationships in bearded manakins (Pipridae: rock pigeon. elife 7, e34803 (2018). Manacus) indicate that male plumage color is a misleading taxonomic marker. Condor 57. e. cadieu, M. W. Neff, P. Quignon, K. Walsh, K. chase, h. G. Parker, B. M. vonholdt, A. Rhue, 103, 248–258 (2001). A. Boyko, A. Byers, A. Wong, d. S. Mosher, A. G. elkahloun, t. c. Spady, c. André, K. G. Lark, 83. J. A. c. Uy, A. c. Stein, variable visual habitats may influence the spread of colourful M. cargill, c. d. Bustamante, R. K. Wayne, e. A. Ostrander, coat variation in the domestic plumage across an avian hybrid zone. J. Evol. Biol. 20, 1847–1858 (2007). dog is governed by variants in three genes. Science 326, 150–153 (2009). 84. e. h. duval, c. L. Fitzpatrick, e. A. hobson, M. R. Servedio, inferred Attractiveness: A 58. N. chaves, Orange- collared Manakin (Manacus aurantiacus), version 1.0, Birds of the generalized mechanism for sexual selection that can maintain variation in traits and World (2020). https://birdsoftheworld.org/bow/species/orcman1/cur/introduction. preferences over time. PLOS Biol. 21, e3002269 (2023). Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 18 of 20 S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e 85. W. A. Boyle, e. h. Shogren, J. d. Brawn, hygric Niches for tropical endotherms. Trends Ecol. 117. d. Bates, M. Mächler, B. Bolker, S. Walker, Fitting linear mixed-e ffects models using lme4. Evol. 35, 938–952 (2020). J. Stat. Softw. 67, 1–48 (2015). 86. R. J. Kennedy, direct effects of rain on birds: A review. British Birds 63, 401–414 (1970). 118. t. L. Pedersen, patchwork: the composer of Plots, version 1.1.2 (2022); https:// 87. S. e. Fick, R. J. hijmans, Worldclim 2: New 1- km spatial resolution climate surfaces for cran.r-p roject.org/web/packages/patchwork/index.html. global land areas. Int. J. Climatol. 37, 4302–4315 (2017). 119. J. haffer, Speciation in colombian forest birds west of the Andes. Am. Mus. Novit. 2294, 88. F. M. S. Muzio, M. A. Rubega, What do we really know about the water repellency of 58 (1967). feathers? J. Avian Biol. 10.1111/jav.03259 , (2024). 120. M. G. harvey, G. A. Bravo, S. claramunt, A. M. cuervo, G. e. derryberry, J. Battilana, 89. R core team, R: A language and environment for statistical computing, R Foundation for G. F. Seeholzer, J. S. McKay, B. c. O’Meara, B. c. Faircloth, S. v. edwards, J. Pérez- emán, Statistical computing (2024); https://R-p roject.org/. R. G. Moyle, F. h. Sheldon, A. Aleixo, B. t. Smith, R. t. chesser, L. F. Silveira, J. cracraft, 90. R. Lenth, emmeans: estimated Marginal Means, aka Least- Squares Means, (2020); https:// R. t. Brumfield, e. P. derryberry, the evolution of a tropical biodiversity hotspot. Science cran.r- project.org/package=emmeans. 370, 1343–1348 (2020). 91. A. Kassambara, ggpubr: “ggplot2” based publication ready plots, (2019); https:// 121. K. J. McGraw, J. hudon, G. e. hill, R. S. Parker, A simple and inexpensive chemical test for cran.r-p roject.org/package=ggpubr. behavioral ecologists to determine the presence of carotenoid pigments in animal 92. c. O. Wilke, cowplot: Streamlined Plot theme and Plot Annotations for “ggplot2,” cRAN tissues. Behav. Ecol. Sociobiol. 57, 391–397 (2005). Repos (2016); https://cran.r- project.org/package=cowplot. 122. K. J. McGraw, e. Adkins- Regan, R. S. Parker, Anhydrolutein in the zebra finch: A new, 93. B. Langmead, S. L. Salzberg, Fast gapped- read alignment with Bowtie 2. Nat. Methods. 9, metabolically derived carotenoid in birds. Comp. Biochem. Physiol. B Biochem. Mol. Biol. 357–359 (2012). 132, 811–818 (2002). 94. h. Li, B. handsaker, A. Wysoker, t. Fennell, J. Ruan, N. homer, G. Marth, G. Abecasis, 123. S. W. Wright, S. W. Jeffrey, R. F. c. Mantoura, c. A. Llewelly, t. Bjornland, d. Repeta, R. durbin; 1000 Genome Project data Processing Subgroup, the sequence alignment/ N. Welschmeyer, improved hPLc method for the analysis of chlorophylls and carotenoids map format and SAMtools. Bioinformatics 25, 2078–2079 (2009). from marine phytoplankton. Mar. Ecol. Prog. Ser. 77, 183–196 (1991). 95. A. v. Zimin, G. Marçais, d. Puiu, M. Roberts, S. L. Salzberg, J. A. Yorke, the MaSuRcA 124. N. h. Barton, Genetic hitchhiking. Philos. Trans. R. Soc. Lond. B Biol. Sci. 355, 1553–1562 genome assembler. Bioinformatics 29, 2669–2677 (2013). (2000). 96. G. A. van der Auwera, B. d. O’connor, Genomics in the Cloud: Using Docker, GATK, and WDL 125. A. S. Y. Lee, P. J. Kranzusch, J. h. d. cate, eiF3 targets cell- proliferation messenger RNAs for in Terra (O’Reilly Media, 2020). translational activation or repression. Nature 522, 111–114 (2015). 97. M. Alonge, S. Soyk, S. Ramakrishnan, X. Wang, S. Goodwin, F. J. Sedlazeck, Z. B. Lippman, 126. S. Wang, S. Rohwer, d. R. de Zwaan, d. P. L. toews, i. J. Lovette, J. Mackenzie, d. irwin, M. c. Schatz, RaGOO: Fast and accurate reference- guided scaffolding of draft genomes. Selection on a small genomic region underpins differentiation in multiple color traits Genome Biol. 20, 224 (2019). between two warbler species. Evol. Lett. 4, 502–515 (2020). 98. d. Bryant, v. Moulton, Neighbor- net: An agglomerative method for the construction of 127. F. Liu, M. visser, d. L. duffy, P. G. hysi, L. c. Jacobs, O. Lao, K. Zhong, S. Walsh, L. chaitanya, phylogenetic networks. Mol. Biol. Evol. 21, 255–265 (2004). A. Wollstein, G. Zhu, G. W. Montgomery, A. K. henders, M. Mangino, d. Glass, v. Bataille, 99. d. h. huson, d. Bryant, Application of phylogenetic networks in evolutionary studies. R. A. Sturm, F. Rivadeneira, A. hofman, W. F. J. van iJcken, A. G. Uitterlinden, Mol. Biol. Evol. 23, 254–267 (2006). R.- J. t. S. Palstra, t. d. Spector, N. G. Martin, t. e. c. Nijsten, M. Kayser, Genetics of skin color 100. h. e. Lischer, L. excoffier, PGdSpider: An automated data conversion tool for connecting variation in europeans: Genome- wide association studies with functional follow- up. population genetics and genomics programs. Bioinformatics 28, 298–299 (2012). Hum Genet 134, 823–835 (2015). 101. t. S. Korneliussen, A. Albrechtsen, R. Nielsen, ANGSd: Analysis of next generation 128. c.- h. chang, t.- X. Jiang, c.- M. Lin, L. W. Burrus, c.- M. chuong, R. Widelitz, distinct Wnt sequencing data. BMC Bioinformatics 15, 356 (2014). members regulate the hierarchical morphogenesis of skin regions (spinal tract) and 102. M. Fumagalli, F. G. vieira, t. S. Korneliussen, t. Linderoth, e. huerta-S ánchez, individual feathers. Mech. Dev. 121, 157–171 (2004). A. Albrechtsen, R. Nielsen, Quantifying population genetic differentiation from 129. R. B. Widelitz, G.- W. Lin, Y.- c. Lai, J. A. Mayer, P.- c. tang, h.- c. cheng, t.- X. Jiang, c.- F. chen, next- generation sequencing data. Genetics 195, 979–992 (2013). c.- M. chuong, Morpho-r egulation in diverse chicken feather formation: integrating 103. A. J. coffman, P. h. hsieh, S. Gravel, R. N. Gutenkunst, computationally efficient branching modules and sex hormone-d ependent morpho- regulatory modules. Dev. composite likelihood statistics for demographic inference. Mol. Biol. Evol. 33, 591–593 Growth Differ. 61, 124–138 (2019). (2016). 130. K. F. Stryjewski, M. d. Sorenson, Mosaic genome evolution in a recent and rapid avian 104. t. R. Booker, S. Yeaman, M. c. Whitlock, variation in recombination rate affects detection radiation. Nat. Ecol. Evol. 1, 1912–1922 (2017). of outliers in genome scans under neutrality. Mol. Ecol. 29, 4274–4279 (2020). 131. R. Schiffer, M. Neis, d. höller, F. Rodríguez, A. Geier, c. Gartung, F. Lammert, A. dreuw, 105. A. Zeileis, G. Grothendieck, J. A. Ryan, J. M. Ulrich, F. Andrews, zoo: S3 infrastructure for G. Zwadlo-K larwasser, h. Merk, F. Jugert, J. M. Baron, Active influx transport is mediated Regular and irregular time Series (Z’s Ordered Observations), version 1.8-1 2, (2023); by members of the organic anion transporting polypeptide family in human epidermal https://cran.r- project.org/web/packages/zoo/index.html. keratinocytes. J. Invest. Dermatol. 120, 285–291 (2003). 106. h. Mi, A. Muruganujan, d. ebert, X. huang, P. d. thomas, PANtheR version 14: More 132. K. Sebastian, S. detro-d assen, N. Rinis, d. Fahrenkamp, G. Müller- Newen, h. F. Merk, genomes, a new PANtheR GO- slim and improvements in enrichment analysis tools. G. Schmalzing, G. Zwadlo- Klarwasser, J. M. Baron, characterization of SLcO5A1/ Nucleic Acids Res. 47, d419–d426 (2019). OAtP5A1, a solute carrier transport protein with non- classical function. PLOS ONE 8, 107. h. Mi, P. thomas, Protein networks and pathway analysis. Methods Mol. Biol. 563, 123–140 e83257 (2013). (2009). 1 33. J. Graf, R. hodgson, A. van daal, Single nucleotide polymorphisms in the MAtP gene are 108. A. dobin, c. A. davis, F. Schlesinger, J. drenkow, c. Zaleski, S. Jha, P. Batut, M. chaisson, associated with normal human pigmentation variation. Hum. Mutat. 25, 278–284 t. R. Gingeras, StAR: Ultrafast universal RNA- seq aligner. Bioinformatics 29, 15–21 (2013). (2005). 109. Y. Liao, G. K. Smyth, W. Shi, featurecounts: An efficient general purpose program for 134. L. campagna, M. Repenning, L. F. Silveira, c. S. Fontana, P. L. tubaro, i. J. Lovette, Repeated assigning sequence reads to genomic features. Bioinformatics 30, 923–930 (2014). divergent selection on pigmentation genes in a rapid finch radiation. Sci. Adv. 3, 110. h. Li, Aligning sequence reads, clone sequences and assembly contigs with BWA-M eM. e1602404 (2017). arXiv:1303.3997 [q- bio.GN] (2013). 135. c. Kiefer, S. hessel, J. M. Lampert, K. vogt, M. O. Lederer, d. e. Breithaupt, J. von Lintig, 111. P. danecek, A. Auton, G. Abecasis, c. A. Albers, e. Banks, M. A. dePristo, R. e. handsaker, identification and characterization of a mammalian enzyme catalyzing the asymmetric G. Lunter, G. t. Marth, S. t. Sherry, G. Mcvean, R. durbin; 1000 Genomes Project Analysis oxidative cleavage of provitamin A. J. Biol. Chem. 276, 14110–14116 (2001). Group, the variant call format and vcFtools. Bioinformatics 27, 2156–2158 (2011). 136. c. dela Seña, J. Sun, S. Narayanasamy, K. M. Riedl, Y. Yuan, R. W. curley Jr., S. J. Schwartz, 112. K. Katoh, d. M. Standley, MAFFt multiple sequence alignment software version 7: e. h. harrison, Substrate specificity of purified recombinant chicken β- carotene improvements in performance and usability. Mol. Biol. Evol. 30, 772–780 (2013). 9′,10′- oxygenase (BcO2). J. Biol. Chem. 291, 14609–14619 (2016). 113. B. Q. Minh, h. A. Schmidt, O. chernomor, d. Schrempf, M. d. Woodhams, A. von haeseler, 137. J. Amengual, G. P. Lobo, M. Golczak, h. N. M. Li, t. Klimova, c. L. hoppel, A. Wyss, R. Lanfear, iQ- tRee 2: New models and efficient methods for phylogenetic inference in K. Palczewski, J. von Lintig, A mitochondrial enzyme degrades carotenoids and protects the genomic era. Mol. Biol. Evol. 37, 1530–1534 (2020). against oxidative stress. FASEB J 25, 948–959 (2011). 114. P. danecek, J. K. Bonfield, J. Liddle, J. Marshall, v. Ohan, M. O. Pollard, A. Whitwham, 138. G. P. Lobo, A. isken, S. hoff, d. Babino, J. von Lintig, BcdO2 acts as a carotenoid scavenger t. Keane, S. A. Mccarthy, R. M. davies, h. Li, twelve years of SAMtools and BcFtools. and gatekeeper for the mitochondrial apoptotic pathway. Development 139, 2966–2977 GigaScience 10, giab008 (2021). (2012). 115. M. N. Price, P. S. dehal, A. P. Arkin, Fasttree 2—Approximately maximum-l ikelihood trees 139. S. d. Berry, S. R. davis, e. M. Beattie, N. L. thomas, A. K. Burrett, h. e. Ward, A. M. Stanfield, for large alignments. PLOS ONE 5, e9490 (2010). M. Biswas, A. e. Ankersmit- Udy, P. e. Oxley, J. L. Barnett, J. F. Pearson, Y. van der does, 116. A. M. Rijke, W. A. Jesser, the water penetration and repellency of feathers revisited. The A. h. K. MacGibbon, R. J. Spelman, K. Lehnert, R. G. Snell, Mutation in Bovine β-c arotene Condor 113, 245–254 (2011). oxygenase 2 affects milk color. Genetics 182, 923–926 (2009). Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 19 of 20 S c i e N c e A d v A N c e S | R e S e A R c h A R t i c L e 140. d. i. våge, i. A. Boman, A nonsense mutation in the beta- carotene oxygenase 2 (BcO2) Museum of Natural history Research Grant (M.J.B. and K.F.P.B.); US department of Agriculture, gene is tightly associated with accumulation of carotenoids in adipose tissue in sheep National institute of Food and Agriculture, hatch Project 1026333 (iLLU- 875-9 84) to K.M.L. and (Ovis aries). BMC Genetics 11, 10 (2010). J.d.B.; National Science Foundation grant iOS 1122180 (L.B.d.); National Science Foundation 141. c.- M. chuong, the making of a feather: homeoproteins, retinoids and adhesion grant iOS 1947472 (M.J.F.); National Science Foundation grant iOS 203774 (G.e.h.); National molecules. BioEssays 15, 513–521 (1993). Science Foundation grant dBi 1457541 (M.J.B. and c.N.B.); University of Maryland devra 142. Y. Mishima, Molecular and biological control of melanogenesis through tyrosinase genes Kleiman Fellowship (K.F.P.B.); University of Maryland BiSi dissertation Fellowship (K.F.P.B.); and and intrinsic and extrinsic regulatory factors. Pigment Cell Res. 7, 376–387 (1994). Smithsonian Biodiversity Genomics Postdoctoral Fellowship (h.c.L.). Author contributions: 143. N. i. Mundy, A window on the genetics of evolution: Mc1R and plumage colouration in conceptualization: M.J.B., h.c.L., K.F.P.B., K.M.L., J.d.B., and G.e.h. data collection: h.c.L., K.F.P.B., birds. Proc. Biol. Sci. 272, 1633–1640 (2005). N.M.J., M.J.P., K.M.L., L.B.d., W.R.L., and t.J.P. data analysis: h.c.L., K.F.P.B., K.M.L., M.J.P., N.M.J., 144. e. M. Wolf horrell, M. c. Boulanger, J. A. d’Orazio, Melanocortin 1 receptor: Structure, J.B.P., S.e.K., P.e.B., and M.J.F. visualization: K.F.P.B., h.c.L., M.J.P., N.M.J., and K.M.L. Supervision: function, and regulation. Front. Genet. 7, 95 (2016). M.J.B., J.d.B., G.e.h., and c.N.B. Writing—original draft: h.c.L., K.F.P.B., M.J.B., N.M.J., M.J.P., G.e.h., and K.M.L. Writing—review and editing: K.F.P.B., h.c.L., M.J.B., M.J.P., N.M.J., G.e.h., P.e.B., J.B.P., Acknowledgments: We gratefully acknowledge J. Wilkinson, S. Mount, c. Machado, S.e.K., L.B.d., t.J.P., K.M.L., J.d.B., M.J.F., and c.N.B. Competing interests: the authors declare K. carleton, h. Fisher, M. Baldwin, J. von Lintig, B. Schlinger, L. Fusani, O. McMillan, that they have no competing interests. Data and materials availability: Specimens and e. Bermingham, e. Braun, F. Muzio, and c. dove for providing valuable discussions and genetic resource materials used to generate data for this study are housed at the US NMNh, comments on study design and drafts of this paper. Sample collection: S. Olson, J. Blake, the Louisiana State University Museum of Natural Science, the Academy of Natural Sciences of R. Brumfield, t. Glenn, and d. Mcdonald. Field assistance: A. Pineda Sr. and A. Pineda Jr. drexel University, and the University of Mississippi. collection access is governed by museum Specimen access: We thank c. Milensky and c. huddleston for facilitating access to specimens policies, which may be found at https://naturalhistory.si.edu/research/vertebrate-z oology/ and samples housed in the Smithsonian National Museum of Natural history Bird collection birds, https://apa.lsu.edu/mns/collections/ornithology.php, and https://ansp.org/research/. and Biorepository, respectively. Permits: We are grateful to the following institutions for data and code necessary to reproduce the analyses are available from NcBi (BioProjects: providing collecting and access permits: Autoridad Nacional del Ambiente and Autoridad del PRJNA1129471, PRJNA893627, and PRJNA321179) and on dryad. dOi: 10.5061/dryad. canal de Panamá, the Guyana environmental Protection Agency, and Guyana Ministry of jsxksn0hn. All other data needed to evaluate the conclusions in the paper are present in the Amerindian Affairs. the following people assisted in permit processing: M. Leone, G. Maggiori, paper and/or the Supplementary Materials. L. camacho, and J. Mate. computational resources: computations of genomic data in this paper were conducted on the Smithsonian high- Performance cluster (Si/hPc), Smithsonian Submitted 2 January 2024 institution (https://doi.org/10.25572/SihPc). We thank R. dikow for providing technical Accepted 17 October 2024 assistance. Funding: this work was supported by Smithsonian institution competitive Pell Published 20 November 2024 Grant (M.J.B. and h.c.L.); Smithsonian institution James Bond Fund (M.J.B. and K.F.P.B.); National 10.1126/sciadv.adn8339 Lim et al., Sci. Adv. 10, eadn8339 (2024) 20 November 2024 20 of 20