## Abstract

The number of genes controlling mimetic traits has been a topic of much research and discussion. In this paper, we examine a mimetic, dendrobatid frog *Ranitomeya imitator*, which harbours extensive phenotypic variation with multiple mimetic morphs, not unlike the celebrated *Heliconius* system. However, the genetic basis for this polymorphism is unknown, and not easy to determine using standard experimental approaches, for this hard-to-breed species. To circumvent this problem, we first develop a new protocol for automatic quantification of complex colour pattern phenotypes from images. Using this method, which has the potential to be applied in many other systems, we define a phenotype associated with differences in colour pattern between different mimetic morphs. We then proceed to develop a maximum-likelihood method for estimating the number of genes affecting a quantitative trait segregating in a hybrid zone. This method takes advantage of estimates of admixture proportions obtained using genetic data, such as microsatellite markers, and is applicable to any other system where a phenotype has been quantified in an admixture/introgression zone. We evaluate the method using extensive simulations and apply it to the *R. imitator* system. We show that probably one or two, or at most three genes, control the mimetic phenotype segregating in a *R. imitator* hybrid zone identified using image analyses.

## 1. Introduction

The analysis of phenotypic and genetic variation in geographical areas where two or more phenotypically distinguishable groups of organisms meet and exchange genes has been of substantial interest to evolutionary biologists [1–3]. The evolutionary dynamics in these zones, referred to as hybrid zones, introgression zones or admixture zones depending on context, provide us a basis to study processes relating to speciation and for understanding the genetic and ecological underpinnings of adaptive traits, including mimetic and aposematic traits. Substantial work has been done on such systems, including the now classical work on the *Bombina bombina* versus *Bombina variegata* hybrid zone [4,5] and the hybrid zones between various species of *Heliconius* butterflies [6–8]. Of primary interest in these studies is to understand the genetic basis of the phenotypic traits, how selection is affecting these traits, and to understand the relative role of population history, gene flow and natural selection in determining the evolutionary dynamics of the hybrid zone. Furthermore, there has recently been renewed interest in mapping the genetic variants associated with reproductive isolation or adaptive traits in the hybrid zone using so-called admixture mapping or mapping by admixture linkage disequilibrium [9–14]. Of special interest to us are admixture zones, exemplified by the previously mentioned examples in *Bombina* and *Heliconius*, in which complex morphological traits such as colour patterns are segregating and are probably of adaptive significance [15,16].

Mimetic traits in admixture zones, or otherwise, have often been hypothesized to be associated with so-called supergenes [17]. Supergenes are tightly linked clusters of genes that control a suite of traits that will allow Mendelian, or close to, Mendelian behaviour of the mimicry trait. The existence of such supergenes could help us explain the strong phenotypic correlation between many different phenotypes required to produce pure mimetic forms. If many traits are needed to produce an adaptive mimetic phenotype, then one would expect selection to favour genetic variants that increase the correlation between these traits. Much discussion has ensued on the existence of supergenes, particularly in relation to mimetic phenotypes in butterflies [18]. A recent paper by Kunte *et al.* [19] shows that a proposed supergene underlying memetic phenotypes in *Papilio* butterflies in fact is a single Mendelian gene that serves as a genetic switch for the mimetic type. For both *Papilio* and *Heliconius*, it appears that the mimetic phenotypes are often controlled by one or a few genes or supergenes that behave in a largely Mendelian fashion. However, the degree to which mimetic phenotypes have a similar genetic basis in other systems is uncertain.

