In this study, we capitalize on the fact that genetic diversity within and among demes (or population subdivisions) is highly dependent on local levels of gene flow after a spatial expansion (15, 16) to get direct estimates of sex-specific migration rates in patrilocal and matrilocal populations. We apply a recently developed Bayesian framework (17) to the problem of estimating the number of male and female immigrants arriving each generation into the six traditional Thai populations mentioned above. These populations are particularly interesting because they present contrasting postmarital residence patterns while sharing similar habitats and living in the same geographic region. They thus provide a unique opportunity to infer gender-specific dispersal rates and to contrast them to social patterns in human populations.

Samples. The population samples consist of three patrilocal populations (Ahka and two groups of Lisu that are 220 km apart) and three matrilocal populations (Lahu, Red Karen, and White Karen) from northern Thailand that were previously analyzed for 300 bp of the mitochondrial control region and for nine Y chromosome short tandem repeats (STRs), as described in ref. 14. Population sample sizes (mtDNA/Y chromosome STRs) were: Akha (91/21), C Lisu (53/9), M Lisu (42/22), Lahu (39/17), Red Karen (39/30), and White Karen (40/20).

Bayesian Estimation Procedure. The recently developed approximate Bayesian computation (ABC) method (described formally in ref. 17) allows the estimation of the parameters of complex evolutionary models for which explicit likelihood cannot be formulated analytically or for which computation is prohibitive (18–21). It consists of simulating a large number of data sets (typically several hundreds of thousands or more) under a given evolutionary model whose characteristics (e.g., number of samples, individuals per sample, and number of loci) match those of some observations. The parameters of the simulated model are drawn from specific prior distributions, and simulated data sets are compared with the observations by means of a Euclidean distance computed from summary statistics capturing various aspects of the data (17). The simulated data sets closest to the observed data and their associated parameters are retained for subsequent parameter estimation procedure; the other data sets are discarded. Finally, a multiple regression of summary statistics on parameters is performed to obtain posterior densities, from which estimators are inferred (17).

In the present study, the evolutionary scenario that is simulated is a range expansion in a two-dimensional stepping stone (2DSS) model. As previously shown (15, 16), it depends on just three parameters: the local deme size N, the age of the expansion T, and the migration rate m between adjacent demes (see below). Note that the current methodology is drastically different from previous approaches, which attempted to estimate migration between sampled populations (22, 23), because, instead, we estimate the average number of immigrants in sampled demes (15, 16) without needing to sample the population of origin of these migrants.

Simulation of a Range Expansion. For our evolutionary scenario, we employ a modified version of the 2DSS range expansion model described in refs. 15 and 24. We assumed that observed samples were drawn from a large subdivided population, simulated on a torus consisting of 50 × 50 subdivisions (demes). We assume that the current subdivided population was restricted to a single (ancestral) deme of size N until T generations ago, at which time a range expansion occurred that led to the instantaneous colonization of all of the 2,500 demes, each one having the same size N. After the range expansion, the neighboring demes are assumed to exchange genes at a rate m, so that a total of M = Nm genes are immigrating into each deme per generation.

The consequence of the range expansion is thus a large increase in the total size of the subdivided population through an increase of the number of colonized demes, but the size of each deme does not change compared with the size of the ancestral deme. The instantaneous range expansion used in this model leads to simulation times that are vastly reduced in comparison with a progressive expansion model used earlier (15), without a major impact on inferred parameters (see Table 2, which is published as supporting information on the PNAS web site). This speed gain is important to allow for the very large number of simulations in this study, a key factor in accurately estimating parameters under an ABC approach. Another feature of the instantaneous expansion is that it makes it unnecessary to specify the geographical origin of the expansion relative to the location of the present samples.

Genetic diversity of samples of DNA sequences and Y chromosome STR markers were then simulated by using a standard backward coalescent process after each range expansion for three samples located randomly on the grid, but not directly exchanging migrants. Mutations were added at a rate μ on the simulated coalescent tree. Simulations were thus based on values of the four parameters N, T, m, and μ, which were randomly and independently drawn from prior distributions for each simulation. For mtDNA and Y chromosome STR data, the following priors have been used: N, log-uniform between 100 and 10,000; T, uniform between 1,000 and 3,333 generations, which represent 30–100 thousand years (ky) when taking a generation time of 30 years, which seems appropriate for humans (25); and m, log-uniform between 0.0005 and 0.5. Different priors were taken for mtDNA sequence and Y chromosome STR mutation rates. The prior for the total mutation rate of the 300 bp of the mtDNA control region was taken as uniform between 0.005 and 0.002, corresponding to ≈7–28% divergence per million years, which is in the higher range of previous phylogenetic estimates (2.5–26%) (26). For Y chromosomes, the prior for the total mutation rate of the nine STRs was taken as uniform between 0.002 and 0.02, corresponding to 2.2–22 × 10^–4 per generation per STR encompassing previous estimates (e.g., ref. 27). After every cycle in which priors were drawn, a demographic simulation was performed, and genetic diversity was simulated, summary statistics were computed on the simulated data and recorded into a file with the associated values of the parameters used for the simulation.

