Assortative mating and epistatic mating‐trait architecture induce complex movement of the crow hybrid zone

Hybrid zones provide a window into the evolutionary processes governing species divergence. Yet, the contribution of mate choice to the temporal and spatial stability of hybrid zones remains poorly explored. Here, we investigate the effects of assortative mating on hybrid‐zone dynamics by means of a mathematical model parameterized with phenotype and genotype data from the hybrid zone between all‐black carrion and gray‐coated hooded crows. In the best‐fit model, narrow clines of the two mating‐trait loci were maintained by a moderate degree of assortative mating inducing pre‐ and postzygotic isolation via positive frequency‐dependent selection. Epistasis between the two loci induced hybrid‐zone movement in favor of alleles conveying dark plumage followed by a shift in the opposite direction favoring gray‐coated phenotypes ∼1 200 generations after secondary contact. Unlinked neutral loci diffused near‐unimpeded across the zone. These results were generally robust to the choice of matching rule (self‐referencing or parental imprinting) and effects of genetic drift. Overall, this study illustrates under which conditions assortative mating can maintain steep clines in mating‐trait loci without generalizing to genome‐wide reproductive isolation. It further emphasizes the importance of the genetic mating‐trait architecture for spatio‐temporal hybrid‐zone dynamics.

potentially overlap in their ranges, provided there is little competition. Incomplete isolation will result in hybrid zones (Bigelow 1965;Hewitt 1988;Harrison 1990;Rieseberg et al. 1999). These zones act as semipermeable barriers retaining loci under divergent selection, while unlinked neutral variants can still introgress (Key 1968;Payseur 2010). When the homogenizing effect of dispersal into the zone is offset by selection against hybrids, the resulting hybrid zone can reach relative stability in the form of sigmoid clines ("tension-zone model"; Barton and Hewitt 1985;McEntee et al. 2018). As a result, unique clines are expected for different traits and the underlying genes, unless coupling elicits a genome-wide response (Barton 1983;Feder et al. 2014;Flaxman et al. 2014). Cline width and position will be influenced by the strength of selection, dispersal distances, and the genetic architecture of the trait subject to selection (May et al. 1975;Endler 1977;Mallet 1986). Importantly, hybrid zones need not be static. Asymmetries in relative fitness between parental populations or a dominant genetic architecture of the loci contributing to divergence may affect cline centers and induce hybrid-zone movement until they are caught in a trough of low population density, often aligning with physical or ecological barriers (Mallet 1986;Secondi et al. 2006;Brodin et al. 2013).
Mathematical models and computer simulations within the framework of cline analyses have led to insights on the effects of selection strength (Barton and Hewitt 1985), the underlying nature of selection (exogenous vs. endogenous) (Kruuk et al. 1999), the genetic architecture of the contributing loci (Mallet 1986;Mallet and Barton 1989;Bürger 2017) and their complex interactions (Barton 1983;Kruuk et al. 1999; Barton and Shpak 2000b). Most models, however, rely on the assumption of random mating, and few have explored the role of isolation arising from mate choice alone (but see Irwin 2020). This gap deserves attention, given the importance attributed to premating isolation, in particular during the early stages of speciation (Mayr 1963;Coyne and Orr 1997;Stelkens et al. 2010).
Yet, one of the major challenges for speciation to occur by sexual section is recombination breaking down the association between mating cues and corresponding preferences (Verzijden et al. 2012;Kopp et al. 2018). This obstacle can be overcome in two ways: physical linkage of the loci underlying genetically determined trait and preference (Merrill et al. 2019;Xu and Shaw 2019), or a matching rule aligning the phenotype of a genetically encoded mating trait with a corresponding preference. Among the mechanisms that can lead to such a matching rule are imprinting on parental phenotypes and self-referencing (Hauber et al. 2000;Sherman 2001, 2003;Kopp et al. 2018). Evidence for imprinting of mating traits has been reported for many taxonomic groups, including birds, mammals, fish, amphibians, and even insects (see Verzijden et al. 2012;Kopp et al. 2018, for references). Hauber and Sherman (2001) review evidence for self-referencing in several avian and mammalian species, for example, found in cowbirds (Hauber et al. 2000, see also Hare et al. 2003, andSherman 2003).
With assortative mating according to matching rules, frequencies of preference and trait are per se correlated, which can result in positive frequency-dependent selection. With some degree of incipient divergence (secondary sympatry, geographic isolation), the correlation between trait and preference has the potential to accelerate speciation (Servedio 2016). Indeed, the preference for sexually selected traits, for example, as an effect of self-referencing or imprinting on parental phenotypes, has been shown to promote phenotypic divergence and maintain polymorphism in the face of gene flow (Brodin and Haas 2006;Verzijden et al. 2012;Yang et al. 2019). Empirical work in hybrid zones regularly identifies sexually selected mating traits and preferences as an important component of reproductive isolation (Seehausen et al. 2008;Schumer et al. 2017;Hench et al. 2019;Yang et al. 2019). However, simulation studies suggest that assortative mating can maintain diverged mating traits, but impede genome-wide gene flow of unlinked genetic variation just slightly (Brodin and Haas 2009;Irwin 2020).
Here, we explore the role of assortative mating in the wellstudied hybrid zone between all-black carrion and gray-coated hooded crows (Corvus (corone) corone and C. (c.) cornix) that presumably arose by secondary contact in the early Holocene approximately 12,000 years (or 2000 crow generations) ago (Meise 1928;Mayr 1942;Parkin et al. 2003;Vijay et al. 2016;Kopp et al. 2018). In this system, there is only limited evidence for natural selection against hybrids (Saino 1990;Saino and Bolzern 1992;Saino and Villa 1992), and an ecological contribution to the hybrid-zone dynamics exceeding local effects can almost certainly be excluded (Rolando and Laiolo 1994;Saino et al. 1998;Randler 2007b;Haas et al. 2010). Instead, there is support for a role of plumage-based assortative mating (Meise 1928;Randler 2007a) and social marginalization of minority phenotypes (Saino and Scatizzi 1991;Londei 2013) maintaining the hybrid zone against gene flow (Poelstra et al. 2014).
Plumage pigmentation patterns are assumed to be used as a signal for mate choice, and constitute the only diagnostic difference between taxa (Blotzheim et al. 1993). Sexes do not differ in phenotype suggesting the possibility of mutual, assortative mate choice based on plumage characteristics. As a consequence, migrant individuals with ill-matched mating traits are expected to experience a rare mating type disadvantage. Assortative mating, which is observed at several locations along the hybrid zone (Meise 1928;Randler 2007a), thus induces partial prezygotic isolation between pure parental phenotypes. Hybrids differing in plumage characteristics from both parental types will likewise experience reduced mating success (Londei 2013). The crow hybrid zone is thus characterized by (partial) prezygotic isolation between "pure" parental populations and postzygotic isolation due to reduced mating success of hybrids displaying aberrant mating traits (Meise 1928;Kryukov and Blinov 1994;Brodin et al. 2013;Vijay et al. 2016).
Plumage differences between carrion and hooded crows are encoded by two epistatically interacting, autosomal pigmentation genes ( Fig. 1) that are subject to divergent selection (Knief et al. 2019). The rest of the genome, however, appears to admix without much resistance (Poelstra et al. 2014(Poelstra et al. , 2015Vijay et al. 2016;Knief et al. 2019). Despite an in-depth understanding of the genetic architecture of the presumed mating trait, the mechanism underlying assortative mating remains elusive. In this manuscript, we focus on parental imprinting or self-referent phenotype matching, which have been suggested for crows (Brodin and Haas 2006;Londei 2013), and constitute widespread, generic learning mechanisms readily available in birds (Verzijden et al. 2012). The contribution of preference genes in close linkage to the pigmentation loci, for example, suggested for Heliconius (Merrill et al. 2019), remains highly speculative in the crow system (Poelstra et al. 2014) and will not be explored.
Here, we formulate a mathematical model invoking assortative mating as a potential driver maintaining phenotypic divergence and fit it to empirical data from the crow system. This merger of mathematical modeling with empirical data reduces the arbitrary parameter space underlying stand-alone simulations. Yet, it still allows exploring central features of general relevance within the confines of a well-investigated empirical system. It thus constitutes a complementary approach to generic simulation studies exploring the role of assortative mating on hybrid-zone dynamics and speciation (Irwin 2020). Our first aim is to investigate whether essential features of the observed hybrid zone and the available genetic data can be explained by assortative mating according to simple matching rules alone without the need to invoke natural selection. We further explore model variants in which the matching rule and the genetic architecture of the mating trait are varied. We also consider the effect of assortative mating in isolation by removing the component of postzygotic isolation. Finally, we study whether the models predict speciation for the crow system. Overall, our analyses provide detailed insight into the processes associated with assortative mating and their impact on spatial hybrid-zone dynamics.

Materials and Methods
In the following, we will first detail the empirical data and then specify the assumptions and setup of the model. We fit the dispersal parameters (see the "Dispersal" section) to data from Siefke (1994), the genotype-phenotype map to data from Knief et al. (2019) collected on adult crows in southern Europe (see the "The genetic basis of phenotypic variation" section), and use genetic data from nestlings sampled in the central-European hybrid zone to fit all other parameters of our model, including those for the mating preferences; see Figure S1 for an overview.

Phenotypic variation
Details on quantification of phenotypic variation can be found in Knief et al. (2019). In brief, the authors sampled crow individuals across the German and Italian part of the European crow hybrid zone (hereafter "central" and "southern" transects, respectively). Transects were chosen such that they were perpendicular to the geographical course of the zone and included phenotypically pure carrion and hooded crow populations (C. (corone) corone and C. (c.) cornix) at either end, as well as hybrids within the hybrid zone. In the southern transect, individuals were sampled as fully feathered adults allowing characterization of individual phenotypes. For a total of 129 individuals, 11 plumage patches were scored for the amount of eumelanin deposited in the feathers (coded as 0, 1, 2 with increasing blackness) and subsequently combined using PC analysis. PC1 explained 78.22% of the phenotypic variance with positive values reflecting a black carrioncrow phenotype, negative values a gray hooded-crow phenotype, and intermediate values being associated with phenotypic variation in hybrids (Knief et al. 2019). In this study, we used these PC1 scores as the metric for a one-dimensional representation of phenotypic variation.

The genetic basis of phenotypic variation
Using genome-wide association mapping, Knief et al. (2019) identified two major-effect loci explaining the vast majority of variation in plumage patterns represented by PC1. Most of the variation in PC1 (74%) was explained by a locus on chromosome 18 (chr18) embedded in a region of reduced recombination. For the purpose of this study, we represent allelic variation at this locus with capital letters as D = dark and L = light. The second locus, which is associated with allelic variation of the NDP gene on chromosome 1, will be denoted in small letters with alleles d = dark and l = light. Under an additive genetic architecture, (chr18) and NDP explained nearly 85% of the overall phenotypic variation. A model including recessive epistasis between the two loci explained additional 2%. It was favored by model selection procedures and will here be chosen as a default for the main model. The third locus (LRP6) explained less than 1 additional percent and will not be further considered here.
In the epistatic model, phenotypic variation is encoded as follows: allelic variation in NDP had no phenotypic effect in allblack chr18 DD individuals, but accounted for most of the residual variation in chr18 DL and chr18 LL (Figs. 1 and S2). Under this model of epistasis, the possible nine genotypic combinations (DDdd, DDdl, DDll, DLdd, DLdl, DLll, LLdd, LLdl, LLll) thus reduce to seven phenotypic states, as DDdd, DDdl, and DDll code for the same black carrion-crow phenotype. Figure 1 illustrates the phenotypic space using representative individuals with the respective genotypic constitution.
To explore the effect of genetic architecture on hybrid-zone dynamics in the model (see "Theoretical models" section), we will consider both the epistatic genetic architecture described above, and assume additivity of the two loci, which yields similar phenotypes for LL and DL genotypes (see Fig. S2), but explained overall less phenotypic variance (Knief et al. 2019).

Sample preparation and genotyping
To quantify assortative mate choice and predict its effect on the genotypic constitution across the hybrid zone, one would ideally sample breeding pairs along a transect. This is unfortunately not feasible at a larger scale in crows. Instead, we sampled nestlings from the central transect at an age at which plumage characteristics cannot be inferred with certainty (Blotzheim et al. 1993; see Figure S1). The sampling scheme thus excludes any form of postzygotic isolation arising by natural selection after the nestling stage. Due to the relatedness of nestlings, their genotypes contain information on the genotypic constitution of their parents and hence on phenotype-dependent mate choice. We therefore used the genotypic distribution of 152 nestlings sampled from 55 nests (median brood size = 3 nestlings) along the transect (Knief et al. 2019) as readout for model fitting (see the "Fitting hybridzone model parameter to genotype data" section). This procedure rests on the assumption that the genetic architecture of plumage patterning is equivalent in both transects. This is well supported by a shared evolutionary history of the parental populations at either side of the hybrid zone (Vijay et al. 2016), congruent landscapes of genetic variation across the central and southern part of the hybrid zone (Vijay et al. 2016), and near-identical selection signatures at the underlying loci (Knief et al. 2019).
We genotyped all individuals from the central transect for a selection of 1152 SNPs spread across the whole genome using the GoldenGate assay (Illumina). A detailed description of the assay design, sample preparation, SNP calling, and quality control procedure is given in Knief et al. (2019). The final dataset comprised all 152 individuals genotyped at 1111 polymorphic loci (average call rate of 99.48%). We further included 65 individuals of the allopatric populations in the GoldenGate genotyping and added five hooded crows that had been sequenced on the HiSeq2000 (Illumina) platform (paired-end libraries; sequence coverage ranged from 7.12× to 13.28×, average = 9.77×, median = 9.83×) and genotyped using the HaplotypeCaller in GATK (v3.3.0;DePristo et al. 2011;Vijay et al. 2016). In those five individuals, we found 104 inconsistencies of 5501 genotypes (1.89%). All individuals were sexed based on their heterozygosity for 114 SNPs located on the sex chromosome Z (excluding the pseudo-autosomal region located at chrZ ≤ 2.56 Mb, N = 15 SNPs). Genotypes of the pigmentation loci (see the "The genetic basis of phenotypic variation" section) were inferred as follows. The genetic factor on chr18 was represented by ancestry-informative diplotypes spanning the 2.8 MB region on chromosome 18. Ancestry (pure carrion chr18 DD ; pure hooded chr18 LL ; mixed chr18 DL ) was inferred using the NewHybrids software (v2.0+ Developmental. July/August 2007; Anderson and Thompson 2002) as described in Knief et al. (2019). Individuals that were assigned as backcrosses with evidence for recombination were excluded (N = 12 individuals in the hybrid zone and N = 11 allopatric individuals). To represent genotypic variation in NDP, we chose the SNP on chromosome 1 at 6,195,380 MB (genome version 2.5, RefSeq Assembly ID according to the National Center for Biotechnology Information: GCF_000738735.1; Poelstra et al. 2014) that explained most of the variation in plumage coloration (Knief et al. 2019) and is linked to the likely causal structural mutation (Weissensteiner et al. 2020).

Paternity assignment
Our model assumes that individuals from the same nest are fullsiblings. Thus, we estimated whether nest mates were full-or half-siblings originating from extra-pair copulations using 752 SNPs that were located on all autosomes except chromosome 18 (for details, see Knief et al. 2020). Full-and half-sibling pairs could reliably be separated based on their kinship coefficients (θ; Weir et al. 2006). For one nest containing three individuals, with two full-siblings (θ = 0.230.26) and one half-sibling (θ = 0.12), relationships could not unambiguously be resolved. We removed this nest (N = 3 individuals) and another six extra-pair young, leaving 52 nests containing a total of 131 individuals in the hybrid zone and 64 allopatric individuals for subsequent analyses (see Fig. S3).

THEORETICAL MODEL
Here we specify our main model. The model parameters were then fit to data from the crow system ("Dispersal" and "Fitting hybrid-zone model parameter to genotype data" sections). For an overview of the different components of our models and the datasets to which we fit them, see Figure S1. In the "Model variants" section, we will specify four variants of our model, which differ from the main model in their assumptions. The source code of the programs that we developed to carry out simulation studies and to fit all four model variants to the data is available from https://github.com/statgenlmu/assortative_crows/.

Temporal dimension
Similar to other European taxa, Eurasian carrion and hooded crows were presumably separated into different refugia during the last ice age (Mayr 1942;Hewitt 1988), which peaked in glacial coverage around 20,000-18,000 years ago (Frenzel 1992;Hewitt 1988). Despite a subsequent period of increasing temperature, it was likely only after intermission of a cooling phase in the younger dryas 12,800 to 11,500 years ago (Goslar et al. 2000;Rasmussen et al. 2014) that rapid warming opened up suitable crow habitat including breeding opportunities in larger shrubs or trees (Giesecke et al. 2017). We therefore set the initiation of secondary contact to approximately 12,000 years before present. Using a generation time of 6 years for the Eurasian crow (Vijay et al. 2016;Kutschera et al. 2020), this corresponds to 2000 generations.

Spatial dimension
In our model, we conceptualize space in a one-dimensional grid of 200 bins. When we fit the model to the crow data, we assume that these bins represent 5-km-wide strips that are parallel to the initial contact line (ICL) of the two color morphs, such that our model covers a range from 500 km to the west to 500 km to the east of the ICL. More precisely, we assume that bins 100 and 101 represent the areas that are directly adjacent to the ICL, up to 5 km to the west and east, respectively; bins 99 and 102 are the areas that are 5 to 10 km western or eastern of the ICL, etc. Generally, bin x for any x ∈ {1, 2, . . . , 200} represents the area that has a distance between |x − 100| · 5 km and |x − 101| · 5 km to the ICL and is west of the ICL if x ≤ 100 or to the east for x ≥ 101. The precise north-south extension of these strips is not important for our model as long as it is long enough to avoid boundary effects in our dispersal model ("Dispersal" section).

Modeling approach for assortative mating
We assume that each bin is populated by the same number of crows and that this number is so large that we can neglect genetic drift and any other source of random fluctuations in genotype frequencies, which means that our model is deterministic. For each genotype g with respect to the two mating-trait loci, let f x,g (t ) ∈ [0, 1] be the frequency of g in bin x at time t. We assume discrete time steps s = 0, 1, 2, . . . corresponding to nonoverlapping generations, with time t = 0 referring to the time of the secondary contact of carrion crows and hooded crows after the ice age. We assume that at time 0, all crows in bins 1-100 were carrion crows and all crows in bins 101-200 were hooded crows. For x ≥ 101, this implies f x,g (0) = 1 if g = LLll, and f x,g (0) = 0 otherwise. For x ≤ 100, we obtain f x,g (0) = 0 for g ∈ {DDdd, DDdl, DDll}, and assume Hardy-Weinberg frequen- , where δ is the initial frequency of allele d in bin x, which we assume to be equal in all bins x ≤ 100 at time 0.
The core of our model is the recursion to compute the generation t + 1 genotype distribution matrix f (t + 1) = ( f x,g (t + 1)) x,g from the genotype distribution matrix f (t ). Implicitly, we assume that the sex ratio is the same for all bins and all genotypes. Following a common mass-action law approach (see, e.g., Otto and Day 2011), we assume that the number of matings between individuals with genotypes g 1 and g 2 taking place in bin x is proportional to the number of potential couples in bin x, weighted by a mating-preference value w g1,g2 (with w g1,g2 = w g2,g1 ; see "Mate choice" section). The latter can be interpreted as the probability that a female and a male crow with genotypes g 1 and g 2 will mate if they meet. Note that in this setup it is not relevant whether the female, the male, or both sexes are taking the mating decision based on the similarity of their phenotypes.
If d z,x is the dispersal probability from bin z to bin x ("Dispersal" section), the number of matings between a female of genotype g 1 and a male of genotype g 2 taking place in bin x is proportional to where the summation variables y and z iterate over all combinations of bins from x − 10 (or 1 if x − 10 ≤ 0) to x + 10 (or 200 if x + 10 > 200), with the range of ±10 bins reflecting a maximum dispersal distance of 50 km ("Dispersal" section). This proportionality could arise, for example, if within each bin randomly moving individuals form a pair at the beginning of the mating season and occupy one of a limited number of territories to reproduce, as is observed in the field by Blotzheim et al. (1993). Individuals that need to search longer for suitable mating partners have smaller chances of finding a vacant spot for their nest. Thus, a form of sexual selection emerges even though we assume monogamy. In accordance with this heuristic consideration, our model implies that individuals that are surrounded by fewer preferred mating partners will on average have fewer offspring. Thus, assortative mating induces a reduction in fitness via frequency-dependent selection in both sexes.
To finally obtain the recursion to compute f (t + 1), let for each triplet (g, g 1 , g 2 ) of genotypes be μ(g | g 1 , g 2 ) the probability that an offspring of a mother with genotype g 1 and a father of genotype g 2 has genotype g, according to Mendelian genetics with free recombination between the two loci. The frequency of genotype g in nestlings of generation t + 1 in bin x is then where the value of c x (t ) ensures that g f x,g (t + 1) = 1.
For the above, sum note that if the summation indices g 1 and g 2 refer to different genotypes, for example, g 1 = DDdd and g 2 = DDdl, then also the inverse combination, g 1 = DDdl and g 2 = DDdd, appears in the summation with the same bins y and z, and the two corresponding summands contribute the same value to the sum.

Dispersal
Natal dispersal distance was assumed to be equal for both sexes and parameterized with data from Siefke (1994 , Table 3 and table caption) who analyzed data of hooded and carrion-crow individuals that were ringed as nestlings and recovered during the breeding period. In brief, Siefke (1994) compiled ring recovery data from 198 individuals that were ringed as nestlings. Of these, 89 were recovered during the potential breeding period of consecutive years (April-August) and classified as likely cases of natal dispersal. Dispersal distances are given as Euclidian distances between the location of birth and recovery (putative breeding location). Kalchreuter (1970) investigated ring recovery data of carrion crows and obtained a similar distribution of dispersal distances. To use the dispersal distances quantified by Siefke (1994) on a two-dimensional map, we first needed to project these distances to our one-dimensional representation of the contact zone. For this we assume that the ICL is a straight line going in south-north direction and decompose each dispersal movement of a crow into its component D x orthogonal to the ICL and its component D y parallel to the ICL. The values of D x or D y are positive, if the direction of movement is west to east or south to north, respectively, and negative otherwise. We assume that the random distribution of the dispersal vector (D x , D y ) is a mixture of two spherical two-dimensional normal distributions. A single spherical normal distribution did not lead to a satisfying fit to the long-tailed distribution of the dispersal distance data (Fig. S4, top). Thus, the dispersal vector can be represented as (D x , D y ) = S · (Z x , Z y ), where S is a random variable with two possible values σ 1 and σ 2 and the vector (Z x , Z y ) has a two-dimensional standard normal distribution. To fit the three parameters σ 1 , σ 2 and p = Pr(S = σ 1 ) to the empirical dispersal distances (Siefke 1994) we numerically maximized the likelihood of (σ 1 , σ 2 , p) using the method of Byrd et al. (1995) as implemented in the R command optim with the option L-BFGS-B (R Core Team 2018). To compute likelihoods, we used that with probability p the squared rescaled dispersal distance (D x , D y )/σ 1 2 or other- x + Z 2 y and thus chi-square distributed with 2 degrees of freedom. The fitted values are σ 1 = 6.07 km, σ 2 = 28.12 km, and p = 0.767.
As the bins in our models represent areas that extend as long stripes parallel to the ICL, only the component D x is relevant. Thus, for the probability that a crow stemming from bin x is found in bin x + d, the dispersal distance d is a discretization of D x , whose distribution consists of 76.7% of a normal distribution with a standard deviation of σ 1 = 6.07 km, with the other 23.3% coming from a normal distribution with a standard deviation of σ 2 = 28.12 km. For computational efficiency we restrict the range of d to ±10 bins, corresponding to a maximum ±50 km, setting Pr(d = −10) = Pr(d = 10) = Pr(D x < −47.5) = Pr(D x > 47.5) (Fig. S4, bottom). Furthermore, the boundaries of our binned space are reflecting. That is, probability weights of dispersal d that would result in x + d < 1 or x + d > 200 are shifted to instead increase the probabilities of bin 1 + |x + d| or 200 − |x + d − 201|, respectively.

Mate choice
In the absence of decisive knowledge on the mechanism underlying mate choice in the crow system, and for computational feasibility, we assume self-referent phenotype matching in our main model. Using alternative models (see "Finite-size population models with self-referencing or imprinting," "Effect of genetic drift," and "Imprinting versus self-reference" sections) we further explore whether this simple matching rule can serve as a proxy for maternal or paternal imprinting (Verzijden et al. 2005;Verzijden et al. 2012), as has also been suggested for the crow system Haas 2006, 2009;Londei 2013). We summarize the phenotype by the values of the first PC derived in Knief et al. (2019) (see "Phenotypic variation" section and its encoding by two-locus genotypes as specified in "The genetic basis of phenotypic variation" section). For simplicity, phenotype values are scaled to a range between 0 (pure hooded crow phenotype, genotype LLll) and 1 (pure carrion-crow phenotype black, genotype DDdd, DDdl, DDll).
Let ϕ be the genotype-phenotype map. The matingpreference value w g1,g2 then is the mating probability of two birds with phenotype values ϕ(g 1 ) and ϕ(g 2 ) relative to the mating probability of two birds of equal color, such that w g1,g2 = 1 if ϕ(g 1 ) = ϕ(g 2 ). We model the mating-preference function w with two parameters η and ζ , where η is the mating preference of crows with very different phenotypes and ζ controls the width of a Gaussian kernel modeling how mating preference w g1,g2 of more similar birds depends on their phenotypic difference ϕ(g 1 ) − ϕ(g 2 ): Figure S5 shows how w g1,g2 depends on ϕ(g 1 ) − ϕ(g 2 ), with parameter values η and ζ set according to maximum-likelihood values in our main model. Examples for resulting mating-preference values for various parameter settings are shown in Figure 2B.
Fitting hybrid-zone model parameter to genotype data Weissensteiner et al. (2020) inferred that the d allele of the NDP gene arose approximately 500,000 years ago and segregated in the ancestor of both carrion and hooded crows. As the three genotypes DDll, DDld, and DDdd lead to the all-black phenotype, we therefore allow that the ancestral frequency of allele d in the west can initially be different from 1. This initial frequency δ of the d allele is the third model parameter besides ζ and η. A fourth parameter is the geographic location of the ICL. For this, all sam-pling sites are projected on a one-dimensional transect as specified in Knief et al. (2019), see also the "Sample preparation and genotyping" section. The position on the transect is measured in kilometers, starting with the most western sampling site as km 0. For the position of the ICL (which is assumed to intersect the transect in a right angle) we allow the range from km 300 to km 700 on the transect.
To numerically optimize all four parameters, we first simulated the hybrid-zone model for any triple of candidate values for the first three parameters. Then we calculated for the grid km 300, km 301, ..., km 700 the probability of the empirical genotype data for the case that the ICL was at that position. To calculate the parameter likelihoods, that is the probability of the data, we rounded the nest locations on the transect to the bins and calculate for each possible combination of parental genotypes the probability that the chick genotypes of a nest stem from such parents (assuming Mendelian inheritance and no segregation distortion as Knief et al. 2020, found no evidence for biased sex ratios in F1hybrids or for postzygotic isolation). The likelihood calculations were then based on the simulated genotype distribution after 2000 generations and mate choice according to the model assumptions. Thus, we optimized the likelihood of the four parameters with respect to the genotype data. For this we used the parallel version of the L-BFGS-B method provided by the R (version 3.5) package optimParallel (Gerber and Furrer 2018). To infer the point of inflection of a cline from allele frequencies y(x) in all bins x ∈ {1, . . . , 200} we first chose the bin x with the largest value of |y(x − 1) − y(x + 1)|, then fitted a third order polynomial to y in the range from x − 2 to x + 2 and calculated the point of inflection of the fitted polynomial.
To assess how well our model is in agreement with the empirical data, we simulated 1000 datasets according to our fitted main model and each time calculated the log-likelihoods. We found that the log-likelihoods of the real data is well within the range of the simulated log-likelihoods, see Figure S22.