The dendrobatid frog *Ranitomeya imitator* [20,21] provides us with a new vertebrate model system that shares many features with the well-known *Heliconius* system. In Peru, there are four distinct colour pattern morphs of *R. imitator* that occupy different geographical regions [22]. In each of these regions, the colour pattern of *R. imitator* clearly resembles that of a co-occurring species of dendrobatid frog [20,21]. Phylogenetic analyses indicate that these co-occurring species generally diverged prior to the divergence between the divergent populations of *R. imitator* [23,24]. Evidence for rapid divergence under selection [21,22], and the similarity of each *R. imitator* colour pattern morph to the more anciently diverged co-occurring species, indicates that *R. imitator* has undergone a mimetic radiation, in which different populations have evolved to resemble distinct colour patterns displayed by the local ‘model’ species [21,22]. Recent analyses of colour pattern variation, genetic structure and gene flow have identified multiple zones of admixture where distinct colour pattern morphs of *R. imitator* come into contact and interbreed [20]. These regions vary in terms of the width of the zone of admixture and the degree of genetic divergence (in neutral markers) found across the zone, making this system useful for comparative analyses of divergence. In one region, the zone of admixture is fairly broad (7 km), and populations in the zone show high variability that appears to include elements of both distinct colour pattern morphs. Note in figure 1 that frogs at one end of the admixture zone, where they are mimetic with *Ranitomeya summersi*, tend to be banded with black and orange legs, whereas frogs, on the other end, where they are mimetic with *Ranitomeya variabilis*, tend to be striped with a reticulated green and black pattern on the legs. The genetic basis of this polymorphism is of primary interest, but given the large genome sizes of dendrobatid frogs (e.g. up to 9 Gb, Camper *et al.* [25]), lack of genetic resources and difficulty of captive breeding, direct mapping of the genes involved is a non-trivial task. An objective of this paper is instead to obtain more information about the genetic basis of this polymorphism, solely using image analyses and limited microsatellite typing. In particular, we will be interested in examining whether the polymorphism is controlled by a single Mendelian gene, perhaps a supergene, or by multiple genes.

In order to address this problem, we will first develop automated methods for describing complex colour pattern phenotypes based on images that can be applied in this system and other systems. The advantage of such methods is that they are not subject to the same biases that may occur when a researcher chooses which traits to measure after having observed the images. In addition, such methods may have the potential for identifying important biological features that were otherwise not readily identifiable.

We will then proceed to develop a method for estimating the number of genes affecting a phenotype in an admixture/hybrid zone. For natural populations, in which controlled crosses are difficult or expensive to carry out, and for which parent–offspring pairs cannot easily be sampled, there are no appropriate methods for determining how many genes affect a trait. In other settings, there has been substantial previous work on this problem. The well-known Castle–Wright estimator [26,27] is based on the amount of segregating variation observed in the offspring of controlled crosses of inbred lines. The objective is to estimate the effective number of loci controlling a quantitative trait, i.e. the number of loci required to explain the variance in the trait if all loci have the same effect. There have been numerous extensions of the method, including the incorporation of linkage and variation in effects size [28,29]. Lande [30] showed that the assumption of complete homozygosity in the parental lines is not necessary and provided an estimator applicable to natural populations, rather than to controlled crosses. Building on the idea, dating back to Pearson [31], that the relationship between the variance in the offspring phenotypic values and midparent value depend on the number of genes controlling the trait, Slatkin [32] provided another estimator applicable to outbred populations.

We are interested in estimating the number of genes affecting a trait in a hybrid/admixture zone. This is a problem that has been considered by Szymura & Barton [4] who, based on theory developed in Barton [33] and Barton & Bengtsson [34], estimated the number of genes contributing to selection against gene flow in the *B. bombina* versus *B. variegata* hybrid zone using comparisons of the amount of linkage disequilibrium at the centre of a hybrid zone to the width of the cline. The method we will develop is in the spirit of the of the Castle–Wright–Lande estimators, but is based on using a genetically inferred admixture proportion in each individual. This method does not require data on controlled crosses. It also does also not rely on any assumptions regarding selection models and processes shaping linkage disequilibrium. It is less ambitious in that it does not attempt to determine the number of genes affecting fitness, but the number of genes affecting an observable phenotype. There is substantially less information regarding the number of loci when controlled crosses have not been performed. However, as we will show, there is still sufficient information to distinguish between hypotheses regarding one, two or several genes affecting the trait.