Summary Statistics. Summary statistics were chosen based on a literature search and previous investigations. For mtDNA data, four within-deme summary statistics were calculated for each sampled deme: number of haplotypes, k; homozygosity, H_o; number of segregating sites, S; and the average number of pairwise differences, π. For Y chromosome STRs, three withindeme summary statistics were calculated for each sampled deme: average number of alleles per locus, a; homozygosity, H_o; and mean variance (across loci) in allele repeat number. The fixation index, F_ST, was calculated as Weir and Cockerham's θ_W (28) among sampled demes for both marker types. For each of the four sets of three population samples sharing a residence pattern (matrilocal mtDNA, matrilocal Y chromosome STRs, patrilocal mtDNA, and patrilocal Y chromosome STRs), we performed 5 million simulations of the range expansion to ensure the greatest possible accuracy of parameter estimates. For each of the four data sets, the 2,000 simulations leading to summary statistics closest to the observed statistics were retained for estimation of the parameters of interest here, which are N, T, and M = Nm. Following exactly the approach of Beaumont et al. (17), we assessed the similarity between observed and simulated summary statistics by means of a Euclidean distance, and we estimated the posterior densities of the parameter values by a multivariate weighted local regression procedure. The only step differing from the ABC estimation procedure described in ref. 17 was the application of the following transformation of the parameter values before the regression:

In agreement with previous work (14), we find that males move more in matrilocal populations than in patrilocal populations and that the reverse is true for females (Table 1). Compared with the previous approach, we can accurately quantify the number of new immigrants of each sex entering the different populations. Our results suggest that, on average, less than one male enters patrilocal populations every generation, whereas approximately four males are incorporated each generation in matrilocal populations, implying that males move ≈8 times more in matrilocal societies than in patrilocal societies (3.9/0.5). Even though this difference appears important, it is not strictly significant because the 95% credible intervals of the posterior distributions overlap. However, the fact that the credible intervals are only slightly overlapping suggests that this difference is real (Table 1 and Fig. 1). Females seem to move only ≈2 times more in patrilocal societies than in matrilocal societies, but the posterior distributions of these quantities here are quite broadly overlapping (Fig. 1).

Table 1. Estimation of the effective deme size N, the age of the range expansion T, and the number M = Nm of male and female immigrants entering populations with different postmarital residence patterns

Fig. 1. Posterior densities of the number of immigrants M received each generation by a single subdivision. Densities estimated in patrilocal populations are represented by solid lines; those of matrilocal populations are shown with dashed lines. Densities estimated (more ...)

Another innovative feature of this study is the possibility to compare male and female movements within and across social organizations. We find that hill tribe patrilocal societies are characterized by an extreme philopatry of males (less than one incoming migrant per generation) and a high number of incoming females (more than seven per generation), resulting in ≈15 times (7.4/0.5) more women than men being exchanged between neighboring populations. This difference is found to be highly significant because the credible intervals of these quantities are nonoverlapping (Table 1), and it is 2 times more than a previous genetic estimate obtained by comparing mtDNA and Y chromosome geographical patterns in Europe (6), where populations are also predominantly patrilocal (29). In contrast, the ratio of female versus male immigrants is much more balanced (0.79) in matrilocal hill tribe societies, where the number of incoming males (3.9) is very similar to that of females (3.1). This result suggests that patrilocality and matrilocality are not just the opposite sides of the same coin but that different regulatory mechanisms are at work in these social systems. The data are consistent with a more important role of men in matrilocal societies than that of women in patrilocal societies concerning the choice of postmarital residence, which would result in a strict limitation on the number of incoming men in patrilocal hill tribes and a much looser control on the immigration of women in matrilocal societies. Although ethnographic data are not extensive for these hill tribes from northern Thailand, available demographic data suggest that the Ahka and Lisu tribes are strictly patrilocal (30, 31), whereas the Lahu and Karen groups are less rigidly matrilocal (30, 32), with a proportion of marriages leading to residence in the village of the male. This pattern would imply the existence of some level of female gene flow among villages that are predominantly matrilocal, which is confirmed here by our results (Table 1).