MODEL VARIANTS
In the following, we formulated deviations from our main model exploring several parameters we deem of importance in the crow system and beyond. These include the specific mode of how phenotypic differences are translated into mating-preference function ("Categorical choice model" section), the impact of the genetic architecture of the mating trait ("Additive model" section), how the system is affected if we remove sexual selection altogether (and thereby postzygotic isolation) while maintaining assortative mate choice ("Neutral model" section), and how the nature of the matching rule affects hybrid-zone dynamics ("Finitesize population models with self-referencing or imprinting" section).

Categorical choice model
In this variant of our model, hybrids are summarized into a single category and have no mating preferences. However, carrion crows prefer carrion crows and hooded crows prefer hooded crows as mating partners. If two crows interact and both are of the same category-carrion crow, hooded crow, or hybrid-their matingpreference value w is 1. If one is a hybrid and the other is a hooded crow or a carrion crow, their mating-preference value w is given by the parameter ψ 0 or ψ 1 , respectively. If one is a hooded crow and the other a carrion crow, their mating-preference value w is the product ψ 0 · ψ 1 . We combined this model with our standard genotype-phenotype map as specified for the main model ("The genetic basis of phenotypic variation" section).

Additive model
This model variant assumes an additive genetic architecture of plumage pigmentation patterns. The corresponding genotypephenotype map was combined with the mating-preference function of the main model ("Mate choice" section). We constructed the additive genotype-phenotype map by taking the same phenotypic values for the LLll, LLdl, and LLdd genotypes as in the epistatic model and the slope between LLdl and DLdl to infer all other additive phenotypic values (Fig. S2). Note that according to this model, only genotype DDdd leads to an entirely black carrion-crow phenotype. As we assume that shortly after the ice age all crows in the west were black, the initial allele frequency of the d allele at locus 2 was 1.0 (see Supporting Information Section G.3; model variant with a relaxation of this assumption).