We will apply the methods developed in this paper to images and genetic data from the aforementioned dendrobatid frog *R. imitator*. Using these methods, we can estimate the number of genes affecting the mimetic phenotype without the use of experimental crosses or mapping approaches.

## 2. Image analysis/quantitative phenotyping

A common way to quantify variation in image analysis is to extract a number of so-called descriptors, combine these into a vector of measurements for each individual and use statistical decomposition methods to condense the collected information. Prior to analysis, all individuals have been warped to a mean shape determined by Procrustes analysis [35]. Manual annotation of 22 anatomical landmarks was used to establish point correspondences.

Descriptors are typically designed to capture elementary characteristics of an image, such as colour or shape. Individually, descriptors are usually too specific, but a well-chosen suite of descriptors can provide us with a rich basis for further analysis.

In our study, we use three different phenotypic descriptors: colour/non-colour ratio, gradient orientation histograms and shape index histograms [36], each of which is defined on the pixel-level and described in detail in the electronic supplementary material, S1. These descriptors collect local zeroth-, first- and second-order information about the image. In the current setting, these three standard descriptors can loosely be thought of as measuring features relating to the proportion of coloured area, the degree to which changes in colour occur along the anteroposterior axis or along the left–right axis (banded patterns versus striped patterns), and the degree to which the pattern consists of stripes/bands as opposed to reticulation, respectively. The quantified information is visualized in figure 2 and in the electronic supplementary material, S1.

All descriptors are extracted on a per-pixel basis and pooled together at four distinct interest points, namely each of the frog's legs, lower back (dorsum) and on the back of the head. An interest point is defined in terms of its coordinates and a radius *r*_{k} > 0. Thus, an average of each of the descriptors is accumulated for the four regions shown in figure 2*c* (see also the electronic supplementary material, S1).

### (a) Revealing a mimicry-related phenotype with sparse discriminant analysis

The collected phenotypic descriptors are here condensed into a single mimicry-related phenotype. This amounts to determining the low-dimensional manifold, in the high-dimensional feature space, describing the phenotype. We have chosen to use sparse discriminant analysis (SDA) by Clemmensen *et al.* [37] for this task. More detail on this procedure can be found in the electronic supplementary material, S1, where alternative methods are also considered.

The composite phenotype is constructed as the linear combination *β* of the *p* descriptors that best describes the mimicry across the hybridization transect, i.e. the direction in the *p*-dimensional space that maximizes the ratio of the between-group variance to the within-group variance under elastic net regularization [38].

We define the mimicry-related phenotype for the *i*th individual as the projection onto the one-dimensional subspace spanned by *β*:
2.1where is the *j*th descriptor value for the *i*th individual. For all *n* individuals, this is equivalent to

## 3. A likelihood method for identifying the effective number of genes

We are interested in estimating the effective number of genes, *K*, affecting a trait, i.e. the number of genes required to explain the observed phenotypic variation assuming all genes have the same effect. We assume we have a sample of *n* individuals from an admixture zone, each with some associated genetic data (e.g. microsatellite data). We will take advantage of the fact that even limited genetic data can be used to infer an admixture fraction for each individual, under the assumption that pure forms exist at each end of the transect in the admixture zone. *f* and 1 − *f* then represent the proportion of an individual's genome that is identical to individuals in the right and left end of the transect, respectively. The method we use for estimating the admixture fractions is described in the electronic supplementary material, S1, and is based on the kernel discriminant analysis (KDA) of Mika & Ratsch [39] with the kernel suggested by Martin [40]. The estimated admixture proportions are shown in the electronic supplementary material, S6.