One significant strength of the ABC procedure is that it can be used in association with models of any complexity given only that data can be simulated under the model (17). Coupling this estimation procedure with the spatially explicit 2DSS model described above thus has the potential to provide greater accuracy of migration and range expansion estimates than methods that do not account for a spatial component (e.g., ref. 6). A modified ABC method was recently used with this 2DSS model and showed accurate estimates of the composite parameters M = Nm and τ = 2μT across a wide range of parameter combinations (33). As shown in Table 2, the estimates of the expansion time T, the composite parameter τ = 2μT, and the number of immigrants M in an instantaneous range expansion are very well recovered, whereas the effective deme size N appears overestimated for n < 500. We also simulated a progressive range expansion (as in ref. 15) and checked whether the parameters of this type of expansion were also correctly estimated under our ABC procedure, assuming that an instantaneous range expansion occurred. As reported in Table 2, the accuracy of the estimations of the expansion time T and the effective deme sizes N are generally good and similar to those inferred under an instantaneous expansion, whereas the number of immigrants M is only slightly underestimated. This underestimation is probably due to an increased rate of coalescence expected during the progressive expansion, because demes at the front of the expansion wave have not reached their carrying capacity (34, 35). We also checked the impact of the size of the grid we used for our simulations by estimating the parameters of a range expansion for the females in matrilocal populations in a torus of 30 × 30 demes instead of a torus of 50 × 50 demes. The result reported in the last row of Table 1 shows that the grid size affects neither our point estimates nor our credible intervals. This result is in keeping with previous results obtained in the finite-island model (16), where it was shown that the genetic diversity within demes was not much affected by the exact number of demes, provided that this number was large (a few hundred).

As expected, the estimates of the timing of the spatial expansion are found to be very similar across genders and social organizations (Table 1) and point to a late Pleistocene expansion ≈1,200–2,000 generations ago (40–60 ky ago, with a generation time of 30 years; ref. 25), in keeping with a previous estimate of expansion signal for non-African populations (1,925 generations; ref. 36). Note that our estimates take into account the uncertainty in the mutation rate of the different markers, because they are inferred directly from the posterior distribution of the expansion times. These estimates are also in relatively good agreement with previous reports of Pleistocene population expansions obtained from mtDNA (>50 ky; refs. 13 and 37) but are slightly larger than those obtained from Y chromosomes (20–35 ky; refs. 18 and 27) or nuclear markers (35 ky; ref. 38).

The estimate of effective deme sizes for males and females in different societies are also reported in Table 1. As expected, effective deme size is found to be much smaller for males in patrilocal societies than in matrilocal societies, and similarly small deme sizes are inferred for females in matrilocal societies. There is thus a clear correlation between effective deme size and dispersal rates, suggesting that our estimates of deme sizes simulated under a discrete stepping-stone model could correspond to neighborhood sizes (39) in a continuous population (see ref. 40 for further discussion on the connection between the coalescents in continuous and 2DSS models). In other words, the different dispersal abilities of males and females imply that the neighborhood size could be a single village for, say, males in patrilocal societies, whereas it could encompass several villages for females in the same ethnic group, explaining the larger inferred deme size for females.

Overall, our study confirms the gender-related asymmetry in levels of gene flow between patrilocal populations (6, 14, 41) and shows the differential strengths of the social regulation of dispersal patterns among patrilocal and matrilocal populations. Given the diversity of social and mating systems among human populations (42) and the many factors controlling postmarital residence choice in tribal groups (43, 44), we would not generalize our results to other ethnic groups and other regions of the world. Moreover, it is still unclear how local patterns of sex-bias dispersal may affect genetic structure at a larger continental scale (11). However, we show here that accurate information on social organization and sex-specific dispersal can be retrieved from genetic data alone. It should prove extremely valuable in the human context, because proper assessment of effective dispersal rates and residence patterns are difficult to obtain from point observations, given the potential changes in residence location throughout man's reproductive life (45). Similar approaches should be applicable to other social species for which long-term and direct observations are difficult to perform.

We thank Mark Beaumont for his comments on a previous version of the manuscript and for his advice. We are also grateful to two anonymous reviewers for their comments and several useful suggestions. G.H. and L.E. were supported by Swiss National Science Foundation Grant 3100A0-100800.