Neutral model
Assortative mating according to matching rules usually induces frequency-dependent sexual selection reducing the fitness of rare phenotypes and can thus contribute to hybrid-zone stability (Kopp et al. 2018;Irwin 2020). If, however, in the range of the hybrid-zone darker crows tend to find mating partners in the west and lighter crows tend to mate more in the east, this may also contribute to the maintenance of clines or at least slow down their disappearance. To examine whether these migratory effects of mate choice alone can explain the observed clines in the genetic data or whether induced sexual selection is crucial for this, we define a "neutral" model, that is, without induced sexual selection or any other form of selection. In this model, sexual selection is eliminated by compensating the fitness cost of preferring rare phenotypes as mating partners, such that all individuals have the same expected number of offspring. Rare phenotypes do not incur any costs if they do not find a mating partner upon dispersal, but can keep on searching without additional costs to find a match. In essence, this version of the model removes the component of postzygotic isolation that otherwise arises as a consequence of frequency-dependent sexual selection.
Given the narrow width of the hybrid zone, costs for surveying larger distances for this highly mobile species appear minor. Factoring out the potential relevance of individual-based recognition and local dominance hierarchies for reproduction in corvids (Chiarati et al. 2010;Holtmann et al. 2019), the assumption of the neutral model thus seems not too far-fetched. Accordingly, crows stemming from areas in which their phenotypes are rare will tend to ultimately mate in areas in which the frequency of their matching mating partners is higher. Thus, effective migration patterns are not symmetric and lead to slight differences in population densities that depend on the phenotype. Therefore, to avoid induction of sexual selection by local competition, we have to assume that local population densities can vary keeping only the global population size constant. This means that the local population densities are not governed by carrying capacities but only by the distribution of mating opportunities of crows according to their mating preferences. This may run counter to ecological intuition, but we make this assumption to obtain a neutral model that we can compare to our main model in which selection is induced. As a consequence, the neutral model posits that the total number of individuals in a bin can vary among the bins, and bins of higher density will produce more emigrants in the next generation. We combined this model with the standard genotype-phenotype map of the main model ("The genetic basis of phenotypic variation" section).
For the precise specification of this model, let f x,g (t ) (still) be the frequency of g in bin x at any time t relative to the initial total frequency g f x,g (0) in bin x (or any other bin), but note that, different from the main model, the sum g f x,g (t ) is in general not 1 for t ≥ 1. The initial matrix f (0) is the same as in "Modeling approach for assortative mating" section. The recursion to calculate f (t + 1) is now where the factors v x,g (t ) compensate induced fitness effects for individuals stemming from bin x and having genotype g, which could reflect, for example, that these individuals intensify or extend their search for mating partners by a factor of v x,g (t ), with the assumption that this is possible without incurring any fitness cost. For this, the values of v ·,· (t ) must simultaneously for all (y, g 1 ) fulfill which we solve numerically (see Supporting Information Section F for details).