We will assume that the phenotypic values, are normally distributed, given the underlying genotype, and that each locus contributing to the phenotype has the same effect and dominance factors, and that the effects are additive among loci. We will also assume that each locus is di-allelic and that the allele favouring the phenotype in the right end of the transect has frequency 1 in the right extreme of the transect and frequency 0 in the left end of the transect. We will also, without loss of generality, denote the alleles favouring the phenotype in the right and left ends of the transect by *a* and *A*, respectively. An individual with admixture proportion *f*, assuming independence among the parental contributions, then has genotype *AA* in any locus with probability (1 − *f*)^{2}.

We consider the phenotype, *z*, of an individual to be a realization of the stochastic variable *Z* with the conditional distribution
3.1where is the environmental variance and is a vector of the *K* genotypes
3.2Three averages are used for the conditional Gaussians *μ* = [*μ*_{0}, *μ*_{1}, *μ*_{2}]^{T} and
i.e. a vector containing fractions of the *K* genes having the genotypes *AA*, *Aa* and *aa* respectively.

So, for example, if *K* = 3, an individual with genotypes *AA*, *AA* and *aa* in the three loci, respectively, will have mean phenotype 2*μ*_{0} + *μ*_{2}.

Thus, in a noise-free scenario, a single gene would be able to explain a trait as a piecewise constant function (of the admixture proportion) with three steps. *K* genes would be able to explain a trait attaining different values. Here, a noise-free scenario would mean no environmental variance in the phenotype *and* no noise caused by the quantification of the phenotype.

To calculate the likelihood, all possible combinations of genotypes must be considered. The set of all possible combinations will be denoted i.e. the *K*th Cartesian power of possible genotypes, where a single tuple from this set will be denoted This set consists of all possible combinations, with replacement, where the order is significant. A total of 3* ^{K}* such combinations exists.

The probability of a certain combination of genotypes given the mixture proportion *f* is

The likelihood of observing the phenotypic trait over the entire population, allowing *K* genes to contribute to the expression of the trait, is modelled as
3.3

However, the estimates of *f _{i}* may be associated with statistical uncertainty. Ignoring this uncertainty could lead to biased estimates. We therefore provide an alternative formulation that incorporates uncertainty in the estimates of

*f*using a bootstrap approach, i.e. we assume that marker loci used for estimation of

_{i}*f*have been bootstrapped to provide a bootstrap distribution The likelihood of observing the phenotypic trait over the entire population, allowing

_{i}*K*genes to contribute to the expression of the trait, is then modelled as 3.4

For a fixed value of *K*, we maximize this function for *μ*_{0}, *μ*_{1}, *μ*_{2} and using the Broyden-Fletcher-Goldfarb-Shanno algorithm [41]. We then repeat this procedure for multiple values of *K* and choose the value of *K* that maximizes this profile likelihood function as our maximum-likelihood estimate of *K*. To increase the probability of converging to a global maximum, we use a scheme with multiple starting points, see the electronic supplementary material, S2 for details.

We evaluate the performance of the method using simulations allowing for varying heritability and uncertainty in the estimates of *f*. The heritability is defined as the fraction of the total phenotypic variance *V*_{P} that can be attributed to genetic variance:
3.5

The average phenotypic value is where *z _{j}* is the phenotypic value determined by the genotype and

*p*is the proportion of individuals with the

_{j}*j*th genotype.

The genetic variance is determined as 3.6

To simulate data for a phenotype determined by *K* genes, *n* mixture proportions are drawn, e.g. from a uniform distribution on the interval [0, 1]. The genotype for each of the *K* loci is then drawn from a multinomial distribution with probabilities as in equation (3.2). Phenotypes are then assigned by simulating from a normal distribution as in equation (3.1). In simulations with noise in the estimate of *f*, we simulate *B* samples from a normal distribution with standard deviation *σ*_{f} around *f _{i}*, such that mixture proportions used for inference

## 4. Image and microsatellite data

We used published microsatellite data from three sources: Twomey *et al.* [20] (92 samples), Twomey *et al.* [42] (36 samples) and Twomey [43]. The final dataset consisted of 285 *R. imitator* individuals from 16 localities in Peru: the 11 localities shown in figure 1 and five localities between Santa Rosa de Chipaota and Achinamisa (i.e. within the banded-striped transition area). For the unpublished microsatellite data, amplification methods follow Twomey *et al.* [20].

