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