Finite-size population models with self-referencing or imprinting
In the models specified above, we neglect random variation such as genetic drift, which is clearly present in real populations of finite size. To assess how robust our results are against random fluctuations of this kind, we performed simulations using individual-based models following a standard Wright-Fisher approach with Poisson-distributed offspring numbers with an average number of surviving offspring of two per nest. In brief, the Central European crow habitat is idealized as a two-dimensional grid of 500 × 1000 territories, each spanning 1 km 2 and containing one crow nest (which lies within the large empirically documented range; Blotzheim et al. 1993). These models consider only genetic drift and do not account for other sources of randomness and spatial or temporal variation of ecological conditions. The same assumption on crow dispersal and mating preferences are used as in the main model (see Supporting Information Section J for details).
The individual-based model approach not only allows comparing the results from the main model to a finite population size model with a self-referent matching rule (Supporting Information Section J.1). Its inherent flexibility also opens the opportunity to examine the effect of imprinting in which females prefer males that resemble the average phenotype of her parents (Supporting Information Section J.2), males that resemble their fathers or, in another variant of the model, their mothers (Supporting Information Section J.3).

MAIN MODEL: PARAMETER ESTIMATES
The main model assumes an epistatic genetic architecture of the mating trait, and sexual selection induced by assortative mating favoring common phenotypes. This model produced a good fit to the empirical genotype frequency distribution, and with a log-likelihood of −175.6 was superior to the alternative models discussed next. (Note that these models have the same numbers of parameters.) The cline for locus 1, which is responsible for the majority of the overall phenotypic variation, was the steepest (Fig. 2). It was closely followed by a shallower, and geographically offset cline of the second, epistatically interacting locus. This is concordant with the clines that have been fitted to the empirical data in Knief et al. (2019). At the most western end of the transect, 34 single individuals were sampled from different nest, such that the allele frequencies shown in Fig. 2B in the west are all 0, 0.5, or 1. The local genotype frequencies are shown in Figure S3.
Consistent with the proposition that the light allele of the second locus already segregated in the ancestral population of carrion crows (Weissensteiner et al. 2020) the maximum-likelihood parameter estimate of the initial allele frequency of the dark allele in carrion crows was below 1 (0.777). As maximumlikelihood parameter values for the mate preference function we found (ζ , η) = (2.17, 0.617) and kilometer 302 on the transect as position of the ICL. Individuals differing most extremely in appearance (i.e., pure hooded crows with genotype LLll vs. pure carrion crows with genotypes DDdd, DDdl, or DDll) were only 61% as likely to choose each other for reproduction compared to pairs of individuals of equal phenotypes. Note that the ultimate mating rates not only depend on this preference matrix but are also contingent on the frequency of each genotype in the local population determining the relative frequency of mutual encounter.

