Testing possible demographic scenarios
We used model-based demographic inference to evaluate support for alternative scenarios describing the evolutionary dynamics of contact zones between subspecies. We built demographic models using the program MOMI2 (Kamm et al., 2019), which uses the joint Site Frequency Spectra (SFS) to evaluate different models that consider population size, divergence time, and gene flow for each clade. We ran the demographic inferences based on a mutation rate of 2.5 × 10-9sites per generation (Nadachowska-Brzyska et al., 2015) and a generation time of 2.2 yr, which is a common generation time for Neotropical songbirds (Luzuriaga-Aveiga et al., 2021). We executed our models in two steps (as in Corbett et al., 2020). First, we built a model that estimated effective population size for the three lineages, based on the topology and time of divergence of the mitochondrial (ND2) phylogenetic reconstruction (Ocampo et al., 2022a), in which the two pied subspecies (i.e., S. c. hicksii, S. c. hoffmanni ) are sister taxa to the black S. c. corvina (Figure 2A). Second, we built demographic models that added bidirectional bouts of gene flow and estimated the strength of the bouts.
Demographic models were built to test for two specific questions: 1) Is the apparent absence of corvina-hoffmanni hybrids at their contact zone consistent with a scenario of no contemporary gene flow between these subspecies? We tested whether gene flow (bouts of 10-25% of the population) occurred between S. c. corvina and S. c. hoffmanni 1000 years ago. 2) Did S. c. corvina diverge in isolation until recent secondary contact with S. c. hicksii or have these subspecies been in constant contact? We tested for the likelihood of constant gene flow by modeling temporal pulses of 5-10% of the population, occurring at 200, 400, and 600 Kyr ago to distinguish between a scenario of historical isolation or constant gene flow. Given that we found evidence for later-generation hybrids at bothcorvina-hicksii and hicksii-hoffmanni contact zones (see Results), we modeled bidirectional bouts of 10-25% of the population 1000 years ago and kept these gene flow events constant across all our models. We tested our two hypotheses in all possible combinations, for a total of four demographic models (Figure 2B-E), estimating the strength of the gene flow pulses. All models were optimized using the “L-BFGS-B” algorithm and a maximum of 100 iterations. With both steps, we carried out 100 runs per model, and selected the most likely set of parameters of each model based on log-likelihood values. We then compared between different demographic models using the AIC to select the best model. We did not fit more complex models (e.g., changes in effective population size or more complex scenarios of gene flow), as these simple models allowed us to test our two main questions and avoided overparameterization (Bocalini et al., 2021).
RESULTS