A narrow window for geographic
cline analysis using genomic data: effects of age, size, and migration
on error rates
Authors: Gaston I. Jofre1,2,3, Gil G.
Rosenthal1,2
1Department of Biology, Texas A&M University, TAMU,
College Station, TX 77843, USA, 2Centro de
Investigaciones Cientıficas de las Huastecas “Aguazarca”, 16 de
Septiembre 392, Calnali Hidalgo 43230, Mexico,3Department of Biology, University of North Carolina,
Chapel Hill, NC
Correspondence
Gaston I. Jofre , Department of Biology, University of North Carolina,
Chapel Hill, NC.
Email: gastonjr@email.unc.edu
Abstract
The use of genomic and phenotypic data to scan for outliers is a
mainstay for studies of hybridization and speciation. Geographic cline
analysis of natural hybrid zones is widely used to identify putative
signatures of selection by detecting deviations from baseline patterns
of introgression. As with other outlier-based approaches, demographic
histories can make neutral regions appear to be under selection and vice
versa. In this study, we use a forward-time individual-based simulation
approach to evaluate the robustness of geographic cline analysis under
different evolutionary scenarios. We measured the effects of drift on
genetic differentiation, and on false positive and false negatives
detection using geographic clines. We modeled multiple stepping stone
hybrid zones with distinct age, deme sizes, and migration rates, and
evolving under different types of selection. We found that in young
hybrid zones, drift increases overall genomic divergence, distorts cline
shapes and increases both false positive and false negative rates. In
old hybrid zones, genomic divergence and cline distortion are higher.
Our results suggest that geographic clines are most useful for outlier
analysis in young hybrid zones with large populations of hybrid
individuals.
Keywords: outlier analysis, migration, drift, ecological genomics,
tension zone
Introduction
Hybridization occurs when two genetically distinct individuals reproduce
and generate viable offspring of mixed ancestry (Barton and Hewitt 1985;
Abbott et al. 2013). After a hybridization event, selection against
introgression maintains species boundaries (Mallet 2005; Coughlan and
Matute 2020). Alternatively, alleles favored by selection can introgress
differentially across a hybrid zone (Muirhead and Presgraves 2016;
Martin and Jiggins 2017). These different types of selection can leave
signatures that can be identified using outlier identification methods.
Cline theory has long been used to detect phenotypic introgression
across wild hybrid zones and draw inferences about selection (e.g
asymmetric sexual selection, Barton & Szymura 1989, Parsons 1993). More
recently, cline theory has been widely employed to the same end to
detect outlier regions of the genome more permeable or resistant to
introgression (Ryan et al. 2017; Stankowski et al. 2017; Lipshutz et al.
2019). Cline theory addresses changes in ancestry-allele frequencies
across a spatial gradient (Figure 1), which may covary with
environmental parameters (Haldane 1948; Endler 1977; Barton and Hewitt
1985). Geographic cline analysis compares cline positions and shapes to
detect loci potentially under selection (Szymura and Barton 1986). At
its simplest, the analysis uses two parameters, the center of the cline,
which indicates the direction of selection acting on each allele, and
the width of the cline, which indicates the maximum rate of change in
allele frequencies (Barton and Gale 1993). Neutral alleles are expected
to pass easily across the cline (neutral diffusion), and changes in
their frequencies across a spatial gradient represents the null
expectation of no selection (Barton and Hewitt 1985; Szymura and Barton
1986; Gay et al. 2008). Cline shapes that deviate from the null are
considered outliers (Payseur 2010; Stankowski et al. 2017).
Specifically, outlier loci that are involved in reproductive isolation
between parental species have reduced widths (Payseur 2010; Sciuchetti
et al. 2018); and outlier loci that indicate high rate of introgression
can have centers shifted towards either end of the cline (Barton and
Hewitt 1985).
Theoretical models show that genetic drift can distort cline shape (see
Figure 1), reducing the width and shifting the center in either
direction (Felsenstein 1975; Slatkin and Maruyama 1975; Polechova and
Barton 2011; Sciuchetti et al. 2018). Felsenstein (1975) and Slatkin and
Maruyama (1975) showed that drift causes a reduction in the expected
width, and a shift in the expected center. Polechova and Barton (2011)
and Sciuchetti et al. (2018) modeled selection against hybrid genotypes
and show that the deviation in cline center and reduction in width of a
given locus can be the result of drift alone instead of a given model of
selection. However, Polechova and Barton (2011) and Sciuchetti et al.
(2018) used one deme size. Thus, the effects of different deme sizes,
migration rates, and models of selection on hybrid zone cline shapes
remains unexplored.
Migration and drift can thus generate spurious outliers not caused by
selection. For example, small deme sizes and low migration rates
increase inbreeding (Hossjer et al. 2016), causing local fixation along
the hybrid zone. In this study, we evaluate how often genetic drift
causes neutral genomic regions to have cline shapes characteristic of
regions under selection, and how often drift modifies the cline shapes
of regions under selection. Specifically, we aim to calculate the rate
of false positive detection and the rate of false negative detection in
geographic clines using simulations in multiple hybrid zones with
different demographic histories. We model multiple stepping stone hybrid
zones under different evolutionary scenarios.
Materials and Methods
Simulation and modeling
To model how hybrid zone properties affected outlier detection, we used
Admix’em (Cui et al. (2016) ; github:
https://github.com/melop/admixem). This program generates
realistic simulations of admixed diploid populations (Schumer and
Brandvain 2016; Massey 2017; Schumer et al. 2018). For all simulations,
we modeled admixture in a cline with 19 demes and two infinitely large
source populations, each at one extreme of the cline (Figure 2;
Supporting Information Table 1). We restricted the number of hybrid
demes to 19 (similar to Felsenstein 1975), to model clines with
boundaries away from their cline centers, where one or the other allele
will be close to fixation (see Slatkin and Maruyama 1975), and while
keeping an adequate computation speed . We assumed that at generation 1
each admixed deme has an ancestry proportion p from population 2
increasing from p =0.05 to p =0.95 along the cline (see
Figure 2). The genomes consisted of ten chromosomes, each with
one randomly-drawn recombination event per arm per meiosis. This
recombination rate is similar to that suggested by Dumont and Payseur
(2008). We tracked ancestry in 11 markers randomly distributed in each
chromosome, 110 markers in total. We allowed random mating for 1,000
generations and sampled individuals in the hundredth (for young hybrid
zones), and thousandth generation (for old hybrid zones). Generations
did not overlap; at the end of each generation parental individuals were
removed after mating.
We varied the effect of drift by modifying the deme size and migration
rate. We restricted the deme sizes to 100, 500, and 1,000 individuals.
Migration rate was symmetrical and restricted to 0.1 and 0.01. We
performed 100 replicated simulations with each parameter combination,
for each type of selection (see Figure 2).
Models of selection in hybrid
zones:
A) No selection:
We ran the simulations as described above with the null expectation of
no selection with respect to ancestry.
B) Directional and underdominance selection
To determine the effects of selection on hybrids, we modeled the effects
of directional selection (with incomplete dominance) and underdominance.
We modeled directional selection (Supporting Information Figure S1, A)
as one allele having lower fitness than the other allele along the
spatial gradient. We modeled directional selection on one locus at the
middle of chromosomes three, denoted as 3.5, and another on chromosome
seven, denoted as 7.5. We modeled underdominance (Supporting Information
Figure S1, B), as heterozygous disadvantage along the spatial gradient.
We modeled underdominance on loci in the middle of chromosomes five and
nine, denoted as 5.5 and 9.5. Although underdominant selection models
are unrealistic in nature (Sedghifar et al. 2016), they portray a simple
scenario in which hybrids are unfit. We assumed additive fitness, where
relative fitness is \(f=1-x_{3.5}-x_{5.5}-x_{7.5}-x_{9.5}\) ,
where \(x\) is a function of an individual’s genotype at the specific
focal loci.
C) Bateson–Dobzhansky–Muller incompatibility
Unlike underdominant selection, Bateson–Dobzhansky–Muller (BDM)
incompatibilities (Dobzhansky 1937; Muller 1940) may be a more common
cause of reduced hybrid fitness (Coyne and Orr 2004) . We modeled the
effects of epistatic interactions by simulating a BDM interaction
between two autosomal loci, one in the middle of chromosomes three (3.5)
and another on chromosome seven (7.5). We assumed codominance in the
incompatible alleles and applied additive fitness, where relative
fitness is \(f=1-x_{3.5}-x_{7.5}\) (see selection coefficients in
Supporting Information Figure S1, C).
Patterns of genomic
differentiation using FST
To compare variation in genomic differentiation between young (100
generations) and old hybrid zones (1000 generations), and between models
of selection, we calculated Wright’s fixation index FST.
Drift generates differentiation among demes, increasingFST, while migration does the reverse (Efremov
2004). To estimate the degree in genomic differentiation, we calculated
per-locus FST among all hybrid populations
following Weir and Cockerham (1984) and using the R package pegas(Paradis 2010). We calculated the mean, median, range, and interquartile
range of the FST values from 100 loci with no
selection, under directional selection, under underdominant selection,
and in BDM interaction.
Geographic cline analysis
To study the effects of migration and drift on cline shape outliers, we
sampled 30 random individuals from each deme in generations 100 and
1,000. Admix’em generates an output file with ancestries from the
defined markers in the input files. We then used custom scripts
(https://github.com/gjofre/simulation_clines) to generate the input
files necessary for geographic cline analysis. The Metropolis–Hastings
algorithm in the R package hzar (Derryberry et al. 2014), was
used to fit geographic clines to allele frequencies for each of the 110
markers genotyped. To fit a cline, we estimated the following five
parameters in hzar : cline center \(\mathcal{y}\), cline width\(\mathcal{w}\), the ends of the cline \(P_{\min}\) and \(P_{\max}\) ,
the rate at which the cline tail decays\(\theta_{P1}\text{\ and\ }\theta_{P2}\), and the size of the cline
tails \(\mathcal{B}_{P1}\text{\ and\ }\mathcal{B}_{P2}\). For
computational speed, we fitted three models with different cline
parameter combinations. Model I estimated only center \(\mathcal{y}\)and width \(\mathcal{w}\) from the data, assumed no tails
(θ = 1and \(\mathcal{B}\) = 0), and included fixed ends
(\(P_{\min}\) and \(P_{\max}\) fixed to zero and one). Model II
estimated \(\mathcal{y}\), \(\mathcal{w}\), and \(P_{\min}\) and\(P_{\max}\) from the data, and assumed no tails. Model III was the same
as Model II, and with tail estimates θ and \(\mathcal{B}\)allowed to vary. Each model parameter was estimated using three
independent chain runs using 500,000 MCMC steps after a burn-in of
100,000 steps. The model with the lowest corrected Akaike information
criterion (AIC, Akaike (1998)) was chosen for each marker. We then used
the average estimates of all the parameters to estimate a genome-wide
average cline for each replicate.
Outlier detection analysis
We then scanned for outliers across the genome. A marker was considered
an outlier if the 95% confidence intervals for the lowest AIC value of
either \(\mathcal{w}\) or \(\mathcal{y}\) did not overlap with the
confidence intervals of the lowest AIC genome-wide average cline model
in hzar. Centers and widths of any two clines were considered coincident
(same \(\mathcal{y}\)) or concordant (same \(\mathcal{w}\)) if their
parameters overlapped with each other’s 95% confidence intervals. This
allowed us to categorize neutral regions as true negatives or false
positives, and regions under selection as true positives or false
negatives (Supporting Information Table 2).
We categorized markers under directional selection as true positives if
they displayed the cline center shifted towards the side of the
geographic gradient with higher fitness and as false negatives
otherwise. We categorized markers under underdominant selection and in
BDM incompatibilities as true positives if their center coincided with
the center of the whole genome hybrid index cline, and if their width
was lower than the width of the whole genome hybrid index cline.
Finally, we categorized markers without selection as false positives if
they were non-coincident or discordant compared to the genome-wide
cline.
We computed an error matrix for each hybrid zone (see Supporting
Information Table S2). For each replicate, we calculated the average
false positive rate (FPR) as the proportion of neutral loci that were
flagged as outliers and the false negative rate (FNR) as the proportion
of loci under each form of selection that were undetected (see SI table
2).
Statistical analyses
We fitted a linear model to observe which hybrid zone parameters
(N , m , generation on admixture, and type of selection) had
a higher effect on FST , FPR, and FNR. We first
log converted the different deme sizes, migration rates, and generations
of admixture. We then analyzed the effects of these parameters onFST , FPR, and FNR with the lm function in R. For
the FNR linear model, we omitted the data from the hybrid zones with no
selection, since there are no true positives and false negatives. We
then performed post hoc comparisons of estimated means ofFST , FPR, and FNR in each hybrid zone model
parameter.
Results
Effect of model parameters in genetic differentiation
To observe the effects of drift on genetic differentiation, we compared
the effects of each model parameter on FST values
across replicated hybrid zones. Underdominance and BDM incompatibilities
created genetic differentiation (Figure 3; Supporting Information Table
3; Supporting Information Figure S2). Deme size,
migration rate, and generation in admixture had also an effect on meanFST. Larger deme sizes and higher migration rates
decreased genetic differentiation. More generations in admixture, by
contrast, increased genetic differentiation.
We found high genetic differentiation without any need for selection.
Average FST was highest in hybrid zones with
small deme sizes, and having a higher migration rate decreasedFST further (Supporting Information Figure S2
A-C). FST values peaked in old hybrid zones with
low migration rates (m= 0.01). Drift can drastically inflate
genetic differentiation in hybrid zones with small population sizes and
generate patterns of selections in neutral genomic regions.
Under directional selection, hybrid zones with small deme sizes had
higher and wider FST distribution than hybrid
zones with bigger deme sizes (Supporting Information Figure S2 D-F). Old
hybrid zones had higher and narrower FST distributions than young hybrid zones, but only when migration rates
were high (m =0.1). Old hybrid zones with low migration rates
(m =0.01), had drastically decreased averageFST values. With reduced migration between
adjacent demes, some markers fixed along the geographic gradient,
resulting in low FST values. These results
suggest that in old hybrid zones with reduced migration rate, markers
can be under directional selection and yet show patterns of low
differentiation.
In hybrid zones under underdominant selection, or with BDMIs, small deme
sizes did not affect mean FST (Supporting
Information Figure S2 G-M). Similar to the hybrid zones under
directional selection, at 1000 generations of admixture theFST distribution became narrower in both hybrid
zones with high (m =0.1) and low migration rates (m= 0.1).
We observed again an effect of drift in old hybrid zones drastically
reducing FST values. However, this effect
happened only at a few loci and only in hybrid zones with small deme
sizes. Selection against hybrid genotypes, either via underdominance or
BDM incompatibilities, showed more consistent patterns of genomic
differentiation.
Effect of model parameters on
false positive detection
Deme size had the strongest effect on the FPR distribution, with the
highest FPR (mean FPR 0.98) in hybrid zones with small deme sizes and
low migration rates (Figure 3; Supporting Information Table 3;
Supporting Information Figure S2). More
generations in admixture also increased FPR. The mode of selection had
no effect on mean FPR values (Figure 3; Supporting Information Table 3).
Hybrid zones with small deme sizes generated the most clines with their
centers drastically shifted, showing spurious patterns of introgression
(Figure 4, A & B ; Supporting Information Video 1), which highly
inflated the FPR (Supporting Information Figure S3, A). Specifically,
old hybrid zones with small deme sizes (N= 100) and low migration
rates (m =0.01) had an average FPR of 0.98. Old hybrid zones with
large deme sizes (N = 1000) had a drastic decrease in FPR, but
only when they had migration rates of m = 0.1. Old hybrid zones
with low migration rates still had average FPR values of 0.58. The best
scenario was in young hybrid zones
with large deme sizes and high
migration rates, with mean FPR of less than 0.1.
These results suggest that when analyzing natural hybridization, young
hybrid zones will have lower risk of generating false positives than old
hybrid zones, especially with population sizes larger than N= 500.
Under the most favorable conditions (young hybrid zone with high
migration and large deme size) the false positive rate still averaged
1.35% across replicates (Figure. 4e). The risk of misidentification in
hybrid zones with small population sizes is high in both young and old
hybrid zones, regardless of migration rates. In old hybrid zones larger
than N= 500, having migration rates larger than m= 0.1 will
indeed decrease the risk of identifying false positives.
Effect of model parameters on false negative detection
Loci with underdominance and BDM incompatibilities had the highest risk
of FNR across our simulations . Hybrid zones with bigger deme sizes
showed lower average FNR, as did fewer generations of admixture.
Migration rate had a minor effect, slightly increasing the overall FNR
(Supporting Information Table 3).
Markers under directional selection showed low overall FNR which was
elevated in young hybrid zones. Drift indeed distorted cline shape, but
centers were rarely shifted to the opposite side of the gradient
(Supporting Information Figure S4 A and B, Supporting Information Video
2). In all hybrid zones with population sizes larger than N= 500,
and in all hybrid zones after 1000 generations, all markers under
directional selection were categorized as true outliers (FNR = 0;
Supporting Information Figure S4). Even when drift modifies the cline
shapes of markers under directional selection, the impact on FNR is low
and only in young hybrid zones with small deme sizes.
By contrast, markers under underdominant selection showed increased FNR
across all hybrid zones (Supporting Information Figure S3 C). In young
hybrid zones, small deme sizes ( N = 100) drastically increased
FNR, regardless of migration rates. Old hybrid zones with small deme
sizes (e.g. N =100) showed the highest FNR, and even with larger
deme sizes FNR was not substantially reduced. For example, in hybrid
zones with deme sizes of N= 1000, regardless of migration rates,
~40% of the markers had drastic shifts in their cline
centers (Supporting Information Figure S5). These results suggest that
in small and young hybrid zones, and in old hybrid zones, drift can
generate clines with patterns similar to markers under directional
selection (e.g. Supporting Information Video 3).
In markers involved in BDM incompatibilities the average FNR was
slightly lower relative to underdominant selection (Figure 3; Supporting
Information Video 4). Young hybrid zones with large deme sizes had the
lowest FNR, sp ecifically with high migration rates
(Supporting Information Figure S3 D, solid squares). On the other hand,
old hybrid zones with high migration rates, showed the highest FNR when
(N =100 and 500). The cline centers shifted to a lesser degree,
but the clines were outside the hybrid index confidence interval
(Supporting Information Figure S6 A &C). In old hybrid zones with low
migration rates, having large deme sizes decreased drastically the FNR.
These results suggest that markers involved in BDM incompatibilities
will generate fewer false negatives, compared to markers under
underdominant selection.
Discussion
Geographic cline theory (Barton and Hewitt 1985) makes powerful
predictions about how allele frequencies should change across a spatial
gradient given a range of selection dynamics. Numerous studies have used
whole-genome data with geographic cline analysis to detect outliers from
average ancestry patterns (Payseur 2010). However, drift within demes
can distort cline shapes from hybrid zones (Polechova and Barton 2011).
Our results are consistent with theoretical (Felsenstein 1975; Slatkin
and Maruyama 1975; Nagylaki and Lucier 1980) and empirical (Polechova
and Barton 2011; Sciuchetti et al. 2018) studies showing that drift can
distort spatial clines in hybrid zones. Our results suggest that drift
makes geographic cline analysis inappropriate for any hybrid zones with
small deme sizes, as well as for older hybrid zones with low migration
rates. We show that, depending on the demography (deme sizes, migration
rates and age) of a hybrid zone, drift increases false positive and
false negative rates. In the absence of selection, drift can generate
spurious outliers with strong patterns of reduced or increased
introgression that could be taken as evidence of selection. Conversely,
drift can act to obscure the effects of selection, making it even less
likely that outlier loci represent robust candidates.
Our comparisons show that the drift causes high variance inFST regardless of the type of selection acting
upon hybrids, creating patterns of differentiation in neutral genomic
regions similar to genomic regions under selection. Genomic
differentiation in turn increases the rate of false positives when
scanning the genome for signatures of selection. i This differentiation
increases even further as the hybrid zone ages. Migration rates plays
also an important role but, ultimately, having large deme sizes
substantially decrease genomic differentiation and produce the least
false positives and negatives.
Hybrid zones with small deme sizes and low migration rates show high
false positive rates, suggesting limitations for cline analyses. Similar
to the results by Polechova and Barton (2011) and Sciuchetti et al.
(2018) we found that drift increases risk of fixation within a deme, and
shifts the centers away from the expected hybrid index cline (Figure 4).
Increasing the population size does decrease the effects of drift,
generating more concordant and coincident clines, and reducing
significantly the rate of false positives. However, this reduction in
FPR is not sufficient in old hybrid zones with low migration rate, where
even with large population sizes, drift can maintain high false positive
rates.
Small deme sizes can also distort our assessment of how selection is
acting, even in cases in which we can correctly identify regions under
selection. Hybrid zones with small deme sizes and restricted migration
have higher risk of distortion caused by drift, masking outliers as a
different type of selection (e.g. adaptive introgression instead of an
ecotone or selection against hybrids). Hybrid zones with underdominance
and BDM incompatibilities, show clines with drastically shifted centers.
Since heterozygous genotypes have lower fitness, drift can increase the
rate of fixation, and shift the center to either side of the gradient
(e.g. Supporting Information Figure S5 and S6; Supporting Information
Videos 3 and 4), which could also be misinterpreted as directional
selection. Additionally, in ~3% of our replicates with
small deme sizes, genomic regions under directional selection, instead
of having centers shifted towards the adaptive side of the geographical
gradient, showed centers shifted towards the opposite (the less fit)
side of the gradient.
Between false positives on neutral and false negatives on underdominant
and BDMI loci, spatial analyses of hybrid zones may be markedly
overestimating the importance of adaptive introgression, or directional
selection on an allele from one parent species. Putative adaptive
introgression has been invoked as evidence of asymmetric sexual
selection (Parsons et al. 1993) and of resilience to climate change
(Hamilton and Miller 2016). Our simulations indicate that apparent
adaptive introgression may instead reflect drift or hybrid
incompatibilities.
Our results suggest that geographic cline analysis is most appropriate
in recently formed hybrid zones with large populations of hybridizing
individuals. Knowing the demographic history of the model organisma priori can reduce the risk of misidentifying neutral regions as
candidate regions under selection in natural populations. Geographic
cline analysis is most helpful with validation of the same outliers
across multiple independent hybrid zones (Schumer et al. 2014; Dufresnes
et al. 2015; Janousek et al. 2015) or temporal analyses of selection on
outlier candidates (Taylor et al. 2014; Ryan et al. 2018).
Geographic cline analysis is one of a family of outlier detection
methods used in evolutionary and ecological genomics, for example
genomic clines (Gompert and Buerkle 2011), ancestry tract lengths
(Sedghifar et al. 2016), and ancestry junctions (Hvala et al. 2018)).
All of these methods are susceptible to increased error rates due to
drift, in these cases by generating excess ancestry, increasing ancestry
tract sizes, and reducing junction densities. Further, we should expect
to see similar patterns when it comes to phenotypic introgression,
particularly when interspecific differences are associated with a small
number of loci.
Asymmetric introgression of traits or candidate genes is often taken as
prima facie evidence of selection. Our results suggest that drift can
generate both false positives and negatives that look like adaptive
introgression, rather than neutral movement of alleles or selection
against recombinant genotypes. Hybridization and introgression are
frequently favored by selection, and many taxa show pervasive evidence
of historical hybridization in the genome, but spatial analysis of
hybrid zones may lead us to overestimate the importance of adaptive
processes therein. Careful consideration of hybrid zone demographics,
combined with independent replication or validation of outliers, can
make geographic cline analysis a powerful framework for studying
genome-wide patterns of introgression.