HYBRID-ZONE MOVEMENT
The maximum-likelihood parameterization of the main model predicted a substantial shift of the hybrid-zone center to the east after 2000 generations. More precisely, the inflection points of loci 1 and 2 were inferred to be 189.6 and 188.3 km to the east of the ICL. To further explore these dynamics, we simulated 6000 generations under the main model that we parameterized with the maximum-likelihood estimators described above. During the first ∼1200 generations, the zone moved eastward from the ICL (in favor of dark morphs). After allele d of locus 2 decreased to a certain frequency, the hybrid zone began to reverse its movement and shifted westward (now in favor of the light morphs, Fig. 3 and Supporting Information Section I.1).
The movement of the hybrid zone can be explained by local selection pressures on the alleles of the two mating-trait loci that are induced by the mating preferences as follows. Due to the epistatic architecture, three different genotypes code for allblack crows (Fig. 1). This induces a high initial frequency of dark morphs among hybrids, which induces positive frequencydependent sexual selection in favor of darker crows. This, in turn, increases the frequency of alleles D and d in the hybrid zone (Fig. 3, arrows A and B), which promotes a shift of the hybrid zone to the east (arrow C). At the same time, the l allele moves into the west. Although allele L is under negative selection in the west, allele l has no fitness effects. It freely introgresses into the DD background, which entails a decrease in allele d (arrow D). In the long run, this also decreases the frequency of d in the proximity of the hybrid zone (arrow E), which reduces the fraction of dark morphs among the hybrids and thus the fitness of dark crows. This eventually brings the hybrid-zone movement to a halt after approximately 1200 generations. As the decrease in allele d continues (arrows D and E), sexual selection provides an advantage to lighter morphs, decreasing the frequencies of D and d in the hybrid zone (arrows F and H). As a consequence, the hybrid zone travels westward (arrow G) and finally drives alleles d and D to extinction after 5500 generations.

Categorical choice model
With a log-likelihood of −199, the model fit the data substantially worse than the main model assuming mate discrimination between all seven phenotypes. In this model, mate choice relevant phenotypes were categorized into pure and hybrid type. Maximum-likelihood estimates for the mating-preference parameters ψ 0 and ψ 1 of hooded and carrion crows, respectively, were 0.318 and 0.645, which translates into a larger difference in mate choice (Fig. 2D). As in the main model, pure morphs had a clear preference for phenotypically matched partners. The probability of mating between carrion and hooded crows, however, was only 21% and thus only about a third of that in the main model (61%). According to the best-fitting model parameter values, also for lo-cus 2 the dark allele d was fixed in the west before secondary contact.
In terms of spatio-temporal dynamics according to this model, the hybrid zone had shifted westward since the initial contact. The inferred inflection points of clines of loci 1 and 2 after 2000 generations were located 142.1 and 141.7 km to the west of the ICL, respectively. According to simulations with more than 2000 generations, the model predicts a continuous westward shift of the hybrid zone until carrion crows go extinct (Fig. S15).