We used JPEG compressed images of six *R. summersi*, seven *R. variabilis* and 304 *R. imitator* individuals from the 11 localities shown in figure 1. The images are 3888 × 2592 pixels of size captured with a Canon EOS Rebel XS SLR. Both microsatellite data and image data were available for 179 of the *R. imitator* individuals. See the electronic supplementary material, S5 for a full overview of the image material.

## 5. Results

### (a) Phenotypic descriptors

The phenotypic descriptors described in §2 were automatically extracted from all 317 images. Different aspects of the patterns in the population are captured by this collection of descriptors, the most dominant being the stripe directionality; for more detail on the phenotypic variance captured by these descriptors, see the electronic supplementary material, S1. For every individual, the suite of descriptors extracted for the four interest points (left leg, right leg, lower back, upper back) are colour/non-colour ratios for each point of interest, gradient orientation histograms binned in two bins at scales *σ* = [2, 7] with tonal range *β* = 1 and shape index histograms in *five* bins at scales *σ* = [4,8] with tonal range *β* = 1. This adds up to a total of *p* = 4 · (1 + 2 · 2 + 5 · 2) = 60 extracted phenotypic descriptors collected in The columns of this matrix are centred and normalized to unit variance prior to further analysis.

### (b) Mimicry-related phenotype

We use SDA to identify the linear combination of phenotypic descriptors that best captures the variation in mimetic phenotypes. Under the assumption that the mimetic phenotype has been under selection to resemble the phenotypes of either *R. variabilis* in one end of the transect, or *R. summersi* in the other, we use images of seven *R. variabilis* individuals to represent one group and six imaged *R. summersi* the other group, as the training set. The *R. imitator* individuals only enter the analysis to influence the choice of regularization parameter; see details in the electronic supplementary material, S1.

In the supplementary material, we also provide results when instead using the most extreme *R. imitator* populations to represent the end populations. There are disadvantages and advantages of both of these approaches. Using the model species amounts to defining the mimicry-related phenotype in terms of similarity to those species. This is desirable when the mimicry-related phenotype is of prime interest. However, it has the disadvantage that the two model species may differ in traits not mimicked by *R. imitator*. Using the most extreme *R. imitator* populations has the disadvantage that some of the individuals may not be pure mimetic forms. We obtained similar results using either of these approaches, or if we pool both the most extreme *R. imitator* populations and the model species individuals (see the electronic supplementary material). In the following, we will refer to the extreme groups as the mimicry-defining groups, independently of how they were defined. Combinations of these different ways of specifying the mimicry-defining groups and the three manifold learning methods used to quantify the phenotype are included as the electronic supplementary material, S1.

The mimicry-related phenotypic value for each individual is obtained by projecting onto the direction *β* according to equation (2.1) and will be denoted *z _{i}* for the

*i*th individual. The values are scaled linearly such that the average value for each of the model species is −1 and 1, respectively.

Grouping the individuals by location and ordering them along the transect from south to north (figure 1), a boxplot summarizing the mimicry-related phenotypic values as a function of location can be seen in figure 3. Note that the first half of the locations tend to have a value similar to *R. summersi*, whereas the other half have values closer to *R. variabilis*. Chipesa and Callanayacu have phenotypic values that are more intermediate and with relative high variances. Note that the ordering of the locations on the *x*-axis does not correspond to the actual geographical distances.

### (c) Simulation studies

