As expected, our simulations show that both He and Fst values differ for Y-linked and mtDNA markers when migration patterns differ in males and females. However, our results also demonstrate complex patterns that would have been difficult to predict without simulations. We will concentrate on the patterns observed for the farmers, since they correspond to the modern populations.
(a) General results across all scenarios
Looking at all samples obtained in a particular generation, we found that as time goes from generation 1300–1600, the set of He values becomes more concentrated in one region (figure 2a–c; see also the electronic supplementary material, figure S1). In other words, we found fewer differences on levels of genetic diversity across modern populations (generation 1600) compared with ancient populations (generation 1300).
At the genetic differentiation level, the Fst values against the starting deme (9_9) increase with distance from it, as expected. This increase can be significantly greater for the NRY compared with mtDNA data (as in figure 2f), or the opposite (as in figure 2e) depending on the scenarios (see below), but the Fst values increase with distance from the starting deme. As generations go, Fst values between the starting and last deme decrease with time (figures 2d–f), in agreement with the fact that He values are increasingly similar among samples. Thus, in these simulations, modern populations are genetically less differentiated than ancient ones.
(b) No admixture scenarios
In the scenarios without admixture (figure 2), the He values decrease along the axis of the expansion, with the highest values being observed in the starting deme (9_9) and the lowest in the last colonized deme (0_0). This is observed for all generations sampled. Moreover, these points typically move, as a group and across generations, from higher to lower NRY diversity and from lower to higher mtDNA diversity. In other words, present-day populations have more mtDNA diversity and slightly less NRY diversity than ancient populations, whichever the PMR pattern. Note that this is true when the set of samples from the diagonal are analysed as a group but not necessarily for each sample individually, owing to the fact that the points are also more compact, as we noted above. For instance, the starting deme loses NRY diversity (figure 2), whereas the last colonized deme actually sees its NRY diversity increase. Another very striking result was that the three scenarios (bi-, matri- and patrilocality) exhibit clearly differentiated patterns (figure 2). This can be seen in the way the points are arranged in ‘parallel lines’ through time.
In bilocality scenarios, we predicted similar mtDNA and NRY data. In fact, the points are arranged in a direction parallel to the solid line corresponding to equal values for the x- and y-axes (i.e. for mtDNA and NRY data). Interestingly, we observed that the bilocality He values were higher for mtDNA than for NRY data (figure 2a), whereas Fst values were quite similar (figure 2d). The difference in He values is probably owing to differences in the mutation rates, higher in mtDNA (see the electronic supplementary material, appendix B). Indeed, when we repeated these simulations by assuming the same mutation rate for the two markers, we found symmetrical results (not shown).
In the matrilocal scenarios, all demes had similar NRY He, hence generating values forming ‘lines’ parallel to the y-axis (figure 2b). As expected, the mtDNA Fst values were higher than the NRY Fst values (figure 2e), that were themselves very similar between demes (i.e. the gene flow between demes was high), generating ‘lines’ near-parallel to the y-axis. On the contrary, in scenarios with farmer patrilocality (and still no admixture), a similar behaviour is seen but inverted for the two markers (i.e. higher NRY Fst), and with almost no variation along the y-axis (figure 2f). Similarly, the He values also show this behaviour (figure 2c). A particularly interesting result was that this trend was parallel to the regression obtained from the real data from modern populations (dashed line in figure 2). This was true for both Fst and He values and was not observed in the other scenarios (matrilocality and bilocality).
(c) Influence of hunter–gatherer post-marital behaviour on the farmers genetic diversity (He)
In the scenarios with admixture between HG and farmers, some significant changes are found on the level of genetic diversity (figure 3; see also the electronic supplementary material, figures S1 and S2). First, compared with the no admixture scenarios, the sets of points are shifted towards lower NRY diversity when γ increases, whereas mtDNA diversity does not change very much or shows a slight increase. Thus, in our simulations, admixture leads to a decrease in NRY diversity in all cases compared with the no admixture scenarios. Second, scenarios with patrilocality in HG populations generated fewer changes relative to the no admixture scenarios in NRY diversity compared with bilocality and matrilocality. This is true whether the farmers were patrilocal, matrilocal or bilocal and can be seen in figure 3 where the points with a P (HG patrilocality) are closer to the points of the no admixture scenarios (figure 3a,d,g) compared with the points with a B (HG bilocality) and even more with an M (HG matrilocality). In other words, in situations where HG could mate with farmers, the post-marital behaviour of HG populations clearly leads to differences in the distribution of farmers genetic diversity, in modern populations that have the same PMR system. Third, when farmers are patrilocal, a higher admixture rate (figure 3i) would tend to blur this effect and make the He pattern almost indistinguishable, whichever post-marital behaviour the HG may have had. However, the simulations that seem to better fit the trend of the observed data are the ones from patrilocality in farmers, whatever the HG's post-marital behaviour is and whatever the admixture rate is.
(d) Influence of hunter-gatherer post-marital behaviour on the farmers genetic differentiation (Fst)
Interestingly, in a 10 × 10 lattice, the modern samples Fst values seem less affected by the PMR system of the HG, than the corresponding He values (electronic supplementary material, figure S3). In particular, the analyses of only the last generation data show that the Fst values were nearly identical across all HG scenarios with and without admixture. Conversely, in the generations that follow the admixture events, there were clear differences between the no admixture and admixture scenarios (electronic supplementary material, figure S4). However, in the scenarios analysed using a larger lattice (30 × 30, i.e. a larger geographical area) it was possible to separate the PMR system of the HG, on the basis of modern Fst values (electronic supplementary material, figure S5).
Altogether, our simulations allowed us to study the effect of (i) variable migration rates in males and females within the HG and farmers layers and (ii) variable admixture between layers, on the patterns of genetic diversity and differentiation in present-day and ancient populations.
(a) Main results: (i) farmers were patrilocal and (ii) different PMR systems have a different impact on human genetic patterns
Patrilocality was the most probable scenario among farmers. It was particularly obvious in the no admixture scenarios, but it was also found in the scenarios with admixture, even though not so obvious in some scenarios. This result agrees with ancient DNA (aDNA) and strontium isotope analyses that suggested patrilocality in Linear  and Corded  Ware Culture burials from Germany. Cultural phylogenetics studies also suggest that patrilocality started to increase after the advent of agriculture [29,39].
Changes in PMR systems were also found to lead to different genetic patterns in present-day populations. In particular, the farmers' He values changed significantly depending on whether the HGs were patrilocal, bilocal or matrilocal. Thus, it appears that even though HG populations disappear as far back as 5000 years before the present in our simulations, they influence present-day patterns in modern-day populations.
(b) Behaviour of summary statistics
Pairwise Fst statistics were much less influenced than He by the PMR pattern of the HG populations. While Fst values were different across scenarios after the start of admixture, this signal disappeared in the modern samples.
However, when a larger lattice (30 × 30) is analysed (corresponding to a larger geographical area), it was possible to distinguish between the HG populations PMR systems. This is compatible with the notion that the degree of genetic differentiation between two demes depends on their geographical distance, on the migration rate between local demes and on the time since the populations started expanding. If migration rates are large and/or enough time has passed, then it may be necessary to use large lattices to avoid this homogenizing effect in Fst values. This Fst statistics' dependence on geographical distances implies that inferences based on local/regional sampling is valid only for the most recent history, while sampling from more distant places may be able to recover older patterns, a point that has been stressed by Wilkins & Marlowe .
Furthermore, the access to the genetic composition of ancient HG populations may be not only useful, but necessary to provide us with significant information on this issue. In other words, aDNA may be required to allow us to determine the post-marital behaviour of European HG populations, before and after the Neolithic transition. Currently, the number of Neolithic transition aDNA studies are slowly increasing [38,40–43] but are unfortunately limited to mtDNA. Our results suggest that obtaining NRY DNA from the same samples would be particularly important, as was done in a recent study .
To identify the most probable scenario, we focused on the trend observed in the statistics of both simulated and real data. Our approach was thus to some extent qualitative. To obtain a better fit, one would also need to consider the spread of populations in both the simulated and real data (not just the regression slope). In theory, simulating different scenarios should allow us to better tune migration rates, and identify the original level of diversity in both HG and farmers populations, that are compatible with the observed modern-day data.
(c) Mutation rates can generate asymmetries between NRY and mtDNA data
Although mtDNA and NRY data are often presented as symmetrical counterparts of the female and male demography, respectively, it is not necessarily that simple. The difference in mutation rates can generate an asymmetry between mtDNA and NRY He values in bilocal scenarios with no sex-related variance in reproductive success (figures 2a and 3a). This is something to keep in mind when analysing differences observed in real mtDNA and NRY data because such differences are often interpreted deterministically in terms of differences in male and female behaviours [28,45,46].
(d) Admixture decreases farmers NRY genetic diversity
We found this surprising at first, as it is usually assumed that regions where populations admix will exhibit higher levels of genetic diversity. However, the underlying assumption is that admixing populations have similar Ne. Several studies have shown that during spatial expansions the expanding population is diluted [10,47,48]. We thus believe that as admixture took place between populations with different sizes (i.e. HG having much smaller populations than the farmers), the incoming population will ‘dilute’ the farmers genetic diversity and lead to the decrease in genetic diversity. The decrease is observed even when both HG and farmers had the same pattern of PMR rules. However, this is not necessarily a general result as mtDNA genetic diversity did not always decrease. Again, the difference in mutation rates between mtDNA and NRY markers may interact in a complex way with demographic parameters leading to asymmetries in present-day data.
(e) Comparison with other sex-biased migration studies
Until now, the inference of patterns of sex-biased migration have relied mainly on the comparison and estimation of dispersal from pairwise NRY and mtDNA Fst values [45,46], and by cultural phylogenetics [29,39] in modern populations.
Hamilton et al.  also used a spatial framework, using NRY and mtDNA data, to study matrilocal and patrilocal groups from northern Thailand. In their study, the authors applied a modified version of SPLATCHE and were not interested in detecting shifts in PMR patterns, which were assumed to be invariant in their simulations. Instead, their aim was to compare male and female migration rates in known patrilocal and matrilocal societies that would explain present-day levels of genetic differentiation and diversity. Here, our aim was to understand how PMR (with or without shifts) interacts with admixture between different societies, to generate differences in maternally or paternally inherited markers analysed jointly.
Model-based approaches have many advantages as they allow us to identify parameters that have a significant impact on the data. However, they also rely on strong assumptions. In our study, it was necessary to make assumptions on the level and patterns of gene flow, carrying capacities and genetic make-up of the founder populations, which suggests that some of the conclusions presented here should not be taken at face value. Hence, we believe that the general trends identified are to some extent robust. For instance, our simulations were performed on a 10 × 10 lattice. But when we repeated the patrilocal scenarios on a 30 × 30 lattice, we found essentially the same results, the main difference being that the power to identify scenarios was increased in the 30 × 30 lattice.
The simulated framework introduced here owes much to the work of Currat et al. , but is sufficiently different to represent an interesting alternative to identify the critical assumptions that are robust and those that are not, and the type of data required to separate scenarios. Altogether, our simulations helped identify important parameters and scenarios, together with data that would be needed to study the Neolithic transition in Europe (NRY aDNA), but much work is still necessary.
There are still very few studies that have dealt with the kind of complex scenarios that involve the characterization of the expansion of two demographically different populations across the same geographical area when migration patterns and admixture levels vary, and those that exist do not deal with sex-biased migration [17,49]. Our work provides some of the first insights into the consequences of complex demographic changes that probably took place during the European Neolithic, on present-day human genetic patterns.