Additive model
Changing the genetic architecture from epistasis between the mating-trait loci to additivity had a strong impact. The best possible fit to the data for the additive model was poor (log-likelihood −228.7), and maximum-likelihood estimates of ρ = 2.878 and η = 0.987 implied only very slight mating preferences overall (Fig. 2). Mating-preference probabilities only differed by 1.3% (0.987 vs. 1.0) indicating very weak sexual selection. It is thus not surprising that clines of both mating-trait loci were shallow. Still, the difference of 1.3% in mating-preference probabilities shielded against gene flow, and clines were more pronounced than by simple diffusion of a neutral locus (see the "Effects of assortative mating on unlinked neutral loci" section). Furthermore, due to additivity the frequency curves of the two alleles were still almost identical after 2000 generations (Fig. 2). With inferred points of inflection of both clines only 0.1 km to the east of the ICL, almost no shift of the clines was predicted. In summary, the additive model could not explain some essential features of the empirical data, highlighting the importance of the genetic architecture of mating traits for hybrid-zone dynamics.
The above results are based on the assumption that, initially, all crows in the west were black, which entails that the dark alleles D and d were fixed in the west. This differs from the main model, where the frequency of the d allele was free to vary. In Supporting Information Section G.3, we report results of further parameter combinations of the additive model relaxing the assumption that all parental carrion crows were entirely black upon secondary contact. Also with this assumption, clines were either not steep enough to fit the empirical data, or the empirical disparity in cline location and inclination between the loci was absent. Overall, none of the models with additive genetic effects fit the data well.

Neutral model
Next, we formulated a model without sexual selection exploring the effect of mate choice alone. All individuals had the same expected number of offspring and did not incur any fitness advantage or disadvantage in relation to their phenotype. With a loglikelihood of −210.3 the fit to the data was substantially worse than in the main model including sexual selection. The bestfitting parameter combination for the neutral model had an initial d allele frequency for locus 2 of 0.619, and showed a slight shift for both loci 1 and 2 of 11 and 7.2 km to the east from the ICL (measured at the inflection points of the clines).
The maximum-likelihood parameter estimates of the mating function ζ = 5.53, η = 0.0237 translate to the most extreme avoidance of non self-phenotypes of all investigated model alternatives. Matings of the opposite morph (carrion with hooded crows) were 50 times less likely than between equivalent phenotypes (Fig. 2). Similar to the categorical model, this may be explained by the requirement of maintaining a narrow hybrid zone over 2000 generations against increased gene flow in the center of the zone by hybrids traveling far to seek an appropriate partner (with no cost in this model).
Although the model somewhat mimicked the empirical disparity in allele frequency clines between loci, it did not capture the steepness of the cline for locus 1. A steeper cline could be achieved by setting the model parameters to values reflecting stronger assortative mate choice. In this case, however, the frequency curve for locus 2 was shifted too far to the east. This would produce many LLdd individuals in the hybrid zone, which were near-absent in the empirical data, and explains the poor fit (compare Fig. S7 to Figs. S8 and S9). Interestingly, the model predicted a decrease in population density in the center of the hybrid zone accompanied by a population increase at a distance of around 200 km on either side (Figs. S7 and S9,bottom). Assuming no cost of mate choice and movement, this shift in densities is due to individual compensation for differences in sexual attractiveness induced by the preference function. As we assume that individuals of a phenotype that is rare in their environment intensify their search for mating partners, these individuals may effectively migrate over larger distances and deplete local population densities.

NEUTRAL LOCI
Using the best-fit parameters of the main model, we ran additional simulations in which we included a third, neutral locus. This biallelic locus was unlinked to the loci associated with the mating trait and had no effect on mate choice or dispersal. We defined the focal allele as the allele with initial frequency of 1 in the west and of 0 in the east. We compared this simulation to a simulation in which all three loci had no effect on the mating trait or dispersal (random mating) and essentially follow neutral diffusion according to our dispersal model ("Dispersal" section). (Note that in the latter case, locus 1 and 3 have the same allelic distribution.) The spatial distribution of alleles from the neutral locus 3 was affected by the mating-trait loci, but this effect was very weak (Figs. 4 and S16). After 2000 generations, the frequency of the neutral allele at locus 3 was 0.620 at the western boundary of the range and 0.398 in the far east in the main model with locus 1 and 2 under sexual selection. With random mating, the frequency of locus 3 was slightly closer to equalized allele frequencies (at 0.5) with 0.597 in the west and 0.403 in the east. Hence, this model predicts that sexual selection does not substantially reduce effective migration at loci unlinked to those determining the mating trait. Even if assortative mating in this model entails that hybrids have reduced fitness, local reduction in gene flow of the mating trait affected other parts of the genome only marginally.

EFFECT OF GENETIC DRIFT
The finite population size models allow for stochastic variation in the spatial distribution of genotypes (Fig. S17), and ultimately cline dynamics (Fig. S20). The finite-population model with the self-referent matching rule, parameterized according to our main model, showed the same type of hybrid-zone cline dynamics as our main model ( Fig. 5B; Supporting Information Section J.1 and Fig. S18). Hence, the outcome was not qualitatively altered by genetic drift.

IMPRINTING VERSUS SELF-REFERENCE
The flexibility of the individual-based models further allowed exploring the effect of the matching rule. Specifically, we assessed whether self-referent phenotype matching would serve as a good proxy for imprinting, as has been suggested before (Verzijden et al. 2005;Verzijden et al. 2012). Imprinting was modeled as female-based mate preference based on the phenotype of one parent (Supporting Information Section J.3) or on the average phenotype of her parents (Supporting Information Section J.2). In the first case, results closely aligned with the outcome based on self-referencing, regardless of which sex served as the imprinting model (Figs. 5B and S20). In case of imprinting on the mid-parental phenotype, cline dynamics were substantially different, Here, the hybrid zone moved eastward until the carrion-crow phenotype became fixed (Fig. 5C). Another differences was that the locus-2 cline was shifted slightly to the east compared to the locus-1 cline, such that directly in the hybrid zone after 1000 generations the locus-2 dark allele was more frequent than the locus-1 dark allele (Fig. S19).

Discussion
In this study, we explored the effects of mate choice on hybridzone dynamics. Specifically, we studied how the interplay of a genetically determined mating trait and preference for the trait based on simple matching rules (self-referencing, imprinting) affects hybrid-zone position and width through time. Analyzing empirical data from the European crow system by simulation-based model fitting allowed us to estimate the strength and mode of mate choice as a function of individual genotype combinations. Clearly, every theoretical model is a necessary simplification of reality (Popper 1934;Box 1979). We chose to neglect ecological factors and did not explicate spatial details of dispersal barriers despite known effects on hybrid-zone dynamics (Barton 1979;Barton and Hewitt 1989). We rather focused on the question of whether intrinsic behavioral properties of the system suffice to explain the phenotypic and genotypic composition across the crow hybrid zone for which rich data have now been collected for nearly a century (Meise 1928;Knief et al. 2020). We reach the conclusion that assortative mating based on learning of phenotypic cues can explain empirical observations to a great extent under realistic assumptions and model parametrization. For a better understanding of the underlying processes, we examined the effects of the mating-preference function, the genetic architecture of the mating trait, the fitness costs induced by frequency-dependent selection and finite population size using model variations. We discuss these aspects in turn with a focus on their relevance for reproductive isolation. First, we home in on the mating-trait loci themselves and then consider genome-wide implications.

