Simple inheritance of color and pattern polymorphism in the steppe grasshopper Chorthippus dorsatus

The green–brown polymorphism of grasshoppers and bush-crickets represents one of the most penetrant polymorphisms in any group of organisms. This poses the question of why the polymorphism is shared across species and how it is maintained. There is mixed evidence for whether and in which species it is environmentally or genetically determined in Orthoptera. We report breeding experiments with the steppe grasshopper Chorthippus dorsatus, a polymorphic species for the presence and distribution of green body parts. Morph ratios did not differ between sexes, and we find no evidence that the rearing environment (crowding and habitat complexity) affected the polymorphism. However, we find strong evidence for genetic determination for the presence/absence of green and its distribution. Results are most parsimoniously explained by three autosomal loci with two alleles each and simple dominance effects: one locus influencing the ability to show green color, with a dominant allele for green; a locus with a recessive allele suppressing green on the dorsal side; and a locus with a recessive allele suppressing green on the lateral side. Our results contribute to the emerging contrast between the simple genetic inheritance of green–brown polymorphisms in the subfamily Gomphocerinae and environmental determination in other subfamilies of grasshoppers. In three out of four species of Gomphocerinae studied so far, the results suggest one or a few loci with a dominance of alleles allowing the occurrence of green. This supports the idea that brown individuals differ from green individuals by homozygosity for loss-of-function alleles preventing green pigment production or deposition.


Introduction
Our simulations explore different inheritance mechanisms that might explain the results of a half-sib-full-sib breeding design implemented with the steppe grasshopper Chorthippus dorsatus. The species occurs in four distinct color morphs: uniform brown (individuals without green areas, B in the following), dorsal green (green dorsally and brown laterally, D), uniform green (clear green coloration on both dorsal and lateral sides, G), and lateral green (green laterally but brown dorsally, L). Morph frequencies at that sampling site are as follows: We implemented a paternal half-sib-full-sib breeding with a total of 46 different fathers ('sires') and 134 mothers ('dams') and scored the color morph of 1,415 offspring. We pooled matings in nine distinct mating combinations, irrespective of which phenotype served as a dam or a sire. The following simulations generate the expected offspring morph frequencies given our specific breeding design and compare the simulation results to the empirically determined offspring frequencies.

Function for estimating allele frequencies for oligogenic models
We use simulations to find allele frequencies that fit field morph frequencies. The function simulates all combinations of allele frequencies in steps defined by AFstep and calculates the sum of squared deviations from field morph frequencies (FieldFreq). Deviations are shown in column Dev of the output, with lower values indicating a better fit. A table holding the possible genotypes and their phenotypes (genSim) needs to be provided as well as vectors of alleles names at up to four loci (allelesL1, allelesL2, allelesL3 and allelesL4, respectively). Note that allele names need to be unique across loci. To save computation time, the range of allele frequencies to be explored can be restricted using the AFmax argument, a vector (across all alleles at all loci) of maximum allele frequency per allele.

Function for estimating genotypic thresholds for polygenic models
We use simulations to find genotypic thresholds that fit field morph frequencies. The function simulates all combinations of genotypic thresholds in quantile steps defined by quantSteps and calculates the sum of squared deviations from field morph frequencies (FieldFreq). Deviations are shown in column Dev of the output, with lower values indicating a better fit. A table holding the possible genotypic values (genSim) needs to be provided as well as vectors of threshold names for up to three traits (ThNamesT1, ThNamesT2 and ThNamesT3, respectively). Note that threshold names should be unique across loci. To save computation time, the range of allele frequencies to be explored can be restricted using the minT and maxT arguments, vectors (across all thresholds at all traits) of minimum and maximum genotypic values, respectively, for each threshold. The option verbose allows that the progess is printed on screen.

Function for plotting and summaries
A function for plotting displays the observed offspring phenotype numbers by mating combination (obs) as grey bars and the simulation results (sim) as dots. The function also calculates summaries for the fit of simulations to observed offspring numbers. First, the proportion of simulations that are as extreme as or more extreme than the observed number of offspring (a p value for the probability that the observations arose by sampling variation given the mating design). Second, TRUE/FALSE values for whether the simulations resulted in all-zero values for offspring morphs that were observed in the data. A TRUE value might indicate that a particular offspring morph cannot be produced by a particular mating combination (given the inheritance rules being tested). Results are returned as a data frame.

