Spotted phenotypes in horses lost attractiveness in the Middle Ages

Horses have been valued for their diversity of coat colour since prehistoric times; this is especially the case since their domestication in the Caspian steppe in ~3,500 BC. Although we can assume that human preferences were not constant, we have only anecdotal information about how domestic horses were influenced by humans. Our results from genotype analyses show a significant increase in spotted coats in early domestic horses (Copper Age to Iron Age). In contrast, medieval horses carried significantly fewer alleles for these phenotypes, whereas solid phenotypes (i.e., chestnut) became dominant. This shift may have been supported because of (i) pleiotropic disadvantages, (ii) a reduced need to separate domestic horses from their wild counterparts, (iii) a lower religious prestige, or (iv) novel developments in weaponry. These scenarios may have acted alone or in combination. However, the dominance of chestnut is a remarkable feature of the medieval horse population.


| | |
This method showed a notable capacity to deal with the complexities of our inferential problem. In contrast, other available methods can't deal with complicated relationships between genotypes and phenotypes, as they assign selection coefficients directly to the genotypes. From the many methods available to estimate selection from time series of genetic data, only two 3,4 are able to deal with dominance/recessive relationships between alleles but none is able to deal with epistatic relationships among genes. It looks also unlikely that the methods described in [1][2][3][4][5] can deal with uncertainty in samples ages which is relevant because it makes the order of the samples uncertain and with them the observed frequencies in time. Another option for inferring selection from time series of allele frequencies consist in using approximate Bayesian computation. It has the same (if not a larger) potential to deal with all the interactions and uncertainties, but it is computationally more demanding and probed incapable for the inference of our system due to the number of parameters.

Noise variables
The simulated part of our model can incorporate necessary parameters for which there is neither available knowledge nor sufficient power in the data to estimate them. So the procedures include them as well as their priors and marginalize the target posterior by integrating over their ranges. The parameters that we treated in this way were the effective population sizes and the generation time: The simulations sample their values from their prior distributions before inserting them in the simulations. It is important to notice that for each new noise variable, the number of simulations should be increased in the same way as if the number of parameters were increased.

Algorithm
We implemented this algorithm in our problem of estimating selection coefficients for nine periods in a sample of 201 horses spread between the Late Pleistocene and medieval times. The algorithm has the next steps: 1. Define the proper parameters, priors and noise variables. Also define, a model and create a program for the simulation of paths of allele frequencies.
2. Sample a set of parameters from their prior distribution .
3. Simulate an allele frequency path, , using Wright-Fisher explicit simulations with the desired level of complexity and employing the values of obtained in the previous step. 4. Calculate analytically the partial likelihood of the i-th simulation | .
5. Use an appropriate kernel for proposing a new set . 6. Simulate an allele frequency path, , using Wright-Fisher explicit simulations with the desired level of complexity and employing the values of obtained in the previous step. 7. Calculate analytically the partial likelihood | .
8. Accept or reject , with probability • (or with a Metropolis-Hasting or another criterion). 9. Go to 2.
10. Repeat 1-9 for a sufficiently large number of steps.
11. Run several chains and evaluate different transition kernels. Evaluate and optimize as a regular MCMC procedure.
The algorithm is also able to incorporate some of the many refinements developed under the theory of MCMC. One of them, a Gibbs sampler could be particularly useful for highly multivariate problems. In our algorithm the Gibbs sampler updated the chains one parameter at a time.

Initial states and introduction of alleles
The different methods that have been designed for estimating selection by means of time series of allele frequencies have dealt with the initial state of a derived allele (the one under selection) in two alternate ways: (i) the derived allele exist since the beginning of the timeframe of interest, requiring the estimation of its initial allele frequency; or (ii) considering that the mutant allele must have appeared by mutation at some time, the estimation of the initial allele frequency is substituted by the estimation of age of the allele.
Considering that our simulations are, in general, more computationally demanding than other approaches (e.g. analytic solutions), our method will generally be circumscribed to the sampling timeframe in order to minimize computational demands. For that reason, the default option will be (i). However, if the mutant allele is absent from the first sample there could be a non-null probability that the allele appeared by mutation inside the timeframe of the simulations. Opting for (i) or (ii) could be solved by performing a Bayesian model comparison between those models. However, it is possible to implement a MCMC that can jump between models in the same way that it jumps between parametric states. This has been done in phylogenetic inference for choosing models of molecular evolution on-the-go 6 .
In our case, the oldest sampling of the derived allele was always much younger than the oldest sample, requiring this hybrid approach.

Programing
We created three versions of the algorithm: one for one-gene-systems, one for basic colours (programing the interactions between the genes ASIP and MC1R) and one with all the eight genes and their interactions.
Our program consisted in:

Bottom line
Our method has several innovative steps but in its general form is analogue to the estimation of demographic parameters by coalescent approaches 9 . It has been successfully tested in a simpler system, and showed more accurate than approximate Bayesian computation while being significantly more efficient (see Figure S7 & 10 ). Such features could make this technique a useful complement of available methods not only for inferring selection but in other types of inference, especially in those in which complexity precludes the use of other approaches. Figure S1: Allele frequencies over time of derived, non-wild-type alleles that are associated with a specific coat color phenotype. Figure S2. Grouping of samples for the temporal test of allele frequencies. Only the samples inside boxes were included in the test. Notice that the temporal scale changes and is discontinuous for the Pleistocene. . The time frame ranges from 12.5 ky BP to present (from left to right); most of Pleistocene was omitted for visualization purposes. Each line corresponds to one simulation and the color black corresponds to the initial simulation, the red to simulations of the burnin stage of the MCMC, and blue to simulations employed for the inference. Figure S4. Relationships among phenotypes and genotypes associated to two dilution genes and two basic colors genes. The arrows mean that a dilution gene can only be expressed in the specific color the arrow are directed to (notice that, for instance, cream heterozygotes, C/cr, cannot be expressed in black horses, ASIP a/a). The t-shape ending of the lines mean that the indicated phenotypes are dominant-repressed by the genotype they come from. The genes that doesn't appear in the figure, KIT13-Tobiano, KIT16-Sabino, MATP-Pearl and TRPM1-LP present phenotypes equivalent regardless the basic color and only present relationships of dominance between alleles of the same gene.