MATING-TRAIT LOCI
Theory predicts that mate choice based on matching rules, such as self-referent phenotype matching or imprinting, readily enhance positive frequency-dependent sexual selection of locally abundant mating-trait alleles (Yeh and Servedio 2015). In the context of secondary gene flow between populations with diverged mating-trait values, this process is expected to result in directional selection on each side of the hybrid-zone pulling in opposite directions (Servedio 2016). Hence, sexual selection mediated by assortative mating may promote hybrid-zone stability in the form of geographic clines (Brodin and Haas 2009;Irwin 2020) or as mosaics, depending on the dispersal kernel (M'Gonigle and FitzJohn 2010). These theoretical predictions are in accordance with empirical observations. Geographically structured polymorphism in phenotypes with relevance to mate choice are not uncommon in nature (McLean and Stuart-Fox 2014), and learning of (self-) recognition cues or imprinting is widespread in the animal kingdom (Pfennig et al. 1983;Ten Cate and Vos 1999;Bereczkei et al. 2004). Matching-rule based mate choice thus features as a prevalent force in maintaining or promoting divergence with ongoing gene flow (Irwin and Price 1999;Verzijden et al. 2012). This prediction applies particularly well to the context of hybrid zones in which mating traits that diverged in isolation come into secondary contact (Randler 2008). Consistent with these predictions, our main model, best fitting the empirical data in the crow system, could explain the maintenance of steep clines over thousands of years for both loci associated with the mating trait (Fig. 2). However, all our models predicted that the hybrid zone will ultimately vanish, either as one color morph takes over, as in our best-fitting model, or as the clines will flatten out.

The nature of mate choice
The evolutionary outcome of sexual selection was contingent on the mating-preference function. Assuming self-reference for all possible seven phenotypes (Fig. 1) in the main model required little variation in mating probabilities (1.7-fold difference between minimal and maximal mating probabilities) to maintain high allelic differentiation between the parental populations. A different outcome was predicted for mate choice based on categorizing crows into three classes of pure parental and hybrid forms, as has previously been suggested for the crow system Haas 2006, 2009). In the categorical model, the lack of discrimination among hybrid phenotypes increased gene flow of alleles from both loci with hybrids serving as a bridge for introgression between the pure parental forms (Irwin 2020). Stronger avoidance (4.8-fold difference) of phenotypically dissimilar matings in the areas dominated by the parental genotypes was necessary to still obtain the empirically observed steep cline in genotypes and alleles of the mating-trait loci. This result demonstrates that the details of the preference function matter, even for a simple trait within a single sensory modality. This has potential implications for the evolution of matching-rule based mating preferences.
In the case of imprinting, many different modes are imaginable, including imprinting of hybrids to only one parental phenotype, to both, to intermediate phenotypes, or to the phenotypes of siblings or even other conspecifics of the local population. Imprinting to parental phenotypes may have similar effects as self-referencing as both result in positive frequency-dependent selection on the mating trait (Kopp et al. 2018). Yet, in simulations addressing this question, we observed that cline dynamics depended strongly on the precise assumptions of imprinting ("Imprinting versus self-reference" section). For the model with imprinting on one parent, cline dynamics mirrored the predictions from our main model with a simple self-referent matching rule. Simulations in which females imprinted on the mid-parent phenotype led to substantially different cline dynamics.
Here, we only explored a small proportion of possible mate choice schemes focusing on simple-self referent matching, which constitutes a reasonable first approximation in birds. Others models, explicitly addressing the interplay between the mating trait and the optimal resolution of preference for the trait seem to be worthy of further exploration. Theoretically, assessing the effect of general avoidance or preference for non-parental phenotypes would, for example, be an interesting avenue to pursue. Empirical investigation of the details underlying the matching rule and the actual mating cues is unfortunately near impossible, as it would require prohibitively large-scale and long-term crossfostering experiments in the crow system.

Genetic architecture of the mating trait
Hybrid-zone dynamics were strongly contingent on the genetic architecture of the mating trait. Gene flow was most readily reduced under an epistatic architecture, which fits the data best. Assuming an additive genetic architecture, we found no parameter combination that would recover cline steepness and offset of cline centers characteristic of the empirical data. Instead, under the best-fitting parameter values, the parental populations were rapidly homogenized by gene flow as illustrated by flattening of the clines across large geographic distances. Brodin and Haas (2009) modeled the southern-Danish hybrid zone between carrion crows in the south and hooded crows in the north. They assumed an additive model of three unlinked loci for the genetic basis of the mating trait without any other selection than assortative mating induced sexual selection. Assuming ongoing influx of hooded crows from the north and carrion crows from the south, their simulation runs led to stable phenotypic clines that resemble the cline observed in the field. However, with additivity of the mating trait it generally appears to be difficult to maintain divergence, as is illustrated by simulations from Irwin (2020) where very strong assortative mating was required for hybrid-zone maintenance (10-fold difference). In these simulations, an increasing number of loci required even stronger assortment. These findings highlight the importance of the genetic architecture underlying the mating trait. Tentatively, these results may suggest that matching rules may be most relevant for evolution where single sensory modalities with a simple, at most oligogenic architecture (such as color; Cuthill et al. 2017), are recruited for mate choice (Nadeau et al. 2007). Whether epistasis generally promotes divergence remains an open question and warrants further theoretical exploration. On the empirical side, genetic investigation of mating traits in systems where preference learning has been shown to play an important role for prezygotic isolation provides a promising avenue (Yang et al. 2019). This likewise applies to sexually selected traits that appear to accelerate speciation rates on macroevolutionary scales (Panhuis et al. 2001;Hugall and Stuart-Fox 2012;Seddon et al. 2013).
The genetic architecture also had a major effect on the spatio-temporal dynamics of the hybrid zone. Hybrid-zone movement is traditionally expected to be induced by dominance or constant asymmetries in allelic fitness, usually neglecting frequency dependence (Mallet 1986;Secondi et al. 2006;Brodin et al. 2013) (for a mock exploration of our data with the model proposed by Barton 1979, see Supporting Information Section H.1). The effect of epistatis on hybrid-zone movement has rarely been explored, least in the context of sexual selection. In this study, we observed that cline centers remained at the initial line of contact in the additive model, but shifted in all models when introducing an (empirically motivated, Knief et al. 2019) epistatic architecture. The combination of epistasis in the mating trait and frequencydependent fitness effects had interesting emergent properties. The best-fit main model predicted an eastward movement of the hybrid zone which reversed its direction after ∼1200 generations. This shift was due to changes in the relative allele frequencies of the two interacting loci. The Central European crow hybrid zone was mapped for the first time with great precision in 1928 (Meise 1928). Reexamination of its location by Haas and Brodin (2005) 80 years (or ∼13 generations) later suggested a slight movement toward the east. Although empirical sampling variance and sensitivity to model assumptions preclude a one-to-one comparison, both our results and those of Haas and Brodin predict that the crow hybrid zone does not conform to an equilibrium or quasiequilibrium state and instead is dynamic to the present day. According to our model predictions, this will ultimately lead to the extinction of one or the other color morph within the time frame of a glacial cycle.

Removing postzygotic isolation
In our main model, assortative mating induced prezygotic isolation between pure carrion and hooded crows. Frequencydependent positive selection likewise induced postzygotic isola-tion in the form of a rare mating type disadvantage of hybrids. To exclusively consider the effects of assortative mating alone, we formulated a neutral model removing any induced fitness effects of sexual selection. Individuals still mate associatively following a simple matching rule, but do not suffer from any costs in finding a suitable partner. A possible interpretation of this model is that there is a global, but no local carrying capacity, breeding sites are not strongly limiting, and the time to find mating partners has no fitness effect. As a consequence, hybrid back-crosses with a phenotype resembling the pure morphs tend to leave the hybrid zone-darker crows to the West, crows with lighter coat to the East. To make this model entirely neutral, we allowed for migration of individuals in search of appropriate mating partners without any constraints due to local population densities. Bins of higher density will produce more emigrants in the next generation, which leads to higher population densities toward the parental ranges and lower population density in the hybridzone center.
Although the neutral model overall poorly supported the empirical geographic distribution of genotypes, it has interesting conceptual implications. In the popular tension-zone model, where hybrid zones are maintained at an equilibrium of migration and selection against hybrids, population density is preset by the local carrying capacity (Barton and Hewitt 1985). Moving tension zones are predicted to be trapped in regions of low carrying capacity or biogeographic barriers to migration (Barton 1979;Barton and Hewitt 1985). On the contrary, our model predicts variation in population density as an emergent property of assortative mating in the absence of any form of selection. Hence, reduced population densities, as observed or predicted for hybrid zones (Barton andHewitt 1981, 1985;Hewitt 1988), need not necessarily result from ecological constraints. The somewhat artificial assumption in our neutral model that population density variations do not have any fitness effects, may be relevant for situations where carrying capacity is not (yet) reached.