We evaluated the accuracy of the method for determining the number of genes underlying a quantitative phenotype presented in §3 on simulated data for different values of the heritability (see Methods section in the electronic supplementary material). The accuracy is evaluated under a variety of scenarios constructed by: (i) varying the true number of genes *K*, (ii) sampling the admixture proportions from a uniform or a bimodal distribution, and (iii) adding white noise to the admixture proportions. The heritability was varied by simulating data for 1000 different values of The graphs in figure 4 show kernel density estimates of the median-likelihood ratios, 5th and 95th percentiles of the hypothesis of *K* = 1, 2, 3 or 4 genes versus the alternative hypothesis of the true number of genes under which the data are simulated. Thus, a negative likelihood ratio means that the correct number of genes is selected. Simulated data and likelihood ratios for a few of these thousands of simulations can be found in the electronic supplementary material, S3.

Generally, the chance of accurate estimation is reduced when: (i) the true number of genes is high, (ii) the heritability decreases, or (iii) the sample size decreases. A measure of confidence in the inference can be obtained by bootstrapping individuals, using the likelihood ratios comparing different hypotheses as statistics. However, if the estimates of *f* are very noisy, there tends to be a systematic bias towards a higher number of genes for intermediate heritabilities. The effect of this can be seen in figure 4*d*. Using a bootstrap test, we find −6.63 and 0.30 as the 5th and 95th percentiles of the likelihood ratio associated with the null hypothesis of *H*_{0} : *K* = 3 versus *H*_{1} : *K* = 4, for the scenario with a heritability of approx. 0.85, despite the fact that *K* = 3 is the true number of underlying genes. Thus, sensitivity to estimation variance in the admixture proportion must be kept in mind when applying this likelihood model.

### (d) Number of genes underlying the mimicry phenotype

The above-described likelihood model was used to estimate the number of genes underlying the quantified phenotype in *R. imitator*. We use 1000 bootstrap replicates to obtain a distribution of likelihood ratios between different alternative models. The bootstrap is performed by sampling individuals with replacement. First, a bootstrap distribution of the mixture proportions for each individual is obtained using the available 285 samples. We take into account uncertainty in the estimation of *f*, by, for each simulation, re-estimating *f* (see the electronic supplementary material, S2) by also bootstrapping microsatellite loci within each individual.

The maximum-likelihood values of *K*, for *K* = {1, …, 5} were then determined for each replicate in a separate bootstrap experiment using the 179 samples with the genetic and phenotypic data available. Figure 5*a* shows a boxplot of the distribution of likelihood ratios associated with the hypothesis *H*_{0} : *K* = *k* for *k* = 1, 2, 3, 4, 5, against the alternative hypothesis of *H*_{A} : *K* = 1 and figure 5*b* the proportion of bootstrap replicates in which each model obtained the highest likelihood value. This proportion can be interpreted as a measure of statistical confidence. In the electronic supplementary material, S4 the full distribution of likelihood ratios associated with the test of *H*_{0} : *K* = 2 against *H*_{A} : *K* = 1 can be seen.

The maximum log-likelihood values are numerically highest for *K* = 1, and in the vast majority of the runs, this model is selected as the most likely. The point estimates of the parameters for the hypothesis of *K* = 1 are [*μ*_{0}, *μ*_{1}, *μ*_{2}, *σ*_{e}] = [0.882,0.071, − 0.855,0.274]. The *p*-value associated with different model comparisons is shown in table 1.

Overall, a model with one (*K* = 1) or two genes (*K* = 2) seems to fit the data best, whereas three genes (*K* = 3) cannot be rejected. The mimetic phenotype, as measured here, is probably mostly influenced by one or two genes of major effect. No combinations of the three different ways of defining the end populations suggest more than three genes. Using alternative, nonlinear, manifold learning algorithms to quantify the mimicry-related phenotype (see the electronic supplementary material, S1), only a single combination (KDA with the model species defining the end groups) cannot reject *H*_{0} : *K* = 4 versus the alternative of *H*_{1} : *K* = 3 with a *p*-value of less than 0.05.

## 6. Discussion

