Alternative forms for genomic clines

Understanding factors regulating hybrid fitness and gene exchange is a major research challenge for evolutionary biology. Genomic cline analysis has been used to evaluate alternative patterns of introgression, but only two models have been used widely and the approach has generally lacked a hypothesis testing framework for distinguishing effects of selection and drift. I propose two alternative cline models, implement multivariate outlier detection to identify markers associated with hybrid fitness, and simulate hybrid zone dynamics to evaluate the signatures of different modes of selection. Analysis of simulated data shows that previous approaches are prone to false positives (multinomial regression) or relatively insensitive to outlier loci affected by selection (Barton's concordance). The new, theory-based logit-logistic cline model is generally best at detecting loci affecting hybrid fitness. Although some generalizations can be made about different modes of selection, there is no one-to-one correspondence between pattern and process. These new methods will enhance our ability to extract important information about the genetics of reproductive isolation and hybrid fitness. However, much remains to be done to relate statistical patterns to particular evolutionary processes. The methods described here are implemented in a freely available package “HIest” for the R statistical software (CRAN; http://cran.r-project.org/).


Introduction
Hybrid zones are natural experiments offering unique insights into evolution (Endler 1977;Hewitt 1988;Harrison 1990;Buerkle and Lexer 2008). In addition, hybridization is common in nature, and constitutes an important phenomenon impacting the evolution of diversity and novelty (e.g., Anderson 1948;Anderson and Stebbins 1954;Whitham 1989;Harrison 1993;Hewitt 2001;Arnold 2006;Arnold and Martin 2009). Hybrids and hybrid zones bring many new combinations of alleles together simultaneously, potentially leading to rapid evolution of multilocus novelties that would be difficult to evolve on a locus-by-locus basis (Rieseberg and Linder 1999;Rieseberg et al. 2003;Arnold 2006;Gompert et al. 2006;Mavarez et al. 2006;Fitzpatrick and Shaffer 2007a). Recombination in hybrid populations allows different loci to evolve differently depending on their linkage relationships and functional interactions with other loci.
An important question is to what extent different genes behave as members of a single "coadapted gene complex" (Dobzhansky 1937;Mayr 1942;Michel et al. 2010), as parts of a few coevolving "genomic islands" (Turner et al. 2005;Nachman and Payseur 2012;Nosil and Feder 2012), or as free agents establishing locus-specific patterns of variation according to their own particular effects on organismal performance (Dawkins 1976;Wu 2001;Morjan and Rieseberg 2004). In nature, hybrid zones are recognizable because many features of the interacting organisms show concordant, narrow clines, consistent with limited exchange between coadapted gene pools (Key 1968;Harrison 1993;Butlin 1998). Early hybrid zone theory suggested different genes would show different patterns of gene flow and introgression (Barton 1979; Barton and Bengtsson 1986). Barton (1983) and Baird (1995) showed how groups of loci can become synergistically "coupled" (depending on the ratio of selection to recombination) to form strong barriers to gene exchange. With increasingly sophisticated molecular tools for evaluating genome-wide patterns of variation, evidence of genomic heterogeneity in patterns of gene flow is accruing at all levels from locally adapted populations to highly differentiated taxa (e.g., Arnold 2006;Yatabe et al. 2007;Nosil et al. 2008;Fitzpatrick et al. 2009;. One promising way to evaluate variation among loci in rates and patterns of gene exchange is statistical analysis of genomic clines. Genomic cline analysis explicitly compares allele or genotype frequencies of each locus (or locus-specific ancestry) to a genome-wide average representing the genomic ancestry of an individual or population (Gompert and Buerkle 2011). The approach was pioneered by Szymura and Barton (1986) as a complementary alternative to geographic cline analysis, where genetic data are evaluated against spatial coordinates or distance. Geographic cline analysis is not always appropriate, for example, in mosaic hybrid zones where hybridizing taxa segregate by habitat at a finer grain than their overlapping geographic ranges (Harrison and Rand 1989;Howard et al. 1993), broadly admixed populations such as humans in North America (Parra et al. 1998), captive livestock herds (Musani et al. 2006), introduced species (Hansen et al. 2001;Fitzpatrick and Shaffer 2007b), or other dynamic hybrid zones with patchy introgression (Machol an et al. 2011).
The basic idea of genomic clines is that, in a hybrid zone or hybrid population between parental populations P1 and P2, each individual (or deme) can be described by their genome-wide mean ancestry S = proportion of nucleotides inherited from P1. If all genes followed a single underlying pattern, then S is their mutual expected value of locus-specific ancestry: p i = proportion of copies of locus i inherited from P1. Note that shared ancestry is not necessarily equivalent to shared state. For a diagnostic marker (fixed for different alleles in P1 and P2), p i is the frequency of P1 alleles. For nondiagnostic markers, inferences about ancestry must be made from observations about shared allelic states. The goal of genomic cline analysis is to evaluate locus-specific deviation from the expectations E(p i ) = S. Particular patterns of deviation can help identify genomic regions experiencing directional selection, hybrid dysfunction, or hybrid vigor (Szymura and Barton 1986;Gompert and Buerkle 2011).
Until now, locus-specific deviations from genome-wide ancestry have been quantified using the polynomial function suggested by Szymura and Barton (1986), and multinomial logistic regression (Lexer et al. 2007;Gompert and Buerkle 2009). Barton's approach is implemented in the package "Analyse" under the name "Concordance" (Barton and Baird 1995). The Barton function has two properties that might be undesirable for genomic clines (Gompert and Buerkle 2011). First, it can be greater than one or less than zeroimpossible values for a proportion or probability. Second, even within the interval [0,1], the Barton cline is not necessarily monotonically increasing: The curve can have an intermediate local maximum and minimum, which is unexpected albeit not impossible in light of population genetics theory. Gompert and Buerkle (2011) proposed to overcome these undesirable features by splicing in flat lines in an ad hoc manner. Logistic regression-based approaches satisfy the challenge of modeling probabilities on the interval [0,1], but are less flexible in terms of the form of the fitted curves and assume the independent variable ranges from negative to positive infinity (McCullagh and Nelder 1989). The latter assumption is explicitly violated in genomic cline analysis, where the independent variable is also a proportion on [0,1]. Here, I evaluate two alternative functions that overcome these problems, one phenomenological and the other derived from population genetic theory. I also implement multivariate outlier detection as an alternative to previous hypothesis testing approaches that confound selection and drift (Gompert and Buerkle 2009;Machol an et al. 2011), and simulate hybrid zone and admixture dynamics to assess the effects of different modes of selection.

Methods
Deriving the logit-logistic cline model Bazykin (1969) used the continuous-space diffusion model for population genetics (Fisher 1937) to show that the expected form for a cline caused by heterozygote disadvantage is the familiar logistic function of Richards (1959): where m i is the cline center for locus i (spatial position of the inflection point) and b i is the slope of the curve at the inflection point, determined in the model by the strength of selection against heterozygotes and the average dispersal distance (Bazykin 1969). Investigators often consider the cline width 1/b as a fundamental description of a cline and the same form has been used to describe clines maintained by different kinds of selection or even neutral clines where the cline width depends only on dispersal and the time since secondary contact (Slatkin 1973;Endler 1977;Barton and Gale 1993;Guedj and Guillot 2011). If we assume the genome average ancestry follows a geographic cline an expression for p i in terms of S arises on rearranging equation 2 as x = (1/b) logit(S) + m and substituting for x in equation 1 to get or, define u i = b i (m i À m) indicating a relative difference in cline position and m i = b i /b indicating the relative slope of p i and This is the "logit-logistic" function describing a genomic cline for a given locus i in terms of the average genome-wide ancestry S. The parameters can be estimated from data (see below) for a sample of loci as a way to describe multilocus clines and potentially identify exceptional loci, even in situations where the geographic model does not literally apply (e.g., in mosaic hybrid zones or within single hybrid populations) or takes a more complicated form than the simple logistic function, for example, when clines are asymmetric or stepped (Barton 1983;Baird 1995;Porter et al. 1997).
The logit-logistic function follows the constraints appropriate for a relationship between two proportions (P i ∊ [0,1] and S ∊ [0,1]). It also conforms to the definition of S as the mean p i , in that when S = 0, all p i = 0, and when S = 1, all p i = 1. The genomic cline for p i deviates from the mean when the ratio of slopes m i deviates from 1.0 and/or the relative cline center u i differs from zero.

Alternative cline forms
The Barton cline (Szymura and Barton 1986) relates expected ancestry at a particular locus p i to genome-wide ancestry S as The coefficients a i and b i describe deviations in cline center and steepness relative to perfect concordance (p i = S when a i = b i = 0). Locus-specific ancestry p i (Barton) = 0 when genome-wide ancestry S = 0 and p i (Barton) = 1 when S = 1. The Barton cline is not strictly monotonic (two local extrema are possible when the absolute value of either coefficient is large) and can take values outside [0,1], undesirable for a function describing probabilities. Gompert and Buerkle (2011) introduced ad hoc splicing of flat lines into the function to remove these undesirable features. Given that local extrema might be real features of a given dataset, I adopt the splicing of horizontal lines at 0 and 1, but allow nonmonotonicity in my analysis (as in Barton and Baird 1995). Among the functional forms compared here, the possibility of nonmonotonicity is a unique feature of the Barton function that might make it the best model for certain datasets.
The regularized incomplete beta function is a strictly monotonic, but otherwise very flexible function that also describes the cumulative distribution function of a beta random variable.
where l and m are location and shape parameters related to the more familiar beta shape parameters as a = lm and b = (1 À l)m (Kruschke 2010). The B(l, m) parameterization is useful for genomic clines because l plays a similar role to u in the logit-logistic and a in the Barton cline and m a similar role to v and b. Like the logit-logistic cline, the beta cline is zero when S = 0 and one when S = 1, it is strictly monotone, and p i (b) is never more than one or less than zero.
Although the logit-logistic, Barton, and beta clines have the same number of parameters to estimate, the procrustean splicing of flat lines onto the Barton function to satisfy range constraints makes it more cumbersome than the others in a way not accounted for by model selection criteria such as AIC. That is, more information is required to specify the spliced Barton model even though it does not produce estimates of more parameters (Rissanen 1978(Rissanen , 2005Gr€ unwald 2004). Moreover, I take it as an underlying principle that the true relationship between locus-specific ancestry probability and the genome-wide average is smooth (continuously differentiable) on the interval (0,1). That is, we do not realistically expect the locus-specific probability to be identically 1.0 except when the genome-wide average is 1.0; and we do not realistically expect an abrupt plateau of identical probabilities over any interval, as implied by the spliced-in flat line. Empirically, the spliced Barton function might be an adequate (or even superior) approximation to the underlying smooth relationship, but all other things being equal, a mathematical model reflecting the underlying principle of smoothness would be preferable.

Regression-based approaches
Multinomial regression has been used to describe the probability of a locus being heterozygous or homozygous as a function of genome-wide ancestry S (Lexer et al. 2007;Gompert and Buerkle 2009). These clines in genotypic frequency differ from clines in allele frequency because multinomial regression captures the heterozygosity component of hybrid genotypes. However, multinomial regression uses sigmoid functions that are somewhat less flexible than the cline models described above, and the estimated probabilities approach zero and one asymptotically instead of reaching limits as S reaches its finite limits at zero and one. Despite these limitations, multinomial regression should be considered because of the potential information gained by considering genotype rather than allele frequencies.
A simple binomial (logistic) regression could be used to describe allele frequencies p i in terms of S. This considers the same variables as the cline models above and estimates the same number of parameters. However, it shares with multinomial regression the problem that it is designed to model probabilities that approach zero and one as the independent variable approaches infinity (positive or negative), and therefore is not entirely appropriate for an independent variable on the finite interval [0,1]. From a practical standpoint, this might be desirable flexibility or an undesirable violation of constraints, depending on one's assumptions. That is, if p i is taken to represent an allele frequency rather than a true ancestry probability, then whatever the value of p i at S = 1 is the predicted frequency of the allele in a "pure" P1 population. However, if ancestry is defined as the expectation (S = E(p i )) for all L loci in the genome, the constraint should not be violated. Approaches for defining and estimating S are evolving (e.g., Alexander and Lange 2011;Gompert and Buerkle 2011), and whether the individuals in the admixture analysis are used to define and/or estimate S and p i might affect the validity of one's assumptions.
In the following analysis, I compare the three cline models (logit-logistic, Barton, and beta) and both multinomial and binomial regression in terms of their goodness-of-fit to real and simulated data and their usefulness in identifying outlier loci potentially linked to loci affecting hybrid fitness.

Fitting models to data
Genetic data are categorical, and often best analyzed as counts. I take the simple approach of considering one allele per locus: A 1 assumed fixed in P1 and absent in P2, so the expected frequency of A 1 in hybrids is the locusspecific ancestry probability P. A more general implementation for multiple or nondiagnostic alleles is feasible, but beyond the scope of this paper, which is focused primarily on the value of alternative functional forms.
The genotype of a diploid individual for the particular locus is represented as a count (number of A 1 alleles: 0, 1, or 2) and can be modeled as a random draw of two from a binomial distribution with probability of success P. Likewise, a sample of n individuals can be modeled as a random draw of 2n from the same binomial distribution. Let x be the number of A 1 alleles in a sample of size 2n. The probability of x given P is the binomial mass x p x ð1 À pÞ 2nÀx . Under Hardy-Weinberg assumptions (e.g., Hardy 1908;Hartl and Clark 1997), P is the population allele frequency of A 1 . But population-level Hardy-Weinberg equilibrium need not be assumed if P can be estimated independently for an individual or group. The genomic clines approach models variation in P among individuals, groups, or sites (e.g., eqs. 1-3), such that the binomial assumption need apply only "locally" for individuals characterized by the same S. Given K samples with varying P, the likelihood of the parameters is proportional to the product of binomials: Now let p k be given by one of the cline functions (eqs. 4, 5, or 6). Then equation 7 can be used to calculate the likelihood of model coefficients given x, n, and an estimate of S for each sample. Note that this is the same likelihood function maximized in binomial regression, with logit = (p k ) = a + bS (McCullagh and Nelder 1989).
I wrote functions in R (R Development Core Team 2011) to search for maximum likelihood estimates for each function described here. I implemented the likelihood search using the native R function "optim." For the regression-based methods, I used the functions "multinom" and "glm" to fit multinomial and binomial models, respectively (Venables and Ripley 2002). Gompert and Buerkle (2011) developed a hierarchical Bayesian method for estimating Barton cline parameters while accounting for uncertainty in parental allele frequencies and explicitly modeling locus-specific effects. This Bayesian method has been extended to account for linkage relationships and uncertainty in genotyping (Gompert and Buerkle 2012;Gompert et al. 2012a,b). Presumably, their approach could readily incorporate the cline functions reviewed here, but developing such computational tools is beyond the scope of this study. My goal is to illustrate the value of considering multiple alternatives for genomic cline analysis. I presume that the relative performance of the functional forms under the likelihood framework predicts relative performance under a Bayesian framework because the additional uncertainty accounted for (owing to genotyping error and uncertainty in parental allele frequency estimates) would be identical, no matter which equation is being fit. This might not be true in the case of the multinomial, which relies on accurate discrimination of heterozygotes and homozygotes; genotype frequency estimates might be prone to greater error than allele frequency estimates.

Outlier detection
The major goal of genomic cline analysis is usually to identify markers affected by selectionthat is, those linked to genes contributing to hybrid dysfunction, hybrid vigor, or local adaptation (e.g., Lexer et al. 2007;Gompert et al. 2012a). However, it is insufficient to identify candidate markers by rejecting the na€ ıve null hypothesis that locus-specific ancestry should match the genome-wide average (p i = S). This is because genetic drift alone will generate real variation among loci. This problem has been acknowledged by previous investigators (e.g., Gompert and Buerkle 2009). But a generally acceptable solution remains elusive.
Traditional hypothesis tests based on likelihood (Barton and Baird 1995;Machol an et al. 2011) or randomization (Gompert and Buerkle 2009) consider whether a model fitted to a given locus deviates from the null hypothesis more than expected from sampling alone. But to identify markers that differ from the average by more than expected from drift and sampling, we need to assess how the true distribution of parameters might be affected by drift. It appears that no appropriate theory has been described for genomic clines, but it would likely depend on difficult to measure factors such as population size, dispersal behavior, and time since secondary contact (Hartl and Clark 1997). Barton (2008) and Polechova and Barton (2011) provide some relevant theory for geographic clines, but no theory-based test for a null distribution of cline parameters has been proposed. Long (1991) provides a sample-wide test for heterogeneity among markers for the special case of a single admixed population (see Fitzpatrick et al. 2009). Gompert and Buerkle (2011) identified outliers by assuming univariate normal distributions for each parameter of the Barton cline. This approach was far less prone to false positives than the methods in INTROGRESS, which test the na€ ıve null hypothesis p i = S, with no effect of drift Buerkle 2010, 2012). I take a similar (but multivariate) statistical approach here, using a traditional, general-purpose method to identify outliers in multivariate data.
Assuming the joint distribution of parameter estimates among loci is multivariate normal, the squared Mahalanobis distance D 2 of each locus is expected to be distributed as a v 2 random variable with degrees of freedom equal to the number of parameters (Johnson and Wichern 1998). A locus with D 2 greater than a specified critical value or visually deviating from a quantile-quantile plot can be declared a statistical outlier (Johnson and Wichern 1998). For automated outlier detection, I used the Bonferroni-adjusted critical P-value, but visual inspection of quantile-quantile plots might be the best recommendation for exploratory analysis of real data.
This approach (or other general-purpose multivariate methods) relies on the observed variation among markers to establish an empirical basis for outlier detection. The advantage is that the distribution evolves over time with genetic drift. So, while genetic drift should make it progressively easier to reject the na€ ıve null hypothesis for any neutral marker with each passing generation, multivariate outlier detection should remain more robust.

Simulated data
To examine how well the cline models and regressionbased methods fit data and reveal outliers, I used stochastic individual-based simulations of secondary contact between two populations under two kinds of population structure. First is a geographically structured hybrid zone where we can use decades of research on geographic clines to inform our expectations with respect to genomic clines. Second is a single admixed population with random mating where geographic clines are not relevant. For each of these scenarios, I consider the effects of immigration from pure parental populations and several kinds of locus-specific selection.
For all model runs, secondary contact was initiated as 500 individuals from parental population P1 on one half of the modeled space, and 500 individuals from P2 on the other half. I kept track of 100 unlinked loci with two alleles (a diagnostic, fixed difference between P1 and P2). Most loci were neutral, but up to four could influence hybrid fitness. These conditions represent "low coupling" in the sense that synergistic effects of many loci are not possible. If many loci affect fitness, we should expect stepped clines and less clear distinction between "normal" and "outlier" markers.

Hybrid zone model
Full details and computer code (written in R) are provided in the publicly available R-package "HIest". Here, I describe the model verbally and explain the range of parameter values investigated. The hybrid zone is modeled as a rectangle in which a diploid individual can occupy any x-y coordinate (space is continuous, 2-dimensional). Individuals are outcrossing hermaphrodites and act as the female parent of a random number of offspring drawn from a Poisson distribution with expected value determined by local density dependence and genetics. Individuals that draw a number greater than 0 draw a mate at random from all other individuals in the population according to a normal density function of their distance from the focal individual. Each offspring draws its x-y location from a bivariate normal distribution centered on its mother. For all simulations analyzed here, space ranged from À3 to +3 in both dimensions and both the dispersal and mating curves had standard deviations of 0.3 (5% of the space).
In the absence of selection, the number of offspring mothered by each individual is a Poisson random variable with expected value given by a Beverton-Holt function (Begon et al. 2006) of local density measured as the sum over all individuals weighted by a normal density function of their distance from the focal individual with standard deviation 0.2. For most simulations I used a baseline growth rate R = 2 and local carrying capacity K = 14, which resulted in a steady state of approximately 1000 individuals distributed evenly across space.
In the presence of selection, the Beverton-Holt function was simply multiplied by the individual's relative fitness (see below). This means selection was limited to the female component of fitness, effectively weakening the intensity of selection.
The hybrid zone could be closed or open. When closed, if an offspring drew a location outside of the defined space, it was reassigned to a position at the nearest edge. That is boundaries were reflecting, and no immigrants were added during the simulation. When open, if an offspring drew a location within 5% of either x-boundary, it was replaced by a pure parental genotype (P1 on one side, P2 on the other). The y-boundary was still reflecting. Thus, for an open simulation, an average of about 5% of the individuals in the hybrid zone on each side (those closest to the ends) were replaced by immigrants from the parental populations.

Admixture model
To maintain as much similarity between models as possible, I simulated a panmictic population using the same continuous space and local density-dependent reproduction, but mates were chosen at random from all other individuals in the population without regard to distance, and offspring drew their x-y coordinates from a uniform distribution covering the entire space. That is, there was no spatial correlation between mates or between parents and offspring. As with the hybrid zone model, the population could be closed or open. Again, in the latter case, an average of about 5% of the individuals were replaced by immigrants from each parental population.

Hybrid fitness
I modeled four kinds of locus-specific fitness effects. One locus could affect fitness according to an environmental gradient, one could have heterozygote disadvantage or advantage, and a pair of loci could have a Dobzhansky-Muller incompatibility (Turelli and Orr 2000;Fitzpatrick 2008). The remaining 96 loci were always neutral.
The environmental phenotype z of an individual was determined by one locus (the "E locus") and could take values of 1.0, 0.5, or 0.0 for P1 homozygotes, heterozygotes, and P2 homozygotes, respectively. The environmental fitness component for an individual was calculated as a Gaussian function of the difference between its phenotype and the value g of the environmental gradient where s E determines the strength of selection (I used values of 2 or 4). This fitness function gives P2 homozygotes the advantage where g = 0, heterozygotes the advantage where g = 0.5, and P1 homozygotes the advantage where g = 1. The environmental gradient was defined as a logistic function of the x-dimension (logit(g) = b(xÀm)) with a slope b of 1 and midpoint m at 1, 2, or -4 from the hybrid zone center (the latter case amounting to a universal advantage for the P1 allele). In the admixture model, the environmental gradient (with midpoints at 1 or 2) favored P2 alleles because most of the area was favorable to P2 homozygous phenotypes. Heterozygote advantage was modeled by assigning fitnesses to genotypes of an "H" locus: This is a simple case of the general model (Turelli and Orr 2000;Gavrilets 2004;Fitzpatrick 2008) with the incompatibility acting as a partial recessive.
selection, or combinations of environmental selection, heterozygote dysfunction, and Dobzhansky-Muller incompatibilities (Table 2). These different causes of hybrid fitness variation can lead to distinct patterns in geographic and genomic clines (Fig. 1). Environmental selection (or universal advantage of an allele from one parental lineage) can cause geographic and genomic clines to be displaced or result in fixation of one allele (Fig. 1A, B, F, G). Heterozygote disadvantage can cause a locus-specific cline to be steeper than average in both geographic and genomic analyses (Fig. 1C, H). Heterozygote advantage causes a relatively shallow geographic cline and an "inside-out" genomic cline (Fig. 1D, I). Dobzhansky-Muller incompatibilities tend to result in displacement of each partner locus in opposite directions in both kinds of analysis (Fig. 1E, J). The same kinds of genomic cline patterns arise under the same kinds of fitness models whether geographic structure is present or not (hybrid zone vs. admixture model), but the expected variance among neutral markers was higher in the hybrid zone model (Figs. 2 and 3).
The examples depicted in Figure 1 also illustrate the hitchhiking effect of each kind of selection; at least for the strong selection simulated here, spatial clines for the unlinked neutral loci were strongly distorted relative to the no-selection case, except for the Dobzhansky-Muller incompatibility case, for which the neutral clines (black  are the slope and midpoint of the environmental gradient; s E is the strength of selection on the environmental phenotype (eq. 8); h is heterozygote disadvantage (eq. 10); s H is heterozygote advantage (eq. 9); d is the strength of the Dobzhansky-Muller incompatibility (Table 1). lines in Fig. 1E) are indistinguishable from the no-selection case (not shown).

An empirical example
To illustrate analysis of a real dataset, I used published data from an expanding hybrid swarm formed in the 1940's when Barred Tiger Salamanders (Ambystoma tigrinum mavortium) from Texas were introduced into California and encountered the native California Tiger Salamander, A. californiense (Fitzpatrick and Shaffer 2007b;Fitzpatrick et al. 2010). The dataset includes 773 salamanders from 58 sites scored for 67 nuclear SNPs. Two of these SNPs are "ringers" having no heterozygotes in the dataset because of technical problems with genotype scoring; they are included here to assess the visibility of heterozygote deficits in genomic cline analyses. Fitzpatrick et al. (2010) showed strong evidence of genomic heterogeneity, with three markers having introgressed 95 km further into the native range than the rest (Fig. 4). Although this striking pattern is hard to miss, overall the dataset is not well suited to geographic cline analysis because there is an abrupt transition from a hybrid swarm in the Salinas Valley, where breeding sites vary widely in mean S without much relationship to Genomic cline analysis offers a satisfying alternative with potential to reveal important patterns of variation that have not sorted out along a geographic gradient. The data consist of individual diploid genotypes of markers presumed to be diagnostic (based on analyses of reference populations). Hybrid Ambystoma can have 0, 1, or 2 introduced (A. t. mavortium) alleles. Analysis of nondiagnostic markers will require an additional step relating observed genotypes to allelic ancestryp i (Gompert and Buerkle 2011).
I analyzed the data using the R-package INTROGRESS (Gompert and Buerkle 2010), using their randomization test (assuming negligible genetic drift since secondary contact) to identify markers potentially linked to loci affecting hybrid fitness. I also used the stand-alone program bgc (Gompert and Buerkle 2012) to fit the spliced Barton function and identify outliers using Gompert and Buerkle's (2011) univariate Bayesian method. Then I fit each of the cline and regression models described here and used multivariate outlier detection to identify candidate markers. Fitting was done using allele or genotype counts and genome-wide average S per locality to account for nonindependence of salamanders sampled from the same breeding pond. Given the likely importance of genetic drift over the past 60 years (20-30 generations), I predicted that the outlier detection methods would be more conservative than INTROGRESS.

Simulated data
Testing the na€ ıve null hypothesis As expected for data influenced by genetic drift, the na€ ıve null hypothesis (p i = S) was relatively easy to reject for neutral markers (Table 3). The parametric test in INTRO-GRESS Buerkle 2009, 2010) flagged 99 to 100% of the markers in no-selection simulations as significantly deviating from the null hypothesis (Bonferroniadjusted critical P ≤ 0.05/100 = 0.0005), while their permutation test was somewhat more conservative. For comparison, I also tested the na€ ıve null hypothesis using traditional likelihood ratio tests for the fitted cline models; results were virtually identical to the Gompert and Buerkle (2009) permutation test (Table 3). As expected, detection rates increased with the influence of drift over time and in closed populations (Table 3). These results strongly caution against na€ ıve null hypothesis tests (assuming zero drift) for identifying candidate markers. As pointed out by Long (1991) in a similar context, the effects of selection and drift are confounded in these tests.

Goodness-of-fit
Comparing models strictly in terms of how well the fitted curves correspond to the data, the Barton and beta clines were often best (Table 4), but multinomial or binomial regression were sometimes better in closed population simulations. This is probably explained by the lack of "pure" (S = 0 or 1) individuals anchoring the curves at each end of the ancestry spectrum. Cline models always fit best in open populations. The relative goodness-of-fit tended to change over time since secondary contact, especially in the hybrid zone scenario, where the beta cline fit best early on, but the Barton cline fit best in later generations when the hybrid zone was open to immigration, and multinomial regression fit best in later generations when the hybrid zone was closed (Table 4). Representative examples are illustrated in Figures 2 and 3.

Outlier detection
On the other hand, comparing models in terms of how well they expose exceptional loci, the logit-logistic had the best combination of precision (best ratio of true positive to false positive results) and sensitivity ( Table 5). The Barton cline had the lowest sensitivity (it missed more selected loci than any other model) and multinomial regression had rather high false-positive rates (low precision). The impact of genetic drift on cline parameters and the robustness of multivariate outlier detection are illustrated in Figure 5. As expected in a closed hybrid zone, locusspecific clines become increasingly variable over time. Because outlier detection uses the empirical distribution of parameter estimates, variation owing to drift alone did not result in statistical outliers in this example. In the example (Fig. 5), there are many loci at generations 50 and 100 that would be outliers if compared to the distribution of clines at generations 10 or 25, but appear statistically normal in their proper contexts.

Ambystoma hybrid swarm
The permutation test from INTROGRESS flagged 65 of 67 markers as deviating from the na€ ıve null hypothesis (Bonferroni adjusted critical P ≤ 0.05/67 = 0.00075). Multivariate outlier detection using the fitted binomial regression, logit-logistic cline, and beta cline models flagged only the previously identified "superinvasive" markers (cm6E11, cm12C11, cm23C6 Fig. 6). Analyses Each number is the average percent of loci best fit by each model (according to AIC c ) across all simulations. Time is the number of elapsed generations since secondary contact, "Immigration" is either closed or open to immigration from pure parental populations, and "Structure" is either the hybrid zone (HZ) or admixture (AD) model. with the Barton cline (using the Bayesian method of bgc or my likelihood implementation in Hiest) failed to detect cm12C11 as an outlier. Outlier detection using multinomial regression identified the three superinvasive markers and the two "ringers" with no heterozygous genotypes. Although the latter detections add nothing to our biological knowledge of the example system, they do underscore the unique ability of the multinomial regression approach to describe genotype frequency variation.

Discussion
Comparative analysis of genomic clines yields important insights into hybrid zones, admixture dynamics, and genomic heterogeneity, particularly when geographic clines are not applicable. Although the genomic cline approach was first proposed over 25 years ago (Szymura and Barton 1986), the approach has suffered from a lack of alternative models, little information on expected signatures of different kinds of selection, and no rigorous method for identifying exceptional markers while accounting for genetic drift (until Buerkle 2011, 2012). I have addressed these problems by introducing the beta and logit-logistic cline models, simulating hybrid zone and admixture dynamics to investigate the effects of different kinds of selection on genomic cline shape, and implementing a well-known multivariate outlier detection method (similar to Gompert and Buerkle's (2011) univariate approach). At least for the conditions examined here, the logit-logistic cline model is the best for identifying markers of interest. The beta and logit-logistic cline models, introduced for the first time here, overcome some theoretical shortcomings of previous approaches. In particular, they appropriately model the relationship between two proportions or probabilities (variables defined only on the finite interval [0,1]) and meet the constraint (imposed by the definition of genomic ancestry) that ancestry probability for a given Figure 5. Effects of genetic drift on genomic cline parameters and multivariate outlier detection. Data are from a simulation of the closed hybrid zone model with strong heterozygote disadvantage (h = 0.8). Black lines and symbols represent 99 unlinked neutral markers; red represents the selected marker. The fitted logit-logistic model is shown, but results were similar for other cline models. The selected locus was a statistically significant outlier in generations 25, 50, and 100. There were no significant outliers in generation 10 (a = 0.0005). The line in each Q-Q plot illustrates the expectation of equality between empirical quantiles of the mahalanobis distance (D 2 ) and quantiles of the v 2 distribution with 2 degrees of freedom. locus must be zero or one (respectively) when the genome-wide ancestry probability is zero or one (respectively). The Barton cline model meets the latter constraint, but allowed the dependent variable to trespass above one or below zero without ad hoc truncation (Gompert and Buerkle 2011). Multinomial regression properly models a dependent variable on [0,1], but assumes the independent variable can stretch from negative to positive infinity.
Although beta distributions can be used to model allele frequency variation in structured populations (Balding and Nichols 1995;Pritchard and Donnelly 2001), the beta cline is phenomenological, as are the Barton cline and regression approaches. In contrast, the logit-logistic cline function arises from simple population genetic theory for geographic clines (Bazykin 1969). It is interesting that the logit-logistic was not commonly the best model for describing data in terms of goodness-of-fit (Table 4). Perhaps this is not surprising given that the simple geographically sigmoid model from which it was derived is inaccurate when several loci affect fitness in the center of a hybrid zone (Barton 1983;Baird 1995). In particular, when several linked loci affect fitness, there is a "coupling" effect where multilocus clines are steeper and more coincident in the center of a hybrid zone (Barton 1983;Baird 1995). This has been described as a "step," where the logistic function describing cline shape in the cline center is discontinuous with more gradual "tails of introgression" on either side (Barton and Bengtsson 1986;Barton and Gale 1993). Further development is needed to treat this, perhaps more realistic model. However, as a practical matter, the logit-logistic was most effective for identifying outliers caused by natural selection in my simulations (Table 5).
As for geographic clines, there is no one-to-one correspondence between genomic cline shape and mode of selection. For example, previous work showed that whether heterozygote disadvantage or epistatic hybrid dysfunction reliably cause sigmoid genomic clines depends on population structure (the influence of dispersal and drift) in addition to selection intensity (e.g., Gompert and Buerkle 2011;Gompert et al. 2012b). A few qualitative generalizations are supported by those studies and the present results. Selection against heterozygotes tends to steepen genomic clines, while heterozygote advantage flattens them out. Dobzhansky-Muller incompatibilities tend to generate complementarily displaced pairs of clines. These might be impossible to distinguish from displacement caused by directional selection or an offset environmental gradient. As noted for geographic clines (Barton and Gale 1993;Kruuk et al. 1999), an environmental gradient can generate clinal patterns indistinguishable from environment-independent selection against heterozygotes. Finally, genomic cline analysis is inherently relativistic; if many markers are associated with fitness in similar ways, they will not be seen as statistical outliers. Nevertheless, genomic clines provide an excellent way to screen for exceptional loci, and joint consideration of alternative cline forms offers valuable perspective on hybrid zone dynamics and patterns of genomic heterogeneity. Future work should incorporate these alternative functional forms with methods accounting for uncertainty associated with nondiagnostic markers.