RELEVANCE TO SPECIATION
The genic view of speciation emphasizes that reproductive isolation will initially be caused by a small number of loci (Wu 2001;Presgraves 2010;Wolf et al. 2010). Divergent or disruptive selection reduces effective migration of alleles at these loci allowing for local allelic differentiation (Ravinet et al. 2017;Wolf and Ellegren 2017). The mating-trait loci considered above constitute such barrier loci subject to selection (direct barrier effects). However, with few barrier loci of moderate effect, reproductive isolation remains confined to the linked local genomic neighborhood. Gene flow of neutral genetic variants disassociated from the selected background by recombination remains high. Indirect barrier effects reducing genome-wide gene flow are only expected to unfold when a sufficient number of barrier loci come into linkage disequilibrium-depending on the strength of selection and the recombination rate (Barton 1983;Bierne et al. 2011;Feder et al. 2014). Here, we explored this indirect barrier effect by introducing a third, unlinked locus serving as a proxy for genome-wide, neutral genetic variation. Adding this locus to the best-fit main model allowed us to quantify genome-wide reduction in gene flow elicited by frequency-dependent sexual selection on a mating trait. Although gene flow was slightly reduced, the effect was overall very small. This is consistent with theoretical models predicting that the effect of selected loci on neutral loci is small if the rate of recombination between these loci is large relative to the selection strength (Barton 1986;Barton and Gale 1993;Barton and Shpak 2000a;Durrett et al. 2000;Sedghifar et al. 2016). Likewise, simulations by Irwin (2020), found that assortative mating had only minor effects on unlinked, neutral loci.
Speciation can be conceptualized as the built up of linkage disequilibrium between loci subject to divergent selection (Felsenstein 1981). Eventually, this requires coupling of alleles across many barrier loci contributing to reproductive isolation (Ravinet et al. 2017;Wolf and Ellegren 2017). In the crow system, however, it seems that only two polymorphic loci are associated with the mating trait and few, if any, other loci appear to be subject to divergent selection. Poelstra et al. (2014) found that only a small number of narrow genomic islands related to the plumage phenotype exhibited resistance to gene flow. The remainder of the genome appears not to experience any barriers to gene flow (see also Vijay et al. 2016). More specifically, Knief et al. (2019) demonstrated divergent selection in the European hybrid zone only for a ∼2 MB region on chromosome 18 and the NDP gene (locus 1 and 2, this study). Geographic clines were narrow for these loci (see Fig. 2), but stretched across several hundred kilometers for the remaining loci in the genome. Overall, these results are consistent with the findings of this study and suggest a plausible history of genome-wide admixture with the exception of few mating-trait loci maintained by assortative mating following a simple matching rule. Divergence in putative matingtrait (and preference) loci against a homogeneous genome-wide background is not restricted to the crow system (Malinsky et al. 2015;Toews et al. 2016;Hench et al. 2019;Wang et al. 2020). The framework presented here thus likely applies more widely, and makes it worthwhile exploring in other organisms tailored to the specifics of each system.
An important element of speciation is reinforcement of assortative mating (Yeh et al. 2018). If hybrids have reduced fitness, even if only in the sense of frequency-dependent sexual selection as explored here, there may be selection in favor of mating strategies that avoid producing hybrids (Felsenstein 1981;Servedio and Kirkpatrick 1997). Thus, reinforcement of assortative mating might be an important factor for speciation (Barton 2013;Butlin and Smadja 2018). It would be interesting to include reinforcement in our theoretical model to analyze how reinforcement could affect the dynamics of the hybrid zone and whether it could lead to speciation before one of the two color morphs becomes fixed. The extent of reinforcement of assortative mating as an evolutionary process, however, is contingent on the amount of variation in the strength of mating preferences in the crow populations and on the genomic architecture of mating-preference behavior. Both are unknown for crows, such that adding reinforcement to our theoretical model would be purely speculative or require extensive behavioral observations and population genomic analyses. As our model neglects the possibility of reinforcement and other eco-evolutionary aspects of the hybrid zone, we cannot rule out the possibility that speciation in carrion and hooded crows takes place. With our analyses, however, we show that the available data can also be explained by a model that will not lead to speciation.

AUTHOR CONTRIBUTIONS
DM and JW conceived of the study. All four authors contributed substantially to the design of the basic rationale of the model, its assumptions, the interpretation of the results and writing of the text.
Further individual author contributions were as follows: DM developed the mathematical model, implemented the simulation software, carried out the model-based data analyses and drafted large parts of the manuscript, especially in the sections on Methods and Results. UK inferred the genotype-phenotype maps. JP designed the figures shown in the main text. JW drafted large parts of the manuscript, especially in the Introduction and Discussion sections.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Figure S1: Overview of how different sources of data (ellipses) are used to fit the components of our model (boxes). Figure S2: Phenotype estimates used in the simulation for an epi static (left) and additive model (right). Figure S3: Observed genotype frequencies in chicks sampled along the Central European transect. Figure S4: Fitting the Gaussian-mixture dispersal model to data of Siefke (1994). Figure S5: Mating preference w g 1 ,g 2 as a function of phenotypic difference ψ(g1)-ψ(g2), with parameter values (ζ , η) = (2.17, 0.617). Figure S6: Genotype frequencies in each bin after 2000 generations in the main model. Figure S7: Genotype frequencies in the neutral model after 2000 generations. Figure S8: Results from neutral model with (ζ , η) = (300, 0.05) and an initial locus 2 dark allele frequency of 0.6 Figure S9: Results from neutral model with (ζ , η) = (300, 0.001) and an initial locus 2 dark allele frequency of 0.6 Figure S10: Results from model with an additive genetic architecture of the coat color trait Figure S11: Results from model with an additive genetic architecture of the coat color trait Figure S12: Genotype frequencies around the hybrid zone after 100 and after 2000 generations after initial contact. Figure S13: Allele frequencies at locus 1 (red) and locus 2 (blue) around the hybrid zone after 100 and after 2000 generations after initial contact. Figure S14: Fitness values of single-locus genotypes at locus 1 (red) and locus 2 (blue) around the hybrid 2000 generations after initial contact Figure S15: According to the categorical model, the hybrid zone first shifts westwards until the hooded crows go extinct Figure S16: Left column: only the third locus (dashed bold line) is neutral, Right column: all three loci are neutral Figure S17: Phenotype distributions in a finite-population self-reference model simulation run after 1, 2000 and 5000 generations after intial contact Figure S18: Clines of the alleles for darker colour for the locus 1 (red) and locus 2 (blue) after 1, 2000 and 5000 generations after initial contact according to a simulation run with the finite-population self-reference model. Figure S19: Clines of the alleles for darker colour for the locus 1 (red) and locus 2 (blue) after 1, 1000 and 2000 generations after initial contact according to a simulation run with the finite-population with female-choice and the imprinting of females according to the average phenotypes of their parents. Figure S20: Clines of the alleles for darker colour for the locus 1 (red) and locus 2 (blue) after 1, 2000 and 5000 generations after initial contact according to a simulation run with the finite-population with female-choice and the imprinting of females according to the phenotypes of their father Figure S21: Clines of the alleles for darker colour for the locus 1 (red) and locus 2 (blue) after 2000 (left column) and 5000 generations (right column) after initial contact according to four independent simulation runs with the finitepopulation with female-choice and imprinting. Figure S22: To get an impression of how well the empirical data fits to our best-fitting model, we compared the log-likelihood of our optimal parameters calculated for the genetic data to log-likelihood values of the same parameter values calculated for data that were simulated according to our model with exactly these parameter values.