We have here developed an automated procedure for characterizing complex phenotypes from images. We believe that this method, or related methods, could be of use in many systems where images are available for complex phenotypes. Automated extraction of phenotypic descriptors reduces the subjective biases that may occur when measurements are taken manually and allows for reproducibility of results. While other biases may be introduced through the choice of an image capturing system, lighting conditions and/or choice of descriptors, we believe these to be easier to identify and overcome. We note that such image analyses open up the possibility for a variety of statistical analyses of phenotypes, and their correlations, not pursued here. In this paper, we use the image analyses to define a quantitative measure of the mimetic phenotype in a transition zone between morphs of *R. imitator*. Using a new method for estimating the effective number of genes affecting this phenotype, we show that the phenotype we measure is likely to be controlled by one or two, or at most three, genes of major effect, and is very unlikely to be affected by many major effect genes. However, there could be substantial phenotypic variation controlled by other genes, but not captured by our quantitative measure of mimetic phenotype.

The fact that we have identified a measure of mimetic phenotype that is controlled by a few genes suggests that future studies aimed at mapping this phenotype have a relatively high probability of succeeding. It is substantially easier to map the genes underlying a phenotype controlled by just one or a few genes, than a phenotype controlled by many genes. The phenotype defined here would be useful for such mapping studies.

We can compare our results with *Heliconius* butterflies where the genetic basis of Müllerian mimicry is better understood. In *Heliconius erato*, the transition between the ‘postman’ and the ‘rayed’ morphs in the well-studied hybrid zone near Tarapoto, Peru is controlled by three loci of major effect, whereas in the co-mimetic *Heliconius melpomene*, the same mimetic shift (postman to rayed) is controlled by five loci [44]. In another example, the polymorphism in *Heliconius cydno alithea* in western Ecuador is controlled by two unlinked loci, one that controls colour (white/yellow) and one that controls pattern (presence/absence of melanin in a specific region of the forewing) [45]. In poison frogs, the genetic basis of colour variation is less well understood. Early crossing studies in *Oophaga pumilio* [46] suggested that pattern is probably controlled by a single locus with a dominant melanin-producing allele, whereas colour may be polygenic or controlled by a single locus with incomplete dominance. However, unlike *O. pumilio*, in which a major axis of variation in pattern is presence/absence of melanin, all known populations of *R. imitator* possess melanin on the dorsum, legs and venter. Thus, a more relevant task in the *R. imitator* system would be identifying the gene or genes that influence the spatial distribution of melanin rather than its presence or absence. Finally, in a field pedigree study [47], it was suggested that the red/yellow polymorphism in a population of *O. pumilio* was controlled by a single locus where red coloration was completely dominant over yellow. Thus, our estimates of one to three genes controlling the mimetic phenotype of *R. imitator* are fairly comparable to other systems.

The method we have developed for identifying the number of genes controlling a phenotype obtains its information from the degree of clustering of phenotypes and from the dependence of variance in the phenotype on the admixture fraction. It can, as illustrated here, be used to distinguish between one, two or several genes, but is not expected to perform well in estimating the exact number of genes, when many genes are involved. In the presence of many genes, the information regarding clustering of phenotypes is lost. We note that the method can be sensitive to the precision in the estimate of the admixture fraction, and results of the method should be interpreted accordingly. Implementations of the presented methods are publicly available at https://github.com/schackv.

## Ethics statement

All research was conducted following an animal use protocol approved by the Animal Care and Use Committee of East Carolina University.

## Funding statement

Field research and collections were supported by an NSF DDIG (1210313) grant awarded to E.T. and K.S., and a National Geographic Society grant awarded to K.S. (8751-10).

## Acknowledgements

Research permits were obtained from the Ministry of Natural Resources (DGFFS) in Lima, Peru (authorizations no. 050-2006-INRENA-IFFS-DCB, no. 067-2007-INRENA-IFFS-DCB, no. 005-2008- INRENA-IFFS-DCB). Tissue exports were authorized under Contrato de Acceso Marco a Recursos Geneticos no. 0009-2013-MINAGRI-DGFFS/DGEFFS, with CITES permit number 003302.

- Received August 7, 2014.
- Accepted March 26, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.