Characterization of hybrid individuals and contact zones
We characterized first-generation and later-generation hybrids to distinguish between current hybridization versus ancient gene flow. Our pairwise Weir and Cockerham’s FST per SNP between pure parental subspecies allowed us to identify informative markers for subspecies ancestry (FST > 0.60). We then used these sets of informative markers to estimate the hybrid index and level of heterozygosity from parental subspecies and their corresponding hybrids on the R package INTROGRESS (Gompert & Buerkle, 2009). This method is commonly used to characterize hybrids, as it does not require fixed markers and broadly distinguishes ongoing hybridization (F1 hybrids with high heterozygosity) from historic hybridization and backcrosses (later-generation hybrids with low heterozygosity, e.g., Scordato et al., 2017).
We also characterized the contact zones between hybridizing subspecies (corvina-hicksii , hicksii-hoffmanni ; see Results) fitting geographic clines (Szymura & Barton, 1986). Geographic cline analyses were performed by fitting genetic and phenotypic clines using the R package HZAR (Derryberry et al., 2014), which allowed us to test how different traits covary along the geographic transition between two subspecies. Because contact zones can roughly align with the Caribbean (corvina-hicksii ; Figure 1A, no. 1) or Pacific slopes (hoffmanni-hicksii ; Figure 1A, no. 2), we set a line along the continental divide following some of the major highland peaks along the Costa Rica–Panama mountain ranges (Rincón de la Vieja, Poás, Chirripó, Barú), throughout the middle of Panama (Santiago, La Yeguada, Gamboa, El Llano, Higueronal, Metetí, Yaviza) to the Panamanian border with Colombia. Then, we assigned samples collected from nearby locations to populations, and extrapolated the average coordinates per population to the closest point on the line, using the “dist2line” function from the GEOSPHERE R package (Hijmans et al., 2017). We estimated geographic distances between each inferred population along this line, starting at Rincón de la Vieja volcano (km 0) using the “trackDistance” function from the R package TRIP (Sumner et al., 2009). This approach allowed us to convert the irregular topography of this region into a single linear axis against which both contact zones can be mapped (Figure S1).
For the phenotypic traits, we carried out cline analyses using only traits that differed between subspecies. First, we performed cline analyses on the male plumage traits that differ by changing from black to white between subspecies (i.e., corvina-hicksii : throat, belly, and rump, hoffmanni-hicksii : throat). For the morphometric traits, we ran a single analysis of variance (ANOVA) per trait and per contact zone, but included sex as a covariable, and tested for its interaction with the subspecies factor. We performed the ANOVA between the most-distant locations (localities Arenal [ARE] and Sarapiquí [SAR] for S. c. corvina ; Garabito [GAR] and Quepos [QUE] for S. c. hoffmanni ; and Cañazas [CAN] and Darién [DAR] for S. c. hicksii ; Table S1). Because we found no significant effect of sex for any of the morphometric traits, we estimated clines for the morphological traits including both males and females. For the phenotypic traits that differed between the most-distant populations, we converted each phenotypic trait into scores that ranged from 0 to 1. Scores values were calculated by applying the equation (X iX min) / (X maxX min), in whichX i corresponds to each trait’s value per individual and X min andX max correspond to the minimum and maximum trait values in each cline comparison, respectively.
For the genetic clines, we modeled 15 possible clines on the average hybrid indices by fitting three different settings for scaling at the ends of the clines (fixed, free, and no-scaling), and five possible fits of the tails of the clines (none, right, left, both, and mirroring), in all possible combinations (per Derryberry et al., 2014). We replicated each model three times with different random seeds, ran three chains per model fit to assess convergence, and selected for the best model based on the Akaike Information Criterion (AIC). For the phenotypic clines, morphometric traits and brightness per body patch, we tested ten different models by applying five different fits to the cline’s tails (none, right, left, both, and mirroring) and two settings for scaling (fixed and free), by fixing the scaling to the maximum and minimum scores. Likewise, for the hybrid index, we ran three chains and three replicates per model, and selected the best model based on the lowest AIC. We compared coincidence and concordance in clines (i.e., center and width) between different traits and hybrid index using the 95% confidence intervals (CI) for each cline parameter.