Overview of models
We explore a number of inheritance models. Models differ in the number of genotypes they produce and in the number of dominance relationships to be determined. Furthermore, models differ in the number of parameters (allele frequencies and allele threshold) to be approximated and the minimum number of rules (i.e., dominance relationships) to be defined. We consider a parsimonious model one that explains the patterns in the data with the lowest number of parameters, fewest rules, and, preferentially, fewest genotypes. The column labelled Dominance shows the cases that are presented in the following section (for details see sections below). The two best-fitting models are Model 7 and Model 5, respectively. Three alleles at a single locus might produce a dominance hierarchy. There are at least two rules to be defined: the order of dominance among the three loci. A potential fit to the data are three alleles in dominance sequence D > B > L, producing G phenotypes from DL genotypes, D phenotypes from DB or DD genotypes, L phenotypes from LL genotypes, and B phenotypes from BB or BL genotypes. There are 3 * 3 = 9 possible genotypes (6 of them unique). gen = data.frame(Genotype=rep(NA,9), Phenotype=NA, Freq=NA) i = 1 for(a1 in c("D", "L", "B")) for(a2 in c("D", "L", "B")) { gen$Genotype[i] = paste0(a1, a2)
We use simulations to find allele frequencies that fit field morph frequencies. The first round scans in allele frequency steps of 5% to explore the full range of possibilities, the second round scans in steps of 2% and the third in steps of 1% to test for the best combination in a restricted range. The output shows a column Dev with the sum of squared deviation between predictions and field morph frequencies, where lower values indicate a better fit (the top 5 results are shown). Now parents are sampled from this population according to the implemented mating design. Offspring genotypes are sampled from parental genotypes, offspring phenotypes determined, and the results tabulated. The process is repeated nruns times. Simulation results are plotted (dots) along with the observed number of offspring (bars) per mating combination. The following table shows the proportion of simulations that are as extreme as or more extreme than the observed number of offspring (a p value for the probability that the observations arose by sampling variation given the mating design). Only mating combinations in which at least one offspring morph shows p < 0.05 are shown. The model predicts offspring morph frequencies that are significantly different from observed numbers for almost all mating combinations. For example for B/B matings, no other offspring morphs than B are possible, for D/D matings L and G offspring morphs are impossible and for G/G matings B offspring morphs are impossible. These patterns are thus incompatible with the observed data.

Model 3: Single-locus models with four alleles
Four alleles at a single locus might produce a dominance hierarchy. There are at least three rules to be defined: the order of dominance among the four loci. A potential fit to the data are four alleles in dominance sequence D > L > G > B. There are 4 * 4 = 16 possible genotypes (4 + 3 + 2 + 1 = 10 of them unique). gen = data.frame(Genotype=rep(NA,16), Phenotype=NA, Freq=NA) i = 1 for(a1 in c("D", "L", "G", "B")) for(a2 in c("D", "L", "G", "B")) { gen$Genotype (substr(gen$Genotype,1,1)=="D" | substr(gen$Genotype,2,2)=="D")] = "D" gen$Phenotype[is.na(gen$Phenotype) & (substr(gen$Genotype,1,1)=="L" | substr(gen$Genotype,2,2)=="L")] = "L" gen$Phenotype[is.na(gen$Phenotype) & (substr(gen$Genotype,1,1)=="G" | substr(gen$Genotype,2,2)=="G")] = "G" gen$Phenotype[is.na(gen$Phenotype) & (substr(gen$Genotype,1,1)=="B" | substr(gen$Genotype,2,2)=="B")] = "B" These genotypes produce phenotypes as follows. We use simulations to find allele frequencies that fit field morph frequencies. The first round scans in allele frequency steps of 5% to explore the full range of possibilities, the second round scans in steps of 2% and the third in steps of 1% to test for the best combination in a restricted range. The output shows a column Dev with the sum of squared deviation between predictions and field morph frequencies, where lower values indicate a better fit (the top 5 results are shown). Now parents are sampled from this population according to the implemented mating design. Offspring genotypes are sampled from parental genotypes, offspring phenotypes determined, and the results tabulated. The process is repeated nruns times. Simulation results are plotted (dots) along with the observed number of offspring (bars) per mating combination. The following table shows the proportion of simulations that are as extreme or more extreme than the observed number of offspring (a p value for the probability that the observations arose by sampling variation given the mating design). Only mating combinations in which at least one offspring morph shows p < 0.05 are shown.  The model predicts offspring morph frequencies that are significantly different from observed numbers for most mating combinations. For example for B/B matings, no other offspring morphs than B are possible, and from G/G, B/G and G/L matings, D offspring will not be produced at all. These patterns are thus incompatible with the observed data. More generally, changes in the dominance order will produce one phenotype (here B) that is only produced by a homozygous recessive genotype, and pure matings involving this genotype (here B/B) will only produce offspring of a single morph. There is no such mating combination in our dataset to produces exclusively one offspring morph. Alternative dominance rankings are thus also incompatible with the observed data.
PunnettSquare(gen, nloci=2) We use simulations to find allele frequencies that fit field morph frequencies. The first round scans in allele frequency steps of 5% to explore the full range of possibilities, the second round scans in steps of 2% and the third in steps of 1% to test for the best combination in a restricted range. The output shows a column Dev with the sum of squared deviation between predictions and field morph frequencies, where lower values indicate a better fit (the top 5 results are shown). Now parents are sampled from this population according to the implemented mating design. Offspring genotypes are sampled from parental genotypes, offspring phenotypes determined, and the results tabulated. The process is repeated nruns times. Simulation results are plotted (dots) along with the observed number of offspring (bars) per mating combination. The following table shows the proportion of simulations that are as extreme as or more extreme than the observed number of offspring (a p value for the probability that the observations arose by sampling variation given the mating design). Only mating combinations in which at least one offspring morph shows p < 0.05 are shown.  The model predicts offspring morph frequencies that are significantly different from observed numbers for several mating combinations. For example for B/B matings, no other offspring morphs than B are possible, for D/D matings G and L offspring are impossible and for B/L matings D and G offspring are impossible. These patterns are thus incompatible with the observed data. More generally, changes in the dominance order will produce one phenotype (here B) that is only produced by a double homozygous recessive genotype, and pure matings involving this genotype (here B/B) will only produce offspring of a single morph. There is no such mating combination in our dataset to produces exclusively one offspring morph. Alternative dominance rankings are thus also incompatible with the observed data.

Model 5: Two-locus models with two and three alleles
There might be two loci involved in color morph determination: one locus G controlling ability to produce green and one locus R that controls the regions in which green occurs. There are 3 * 6 = 18 unique genotypes and at least two rules to be defined: the order of dominance at two loci. A potentially fitting mechanisms has two alleles, G and b, at a G locus where the ability to produce green (G) is dominant over brown (b). For the region that controls green, there could be three alleles, D for dorsal green, L for lateral green and n for no green, with D and L being dominant over n, and D and L co-dominant, such that if the ability to produce green is given (by the G locus), the DL genotype produces fully green, the Dn genotype dorsal green, the Ln genotype lateral green and the nn genotype brown. gen = data.frame(Genotype=rep(NA,36), Phenotype=NA, Freq=NA) i = 1 for(d1 in c("G", "b")) for(d2 in c("G", "b")) for(l1 in c("D", "L", "n")) for(l2 in c("D", "L", "n")) { gen$Genotype We use simulations to find allele frequencies that fit field morph frequencies. The first round scans in allele frequency steps of 5% to explore the full range of possibilities, the second round scans in steps of 2% and the third in steps of 1% to test for the best combination in a restricted range. The output shows a column Dev with the sum of squared deviation between predictions and field morph frequencies, where lower values indicate a better fit (the top 5 results are shown). Now parents are sampled from this population according to the implemented mating design. Offspring genotypes are sampled from parental genotypes, offspring phenotypes determined, and the results tabulated. The process is repeated nruns times. Simulation results are plotted (dots) along with the observed number of offspring (bars) per mating combination. The following table shows the proportion of simulations that are as extreme as or more extreme than the observed number of offspring (a p value for the probability that the observations arose by sampling variation given the mating design). Only mating combinations in which at least one offspring morph shows p < 0.05 are shown.  The model predicts offspring morph frequencies that appear to be an overall good fit to the observed numbers. However, the model does not allow from G or L offspring from D/D matings, since both parents lack the L allele required to produce a uniform green phenotype. The number of G offspring from G/G matings is slightly underpredicted, for B/L matings the number of G and D offspring are slightly overpredicted and the number of L offspring underpredicted, for D/L matings the number of G offspring is slightly overpredicted and the number of L offspring underpredicted, and for D/G matings the number of D offspring from is slightly overpredicted. However, these differences seem more quantitative rather than qualitative.
Alternative allele frequencies of p G = 0.32, p b = 0.68, p D = 0.48, p L = 0.03 and p n = 0.48 show the same issues with D/D matings, an significant underprediction of G offspring from B/G matings and a significant overprediction of D offspring from D/G matings. These patterns are thus incompatible with the observed data.
PunnettSquare(gen, nloci=3) We use simulations to find allele frequencies that fit field morph frequencies. The first round scans in allele frequency steps of 5% to explore the full range of possibilities, the second round scans in steps of 2% and the third in steps of 1% to test for the best combination in a restricted range. The output shows a column Dev with the sum of squared deviation between predictions and field morph frequencies, where lower values indicate a better fit (the top 5 results are shown).
The process is repeated nruns times. Simulation results are plotted (dots) along with the observed number of offspring (bars) per mating combination. The following table shows the proportion of simulations that are as extreme as or more extreme than the observed number of offspring (a p value for the probability that the observations arose by sampling variation given the mating design). Only mating combinations in which at least one offspring morph shows p < 0.05 are shown. The model predicts offspring morph frequencies that are significantly different from observed numbers for several mating combinations. In particular the model does not allow for L offspring from D/D and D/G mating combinations (due to lack of an N allele) and it does not allow for D offspring from G/G and G/L mating combinations (due to lack of an U allele). Several other numbers are significantly over-or underpredicted. These patterns are thus incompatible with the observed data.

## GDL GDn GuL Gun bDL bDn buL bun ## GDL
We use simulations to find allele frequencies that fit field morph frequencies. The first round scans in allele frequency steps of 5% to explore the full range of possibilities, the second round scans in steps of 2% and the third in steps of 1% to test for the best combination in a restricted range. The output shows a column Dev with the sum of squared deviation between predictions and field morph frequencies, where lower values indicate a better fit (the top 5 results are shown). Now parents are sampled from this population according to the implemented mating design. Offspring genotypes are sampled from parental genotypes, offspring phenotypes determined, and the results tabulated. The process is repeated nruns times. Simulation results are plotted (dots) along with the observed number of offspring (bars) per mating combination. The following table shows the proportion of simulations that are as extreme or more extreme than the observed number of offspring (a p value for the probability that the observations arose by sampling variation given the mating design). Only mating combinations in which at least one offspring morph shows p < 0.05 are shown.

## GDl GDN Gul GuN bDl bDN bul buN ## GDl
We use simulations to find allele frequencies that fit field morph frequencies. The first round scans in allele frequency steps of 5% to explore the full range of possibilities, the second round scans in steps of 2% and the third in steps of 1% to test for the best combination in a restricted range. The output shows a column Dev with the sum of squared deviation between predictions and field morph frequencies, where lower values indicate a better fit (the top 5 results are shown). Now parents are sampled from this population according to the implemented mating design. Offspring genotypes are sampled from parental genotypes, offspring phenotypes determined, and the results tabulated. The process is repeated nruns times. Simulation results are plotted (dots) along with the observed number of offspring (bars) per mating combination. The following table shows the proportion of simulations that are as extreme as or more extreme than the observed number of offspring (a p value for the probability that the observations arose by sampling variation given the mating design). Only mating combinations in which at least one offspring morph shows p < 0.05 are shown.  The model predicts offspring morph frequencies that are significantly different from observed numbers for most mating combinations. In particular the model does not allow for D offspring from G/G and G/L mating combinations (due to lack of an N allele), it overpredicts G offspring from B/D and D/G mating combinations, it underpredicts G offspring from B/G and D/G mating combinations and it underpredicts L offspring from D/L mating combinations. The results are thus not compatible with the observed data.

## GMR GMr GwR Gwr bMR bMr bwR bwr ## GMR
We use simulations to find allele frequencies that fit field morph frequencies. The first round scans in allele frequency steps of 5% to explore the full range of possibilities, the second round scans in steps of 2% and the third in steps of 1% to test for the best combination in a restricted range. The output shows a column Dev with the sum of squared deviation between predictions and field morph frequencies, where lower values indicate a better fit (the top 5 results are shown). Now parents are sampled from this population according to the implemented mating design. Offspring genotypes are sampled from parental genotypes, offspring phenotypes determined, and the results tabulated. The process is repeated nruns times. Simulation results are plotted (dots) along with the observed number of offspring (bars) per mating combination. The following table shows the proportion of simulations that are as extreme as or more extreme than the observed number of offspring (a p value for the probability that the observations arose by sampling variation given the mating design). Only mating combinations in which at least one offspring morph shows p < 0.05 are shown.

Four-locus models Model 10: Four-locus models with two alleles (2 dominant alleles for L)
There might be four loci involved in color morph determination: One locus G with two alleles (G and b) might be responsible for the ability to produce green color, with G dominant over b. One locus D with two alleles (D and u) where the recessive allele u suppresses green color dorsally. Two loci L and M with two alleles each (L and n and M and w, respectively) where the recessive alleles n and w suppresses green color laterally when either one is in homozygous state. gen = data.frame (Genotype=rep(NA,256), Phenotype=NA, Freq=NA) i = 1 for(g1 in c("G", "b")) for(g2 in c("G", "b")) for(d1 in c("D", "u")) for(d2 in c("D", "u")) for(l1 in c("L", "n")) for(l2 in c("L", "n")) for(m1 in c("M", "w")) for(m2 in c("M", "w")) { gen$Genotype[i] = paste0(g1, g2, d1, d2, l1, l2, m1, m2) These genotypes would produce phenotypes as follows.

## GDLM GDLw GDnM GDnw GuLM GuLw GunM Gunw bDLM bDLw bDnM bDnw buLM buLw bunM
We use simulations to find allele frequencies that fit field morph frequencies. The first round scans in allele frequency steps of 5% to explore the full range of possibilities, the second round scans in steps of 2% and the third in steps of 1% to test for the best combination in a restricted range. The output shows a column Dev with the sum of squared deviation between predictions and field morph frequencies, where lower values indicate a better fit (the top 5 results are shown).  = c(G=0.36, b=0.64, D=0.44, u=0.56, L=0.18, n=0.82, M=0.12, w=0.88 Now parents are sampled from this population according to the implemented mating design. Offspring genotypes are sampled from parental genotypes, offspring phenotypes determined, and the results tabulated. The process is repeated nruns times. Simulation results are plotted (dots) along with the observed number of offspring (bars) per mating combination. The following table shows the proportion of simulations that are as extreme as or more extreme than the observed number of offspring (a p value for the probability that the observations arose by sampling variation given the mating design). Only mating combinations in which at least one offspring morph shows p < 0.05 are shown. The model predicts offspring morph frequencies that appear to be an overall good fit to the observed numbers. The number of G offspring is slighly underpredicted for B/G matings and the number of L offspring is slightly underpredicted for D/L matings. The biggest misprediction occurs for D/G matings where D offspring are over and G offspring underpredicted. Notably, there is no observed offspring that is impossibly produced by the model.

Threshold-based multi-locus models Model 11: One-trait model with morph genotypic values (D > L > G > B)
There might also be many loci influencing a single trait with trait-value thresholds determining morph identities. We here assume five independent loci (or haplotypes) with equal effect sizes to allow for some sampling variation (though there could be many more loci involved). First, we generate a large pool of individuals (controlled by nind), each with 10 alleles (2 at 5 loci). Genotypic values (GV1 ) are calculated as the average across alleles. The simulation is set up for more traits to be generated later, though genotypic values for (silent) traits 2 and 3 (GV2 and GV3 ) are constraint to zero here. nind=100000 loci = data.frame(matrix(c(rnorm(nind*10,0,sqrt(10)), rep(0,nind*20)), ncol=30)) gen = data.frame(loci, GV1=rowMeans(loci[,1:10]), GV2=0, GV3=0) gen$Phenotype = NA phenFnc = function(mydat, ThreshT1, ThreshT2, We use simulations to find genotypic thresholds that fit field morph frequencies. The first round scans in allele frequency steps of 5% to explore the full range of possibilities, the second round scans in steps of 2% and the third in steps of 1% to test for the best combination in a restricted range. The output shows a column Dev with the sum of squared deviation between predictions and field morph frequencies, where lower values indicate a better fit (the top 10 results are shown). Now parents are sampled from this population according to the implemented mating design. Offspring genotypic values are sampled from parental genotypes, offspring phenotype determined and the results tabulated. The process is repeated nruns times. Simulation results are plotted (dots) along with the observed number of offspring (bars) per mating combination. The following table shows the proportion of simulations that are as extreme as or more extreme than the observed number of offspring (a p value for the probability that the observations arose by sampling variation given the mating design). Only mating combinations in which at least one offspring morph shows p < 0.05 are shown.

Model 12: One-trait model with morph genotypic values (G > L > D > B)
There might also be many loci influencing a single trait with trait-value thresholds determining morph identities. We here assume five independent loci (or haplotypes) with equal effect sizes to allow for some sampling variation (though there could be many more loci involved). First, we generate a large pool of individuals (controlled by nind), each with 10 alleles (2 at 5 loci). Genotypic values (GV1 ) are calculated as the average across alleles. The simulation is set up for more traits to be generated later, though genotypic values for (silent) traits 2 and 3 (GV2 and GV3 ) are constraint to zero here. nind=100000 loci = data.frame(matrix(c(rnorm(nind*10,0,sqrt(10)), rep(0,nind*20)), ncol=30)) gen = data.frame(loci, GV1=rowMeans(loci[,1:10]), GV2=0, GV3=0) gen$Phenotype = NA phenFnc = function(mydat, ThreshT1, ThreshT2, We use simulations to find genotypic thresholds that fit field morph frequencies. The first round scans in allele frequency steps of 5% to explore the full range of possibilities, the second round scans in steps of 2% and the third in steps of 1% to test for the best combination in a restricted range. The output shows a column Dev with the sum of squared deviation between predictions and field morph frequencies, where lower values indicate a better fit (the top 10 results are shown). A reasonable fit to field morph frequencies can be achieved with thresholds of T B = 0.15, T D = 1.34, p L = 1.41 (though these thesholds overpredict the number of G morphs). Also note slight sampling variation between runs. ThreshT1 = c(B=0.15, D=1.34, L=1.41) gen$Phenotype = phenFnc(gen, ThreshT1, NA, NA) hist(gen$GV1, nclass=30, main="Single color-morph trait") abline(v=ThreshT1, lwd=2) text(c (-3,-0.1,1,3), 18000, c("Uniform\nbrown", "Dorsal\ngreen", "Lateral\ngreen", "Uniform\ngreen")) Now parents are sampled from this population according to the implemented mating design. Offspring genotypic values are sampled from parental genotypes, offspring phenotype determined and the results tabulated. The process is repeated nruns times. Simulation results are plotted (dots) along with the observed number of offspring (bars) per mating combination. The following table shows the proportion of simulations that are as extreme as or more extreme than the observed number of offspring (a p value for the probability that the observations arose by sampling variation given the mating design). Only mating combinations in which at least one offspring morph shows p < 0.05 are shown.

Model 13: Two-trait models
Instead of being explained by a single 'trait', morph could be considered to consist of two 'traits'. Hence, there could be many loci influencing two traits with the combination of three trait-value thresholds determining morph identities. These traits could be green on the dorsal side and green on the lateral side. We here assume five independent loci (or haplotypes) with equal effect sizes per trait to allow for some sampling variation (though there could be many more loci involved). First, we generate a large pool of individuals (controlled by nind), each with 20 alleles (2 at 5 loci, for 2 traits). Genotypic values (GV1 and GV2 ) are calculated as the average across alleles. The simulation is set up for more traits to be generated later, though genotypic values for (silent) trait 3 (GV3 ) are constraint to zero here. nind=100000 loci = data.frame(matrix(c(rnorm(nind*20,0, sqrt(10)), rep(0,nind*10)), ncol=30)) gen = data.frame(loci, GV1=scale(rowMeans(loci[,1:10])), GV2=scale(rowMeans(loci[,11:20])), GV3=0) gen$Phenotype = NA phenFnc = function(mydat, ThreshT1, ThreshT2, We use simulations to find genotypic thresholds that fit field morph frequencies. The first round scans in allele frequency steps of 5% to explore the full range of possibilities, the second round scans in steps of 2% and the third in steps of 1% to test for the best combination in a restricted range. The output shows a column Dev with the sum of squared deviation between predictions and field morph frequencies, where lower values indicate a better fit (the top 10 results are shown). A reasonable fit to field morph frequencies can be achieved with thresholds of T 1 = 0.28, T 2 = 1.75. Note slight sampling variation between runs. ThreshT1 = 0.28 ThreshT2 = 1.75 gen$Phenotype = phenFnc(gen, ThreshT1, ThreshT2, NA) par(mfrow=c(1,2)) hist(gen$GV1, nclass=30, main="Trait for color of dorsal side") abline(v=ThreshT1, lwd=2) text(c(-3,3), 18000, c("Brown", "Green")) hist(gen$GV2, nclass=30, main="Trait for color of lateral side") abline(v=ThreshT2, lwd=2) text(c (-3,3), 18000, c("Brown", "Green")) Trait for color of dorsal side gen$GV1 Now parents are sampled from this population according to the implemented mating design. Offspring genotypic values are sampled from parental genotypes, offspring phenotype determined and the results tabulated. The process is repeated nruns times. Simulation results are plotted (dots) along with the observed number of offspring (bars) per mating combination. The following table shows the proportion of simulations that are as extreme as or more extreme than the observed number of offspring (a p value for the probability that the observations arose by sampling variation given the mating design). Only mating combinations in which at least one offspring morph shows p < 0.05 are shown. nruns=1000 res = simQuantFnc(gen, dat, FemaleUq, MaleUq, nruns, phenFnc, ThreshT1, ThreshT2, NA) pVals = plotFnc(obs, res)

Model 14: Three-trait models
Instead of being explained by a single 'trait' or two 'traits', morph could be considered to consist of three 'traits'. Hence, there could be many loci influencing three traits with the combination of three trait-value thresholds determining morph identities. We here assume five independent loci (or haplotypes) with equal effect sizes per trait to allow for some sampling variation (though there could be many more loci involved). First, we generate a large pool of individuals (controlled by nind), each with 30 alleles (2 at 5 loci for 3 traits). Genotypic values (GV1, GV2, and GV3 ) are calculated as the average across alleles. nind=100000 loci = data.frame(matrix(rnorm(nind*30, 0, sqrt(10)), ncol=30), Phenotype=NA) gen = data.frame(loci, GV1=scale(rowMeans(loci[,1:10])), GV2=scale(rowMeans(loci[,11:20] We use simulations to find genotypic thresholds that fit field morph frequencies. The first round scans in allele frequency steps of 5% to explore the full range of possibilities, the second round scans in steps of 2% and the third in steps of 1% to test for the best combination in a restricted range. The output shows a column Dev with the sum of squared deviation between predictions and field morph frequencies, where lower values indicate a better fit (the top 10 results are shown). A reasonable fit to field morph frequencies can be achieved with thresholds of T B = 0.08, T G = 1.34, p L = −1.17. Note slight sampling variation between runs. ThreshT1 = 0.08 ThreshT2 = 1.34 ThreshT3 = -1.17 gen$Phenotype = phenFnc(gen, ThreshT1, ThreshT2, ThreshT3) par(mfrow=c(1,3)) hist(gen$GV1, nclass=30, main="Trait for potential to be green") abline(v=ThreshT1, lwd=2) text(c(-3,3), 18000, c("Brown", "Green")) hist(gen$GV2, nclass=30, main="Trait for color of dorsal side") abline(v=ThreshT2, lwd=2) text(c (-3,3), 18000, c("Brown", "Green")) hist(gen$GV3, nclass=30, main="Trait for color of lateral side") abline(v=ThreshT3, lwd=2) text(c (-3,3), 18000, c("Brown", "Green")) Trait for potential to be green gen$GV1 Now parents are sampled from this population according to the implemented mating design. Offspring genotypic values are sampled from parental genotypes, offspring phenotype determined and the results tabulated. The process is repeated nruns times. Simulation results are plotted (dots) along with the observed number of offspring (bars) per mating combination. The following table shows the proportion of simulations that are as extreme or more extreme than the observed number of offspring (a p value for the probability that the observations arose by sampling variation given the mating design). Only mating combinations in which at least one offspring morph shows p < 0.05 are shown. nruns=1000 res = simQuantFnc(gen, dat, FemaleUq, MaleUq, nruns, phenFnc, ThreshT1, ThreshT2, ThreshT3) pVals = plotFnc(obs, res)