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 i – X min) /
(X max – X 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.