# Detectability of varied hybridization scenarios using genome-scale hybrid
detection methods

Marianne Bjørner
Department of Computer Sciences
Wisconsin Institute for Discovery
University of Wisconsin-Madison
Madison, WI

&Erin K. Molloy
Department of Computer Sciences
University of Maryland
College Park, MD
&Colin N. Dewey
Department of Biostatistics
and Medical Informatics
University of Wisconsin-Madison
Madison, WI
&Claudia Solís-Lemus
Department of Plant Pathology
Wisconsin Institute for Discovery
University of Wisconsin-Madison
Madison, WI
Corresponding author<EMAIL_ADDRESS>
###### Abstract

Hybridization events complicate the accurate reconstruction of phylogenies, as
they lead to patterns of genetic heritability that are unexpected under
traditional, bifurcating models of species trees. This phenomenon has led to
the development of methods to infer these varied hybridization events, both
methods that reconstruct networks directly, as well as summary methods that
predict individual hybridization events from a subset of taxa. However, a lack
of empirical comparisons between methods – especially those pertaining to
large networks with varied hybridization scenarios – hinders their practical
use. Here, we provide a comprehensive review of popular summary methods: TICR,
MSCquartets, HyDe, Patterson’s D-Statistic (ABBA-BABA), D3, and Dp. TICR and
MSCquartets are based on quartet concordance factors gathered from gene tree
topologies and Patterson’s D-Statistic, D3, and Dp use site pattern
frequencies to identify hybridization events. We then use simulated data to
address questions of method accuracy and ideal use scenarios by testing
methods against complex networks which depict gene flow events that differ in
depth (timing), quantity (single vs. multiple, overlapping hybridizations),
and rate of gene flow ($\gamma$). We find that deeper or multiple
hybridization events may introduce noise and weaken the signal of
hybridization, leading to higher relative false negative rates across all
methods. Despite some forms of hybridization eluding quartet-based detection
methods, MSCquartets displays high precision in most scenarios. While HyDe
results in high false negative rates when tested on hybridizations involving
ghost lineages, HyDe is the only method to be able to separate hybrid vs
parent signal in hybridization events. Lastly, we test the methods on a
dataset of ultraconserved elements from the bee subfamily Nomiinae, finding
the possibility of hybridization events between clades which correspond to
regions of poor support in the species tree estimated in the original study.

##  1 Introduction

Phylogenetics studies the evolutionary history between organisms. In many
popular phylogenetic inference models, these relationships are assumed to be
best represented as a binary tree, where each child node arises from only one
direct parent [31, 7, 22]. However, a binary tree model ignores the
possibility of a reticulation or gene flow events. Gene flow occurs when
members of one population reproduce or otherwise exchange genetic information
with another population, forming admixed populations or hybrid species [4].
These reticulation events transform bifurcating phylogenetic trees into
network structures, wherein taxa affected have more than one parental lineage
[20]. Non-tree-like evolution is common across the tree of life, found in
groups such as insects [33], plants [11] and mammals [25].

The study of gene flow events in the tree of life have been aided by recent
advances in sequencing technology, granting evolutionary researchers access to
genome-scale information. This abundance of information can be used to infer
reticulation events such as introgression, hybrid speciation, and horizontal
gene transfer. Each mechanism for gene flow leaves behind various traces in a
population’s genetic information, and may be identified through hybridization
detection methods that leverage gene trees or sequence information [12].

Many pre-existing methods to infer phylogenies are based on a binary tree,
where they do not account for reticulations between taxa. While this
assumption limits the search space to only trees, it might not be a safe
assumption, especially for populations where gene flow is common or expected.
Methods to infer phylogenetic networks [29, 40, 35] are becoming increasingly
popular for their ability to overcome the strictly bifurcating assumption.
However, these methods are very computationally expensive and are thus not
applicable in real-life applications with increasingly large datasets.
Furthermore, most network methods require specification of the number of
expected hybrid events, as any increase in number of reticulations
artificially increases the likelihood of the network [17].

Alternatively, to find evidence for individual hybridization events, summary
methods analyze subsets of triples or quartets of taxa – which is an
intrinsically more scalable endeavor than the search in network space –
without any predetermination of the total of number of hybrids in the
phylogeny [12]. Despite these advantages, summary methods still require
comparisons to each other in order to address questions of method accuracy.
The hybrid detection methods compared in this simulation study (Table 1) are
TICR [32], MSCquartets [18, 19], HyDe [15], Patterson’s D-Statistic [24], also
known as the ABBA-BABA test, as well as methods derived from the D-Statistic
such as Dp [10] and D3 [9]. While TICR tests for how well a binary population
tree fits the data, the D-Statistics, MSCquartets and HyDe are used to
identify specific hybrid relationships within either quartets or triples of
taxa. Here, we address questions of method accuracy and ideal use scenarios by
testing these five summary methods against complex networks which depict gene
flow events that differ in depth (timing), quantity (number of hybridization
events: single vs. multiple), and proportion of genes transferred through the
hybridization event (inheritance probability $\gamma$). We note that we are
not treating hybridization events as continuous flow of genes over a period of
time. Instead, we represent a hybridization event by a single arrow (see
Figure 3). The hybridization scenarios depicted also differ in terms of time
consistency, a characteristic of reconstructed networks that, when violated,
arises in the presence of incomplete sampling or ghost lineages [20, 34, 23].
Finally, we use hybridization detection methods to analyze published empirical
data from the bee subfamily Nomiinae, a dataset complete with sequences, and
proposed gene and species trees.

##  2 Materials and Methods

###  2.1 Methods based on gene trees

TICR [32, 8] and MSCquartets [2, 18] both rely on frequencies of quartet gene
tree topologies to conduct tests into how well an expected tree fits the data
and whether there is evidence for the presence of hybridizations. These tests
are based on the multispecies coalescent model (MSC) which defines expected
distributions of gene tree topologies under incomplete lineage sorting (ILS)
[1]. ILS, also called deep coalescence, occurs when individual gene histories
fail to coalesce at the same time as their given species history [16]. An
example is shown in Figure 1.

Figure 1: Concordant (left) and discordant (right) gene and species histories.
Solid black lines denote species history, and dotted lines denote an
individual gene history. On the right, there is an example of incomplete
lineage sorting (ILS) when the gene history fails to coalesce with the
ancestral speciation event. The tree on the left corresponds to the
bipartition $12|34$ and is denoted a "major tree" as it agrees with the
species tree. Under the MSC model, this tree is expected to have a frequency
of 76% in a given sample of gene trees given the internal branch of $1.0$. The
tree on the right corresponds to the bipartition $14|23$ and is denoted a
"minor tree" as it is discordant to the species tree. Under the MSC model, it
has an expected frequency of 12%.

Samples of gene trees are comprised of topologies concordant and discordant
with the species tree. Together, these topologies form quartet count
concordance factors (qcCFs), where the CFs of concordant topologies are called
major CFs, as they align with the major tree, or species tree, and the CFs of
discordant topologies are called minor CFs. Minor CFs are expected in equal
frequency under ILS, and deviation from this expectation violates an ILS-only
model of genetic inheritance. These qcCFs are gathered calculating the
frequency of each of three possible resolved quartet trees across a set of
gene trees [2] (Figure 2).

####  2.1.1 MSCquartets

MSCquartets [2] is a package implemented in R [27]. From gene trees, it
tabulates counts for each set of four taxa as they form four different tree
structures (Figure 2): $12|34,13|24,14|23$, and $1234$ representing the
unresolved star tree [18].

Figure 2: Four possible quartet relationships, used by MSCquartets [2] and
TICR [32] depicted from left to right: $12|34$, $13|24$, $14|23$, and star
tree.

Unresolved (star) trees are removed, and the counts of each quartet topology
(qcCFs) are compared to expected values under a distribution defined in [18]
derived from the multispecies coalescent model to test the hypothesis of
whether a specific 4-taxon subset follows a tree-like pattern (in agreement
with the MSC model) or not. Each hypothesis test then produces a p-value. The
test builds on the expectation that minor CFs, or the concordance factors that
are discordant with a species tree, should be equal in probability or sample
proportion under the ILS-only model. When this is violated, it leads to low
support for a tree-like species history, subject to ILS-only, between the four
taxa and thus, there is evidence of hybridization in the quartet. In other
words, significant results fail to support an ILS-only model of evolution.
When the distribution of qcCFs between four taxa fall out of expectation, a
hybridization event becomes a possible explanation for the imbalanced
relationship between quartets of taxa.

In terms of efficiency, as this is a quartet method, the number of gene trees
and number of taxa are the two primary factors that influence the speed of the
method. Let $n$ be the number of taxa and $g$ is the number of gene trees, the
function scales with respect to $g\frac{n!}{4!(n-4)!}$, or $O(gn^{4})$. Thus,
if MSQquartets takes 1 second to analyze 300 gene trees from 10 taxa, it could
take roughly 5 hours to analyze 300 gene trees from 100 taxa.

####  2.1.2 TICR

In its original implementation, TICR – Tree Incongruence Checking in R – [32]
was used as part of a pipeline that begins with an alignment, calculates gene
tree topologies, calculates concordance factors from gene trees, and finally
calculates a population tree from concordance factors using the software
Quartet-Max-Cut [32, 28]. However, the only inputs required by TICR are qcCFs,
as observed from gene trees, and expected CFs as calculated from a
hypothesized population tree. Recently, this method was extended to test
goodness-of-fit on a given population network, rather than population tree
[8].

Here, we use the TICR version implemented in the Julia package called
QuartetNetworkGoodnessFit.jl [8] though we restrict our tests to the case of
population trees, not networks. Given branch lengths between taxa on a
population tree, as measured in coalescent time ($t$), the expected
probabilities of observing quartet relationships can be directly computed. For
example, in Figure 1, for an internal branch of $t=1.0$ coalescent units, the
probability of the major gene tree is given by $1-\frac{2}{3}e^{-t}=0.76$.

From expected and observed proportions, TICR is able to compute a $\chi^{2}$
goodness-of-fit test statistics that evaluates the fit of the observed quartet
qcCFs to the expected CFs under the ILS-only model. TICR can also be used to
test for the likelihood of panmixia, or a star tree which occurs when all taxa
arise from the same common ancestor and diverge at the same time. TICR uses
the p-values of the individual tests to form an overall test that inspects
whether the distribution of observed qcCFs falls within the expected qcCFs of
the input tree or network. This overall test indicates whether the proposed
population tree fits the observed quartet concordance factors. While TICR does
not directly test for the presence or absence of specific hybridizations, by
failing to reject a specific tree model, it provides lack of evidence for
hybridization events.

In terms of efficiency, the majority of the computational complexity occurs in
preprocessing the data to produce the input to TICR: aligning sequences,
estimating gene trees, and calculating concordance factors. With these
precalculated, for a set of $n$ taxa, there are $\frac{n!}{4!(n-4)!}$ quartets
possible for comparison which results in a time complexity of $O(n^{4})$.

###  2.2 Methods based on sequences

HyDe [15], Patterson’s D-Statistic [24], D3 [9], and Dp [10] are all methods
that use site pattern frequencies or pairwise differences to test the null
hypothesis of tree-like evolutionary patterns (ILS-only). This eliminates the
need for the estimation of gene trees, such that sequences can be used
directly as input. While HyDe, Patterson’s D-Statistic, and Dp use rooted
triples with an outgroup, D3 uses an rooted triple without an outgroup. This
is advantageous as a poorly chosen, or distant outgroup can result in
inclusion of ghost hybridizations and lead to false interpretations [34, 9].
Ghost hybridizations are defined as hybridization events when one (or both) of
the parent lineages that provide genetic material to the hybrid node are
either extinct or unsampled. Note that the Patterson’s D-statistic, D3, and Dp
are intended as tests for introgression between two species, but do not
indicate directionality. Additionally, they require that the species
relationship between the triple and its outgroup (if any) is known.

####  2.2.1 HyDe

Distinct from other methods which evaluate for the overall presence of
hybridization, HyDe [5] identifies a singular parent-hybrid relationship
between a triple, given its outgroup. It can also be used to estimate a mixing
parameter $\gamma$, depicting the proportion of genetic material contributed
by each parental lineage.

HyDe is based on phylogenetic invariants, or a function of site pattern
probabilities, which evaluate to zero when consistent with given associated
tree models [1, 15]. The linear invariants ($f_{1}$ and $f_{2}$) are
themselves representative of mixing parameters $\gamma$ and $1-\gamma$,
respectively,

$\displaystyle f_{1}$
$\displaystyle=p_{iijj|(S\gamma,\tau)}-p_{ijij|(S\gamma,\tau)}$ $\displaystyle
f_{2}$ $\displaystyle=p_{ijji|(S\gamma,\tau)}-p_{ijij|(S\gamma,\tau)}$

where $p_{iijj|(S\gamma,\tau)}$ represents the probability for the site
pattern $iijj$ under the species triple tree $S\gamma$ with $\gamma$ as
inheritance probability and $\tau$ as divergence times. As a result, HyDe can
also be used to estimate the mixing parameter $\gamma$ between two taxa that
are putative parent lineages of a proposed hybrid as the ratio of
$\frac{f_{1}}{f_{2}}$ represents $\frac{\gamma}{1-\gamma}$. When there is no
admixture, and $\gamma$ is 0, the ratio is expected to be zero.

Site pattern probabilities observed in the sample are used to form estimates
of $f_{1}$ and $f_{2}$, along with means and variance. These are then
rearranged with the Geary-Hinkley transformation to form the Hils statistic:

$H:=\frac{\hat{f}_{2}(\frac{\hat{f}_{1}}{\hat{f}_{2}}-\frac{\mu_{f_{1}}}{\mu_{f_{2}}})}{\sqrt{{\hat{\sigma}^{2}}_{f_{2}}(\frac{\hat{f}_{1}}{\hat{f}_{2}})^{2}-2\hat{\sigma}_{f1,f2}\frac{\hat{f}_{1}}{\hat{f}_{2}}+\hat{\sigma}_{f1}^{2}}}$

which follows the normal distribution with $N(0,1)$ for large $n$ under the
null hypothesis of ILS-only, or no hybrid speciation [15]. This allows direct
interpretation of HyDe test results without the need for resampling by
bootstrapping. However, significance levels should be adjusted with a
Bonferroni correction due to multiple hypothesis testing, as HyDe considers
all possible combinations of sets of three taxa, where one hybrid taxon is
tested for every two distinct parent taxa.

As a C-backed python package, HyDe is designed for use with multiple
individuals per taxon, from a phylogeny where the outgroup is specified, and
aligned sequence information is provided in PHYLIP format. Power increases
with increasing sequence length, with a recommendation that sequence length is
at minimum 50 kbp [5, 14]. In addition, HyDe outputs counts of site pattern
observations, namely AABB, ABBA, AABC,… that can be used for calculations of
other statistics based on site pattern frequencies such as the Patterson’s
D-Statistics (ABBA-BABA) [24], Dp [10], and can be rearranged for use in
pairwise distance metrics, such as D3 [9].

In terms of efficiency, HyDe "operates in less than 30 seconds" according to
[15] for up 20 taxa and 100,000bp. In practice, the method depends upon an
already aligned set of of sequences for SNP comparison, so computational
workload is based on the calculation of site pattern frequencies, which is a
function of the number of taxa to be compared, as well as the length of the
input sequences. For each set of three taxa $(t_{1},t_{2},t_{3})$, there are
three possible hybrid relationships:

  1. 1.

$t_{1}$ and $t_{2}$ are parents of hybrid $t_{3}$

  2. 2.

$t_{1}$ and $t_{3}$ are parents of hybrid $t_{2}$

  3. 3.

$t_{2}$ and $t_{3}$ are parents of hybrid $t_{1}$

Sites from each taxa $t_{1}$, $t_{2}$, and $t_{3}$ are then compared to sites
from an outgroup to calculate the Hils statistic from $f1$ and $f2$. For a set
of $n+1$ taxa, (where $n$ is all taxa not including the outgroup), where each
has a sequence length of size $L$, HyDe performs $3\frac{n!}{3!(n-3)!}$ tests,
where each test involves comparison of $L$ sites, across each of the three
taxa for hybrid testing, plus the outgroup resulting in a complexity of
$O(Ln^{3})$.

####  2.2.2 Patterson’s D-Statistic (ABBA-BABA)

Patterson’s D-Statistic [24], much like HyDe, involves the calculation of site
pattern frequencies from SNPs and scales in a similar manner with respect to
sequence length and number of taxa to compare. Patterson’s D-statistic is
defined as $D=\frac{ABBA-BABA}{ABBA+BABA}$ where $A$ and $B$ are either
distinct alleles or sites. Each of the four positions in these nucleotide
sequences belong to one of four taxa: three of these taxa are compared to each
other for hybrid identification, and the fourth is used as an outgroup. Under
the null hypothesis of no hybridization, the D statistic is expected to be
zero, as the frequency of “ABBA" and “BABA" occurring should be equal.
Significant deviation from 0 implies gene flow events (Figure 3). Note that
the Patterson’s D-statistic does not indicate directionality of this
introgression event.

Figure 3: Center: Population tree of the form (((T1,T2),T3),O). Under each
taxa, ABBA, BABA, and BBAA are three possible site patterns, with positions in
the site pattern sequence corresponding to particular taxa. Under ILS, the
proportion of sites where T3 and T1 match should be equal to the proportion of
sites where T2 and T3 match, given that T3 is equidistant on the species tree
from T1 and T2. Left: An excess of ABBA, corresponds to a gene flow event
between T2 and T3. Right: An excess of BABA, corresponds to a gene flow event
between T1 and T3.

In terms of efficiency, testing on all permutations of three taxa, with a
specified outgroup, computing time of Patterson’s D-Statistic on an alignment
of length $L$, containing $n$ taxa scales with respect to
$4L\times\frac{n!}{(n-3)!}$. However, Patterson’s D-statistic is symmetric,
such that the resulting D-statistic of the topology ((a,b),c) is equal to the
negative D-statistic given the topology ((b,a),c). Further reduction of the
number of tests occurs when provided the correct topology of the underlying
species tree. As a result, only one test is necessary per triple $O(Ln^{3})$.

####  2.2.3 D3

Motivated by the original Patterson’s D-Statistic, D3 is a "faster, but less
powerful" alternative according to [9]. This test is SNP-focused and based on
pairwise differences which likely has limitations when applied to larger
genomes and indels. In addition, the test can only distinguish hybridization
between non-sister lineages [9, 12]. D3 uses only three taxa, and omits the
outgroup used in Patterson’s D-statistic. It uses pairwise distances instead
of site patterns frequencies, where the pairwise distances between taxa
$t_{3},t_{2}$ and between taxa $t_{3},t_{1}$ are expected to be equal, given a
known species tree relationship, where $t_{1}$ and $t_{2}$ are sisters (such
as in $((t_{1},t_{2}),t_{3})$).

D3 can be calculated as a ratio of pairwise differences between three
sequences $D3=\frac{d_{BC}-d_{AC}}{d_{BC}+d_{AC}}$ where $d_{ij}$ corresponds
to the distance between taxa $i$ and $j$. Here, significant deviation of D3
from 0 may imply gene flow between taxa C and B, in the case of a positive
result, and between A and C, in the case of a negative result. This operates
much like the original D-statistic to test for the presence, but not
directionality, of gene flow events.

In terms of efficiency, unlike the Patterson’s D-Statistic and Dp, D3 does not
include comparison to an outgroup. As a result, the efficiency of the method
is increased, as there is 25% less sequence information to analyze and compare
to, as there are now three sequences instead of four. As with the Patterson’s
D-Statistic, prior knowledge of the species topology must be known, and due to
symmetry, only one test per combination of three taxa is necessary
$O(Ln^{3})$.

####  2.2.4 Dp

Dp adds the site pattern frequency BBAA to the denominator of the original
Patterson’s D-Statistic in order to estimate the net proportion of the genome
resulting from introgression [10]. This feature provides comparability with
HyDe’s computation of $\gamma$. The denominator in Dp accounts for the total
number of variable sites: $D_{p}=\frac{|ABBA-BABA|}{BBAA+ABBA+BABA}$.

In terms of efficiency, as with other forms of the D-test, Dp tests each
combination of $n$ taxa given an alignment of length $L$ and a specified
topology. Its complexity is equivalent to that of Patterson’s D-Statistic:
$O(Ln^{3})$.

Method |  Complexity |  Input |  Test |  Output
---|---|---|---|---
TICR |  $O(n^{4})$ |  Population tree and observed qcCFs from gene trees |  If observed qcCFs align with expected CFs under ILS-only MSC |  p-values from individual quartet tests and overall p-value reporting goodness-of-fit to population tree
MSCquartets |  $O(gn^{4})$ |  Gene trees |  If minor qcCFs align with expected CFs under ILS-only MSC |  p-values from individual quartet tests
HyDe |  $O(Ln^{3})$ |  Sequences from taxa and outgroup |  Hils statistic evaluates triples for site pattern frequencies against phylogenetic invariants |  List of all triples with 1) proposed hybrid speciation (parental and hybrid lineages), 2) a estimated $\gamma$, and 3) p-value
Patterson’s D |  $O(Ln^{3})$ |  Sequences from taxa and outgroup |  Excess in site pattern frequencies of $ABBA$ and $BABA$ when expected to be equal |  A D-statistic per triple where significant deviation from 0 indicates gene flow
D3 |  $O(Ln^{3})$ |  Sequences |  If triples reject ILS-only assumption (measured as deviation of D3 from 0) |  Distribution of D3 values of triples mapped to z-distribution
Dp |  $O(Ln^{3})$ |  Sequences from taxa and outgroup |  Excess in site pattern frequencies of $ABBA$ and $BABA$ when expected to be equal |  A Dp value per triple where significant deviation from 0 indicates presence and extent of gene flow
Table 1: Summary Table of Hybrid Detection Methods

###  2.3 Simulations

All methods were tested on the same proposed networks and compared in their
ability to test for the presence of hybridization events within either the
triples or quartets of taxa used as input. Networks used for simulation are
illustrated in Figures 4, 5, and 6. These vary in number of reticulations,
number of taxa, depth of reticulations and their mixing parameter $\gamma$
which denotes how much ancestral DNA is passed from the minor hybrid edge to
the hybrid node. Size ranges from 4–25 taxa, where the number of reticulation
events for networks with 10, 15, and 25 taxa are 20% of the number of taxa.
Each reticulation event can have singular or multiple affected taxa downstream
of the hybridization event. We name the networks based on the number of taxa
($n$) and number of hybridizations ($h$). For example, the first network in
Figure 4 is denoted n4h1 as it has four taxa ($n=4$) and one hybridization
($h=1$). Six out of the twelve networks are replicated from earlier studies:
1) the four networks (n4h1, n4h1${}_{\textrm{introgression}}$, n8h3 and n5h2)
in Figure 4 were used in [14] and 2) the networks n10h2 and n15h3 (in Figures
5 and 6) were used in [29].

We separate networks labeled n10h2 (with two hybridization events) and n15h3
(with three hybridization events) into singular reticulation events
(n10h1${}_{\textrm{deep}}$, n10h1${}_{\textrm{shallow}}$,
n15h1${}_{\textrm{deep}}$, n15h1${}_{\textrm{intermediate}}$,
n15h1${}_{\textrm{shallow}}$ in Figures 5 and 6). These single-hybridization
networks allow us to measure the ability of methods to detect singular
hybridization events without the possible influence of overlapping
hybridizations, and to compare the performance of the methods on shallow vs
deep hybridizations as it has been reported that deep hybridizations are more
difficult to detect [12]. Note that here we use the term "overlapping
hybridizations" not to refer to hybridizations that share edges (e.g. level-2
networks), but to hybridizations that affect the same set of taxa downstream.
Finally, to represent how well methods perform at a larger scale, both in
terms of efficiency and accuracy, we evaluate their performance on a larger
network labeled n25h5, shown in Figure 5.

We simulated gene trees under each network with the software ms [13], with a
single individual per taxon. Note that this approach reduces the power of
HyDe, when compared to simulations with multiple individuals per population
[14, 15]. The software ms allows hybridization events to be modeled with -es t
i p and -ej t i j events which correspond to population admixture and
population splitting, respectively [13]. With this approach, we circumvent the
alternative of decomposing each network into $2^{r}$ trees as in [14], where
$r$ is the number of reticulations, and sampling a proportion of each tree to
represent the mixing parameter $\gamma$. Directly modeling reticulation events
from a network also offers the advantage that gene trees produced follow their
underlying network structure, which has a different probability density than
the combination of individual trees, especially under complex reticulation
events [37]. The decomposition to individual trees may produce reticulate gene
tree patterns that are artificially clearer, and may not be as accurately
representative of the timing of natural gene flow events.

For TICR and MSCquartets, we simulate gene trees for
$g\in\\{30,100,300,1000,3000\\}$ to be used directly as input. For methods
which require sequences, ms generates $g$ unlinked gene trees where
$g\in\\{30,100,1000,3000,10000\\}$ and from each gene tree, seqgen [26] is
used to generate sequences with 100 base pairs generated per gene tree,
similar to the approach used in [14] and then concatenated to form sequences
of total length $L$. The seqgen parameters used to generate these base pairs
are -mHKY -s0.036 -f0.300414,0.191363,0.196748,0.311475 -n1 -l100. Full
simulation details including ms commands and newick structures of these
networks can be found on the GitHub repository
https://github.com/mbjorner/hybrid-detection-comparison. 30 trials were
simulated for each combination of network and gene tree number or sequence
length. A pipeline of this simulation is described in Figure 7.

Note that the D-derived tests rely on pre-specification of topology and thus,
we can expect increased false positives when testing on inputs of
$(((t_{1},t_{3}),t_{2}),O)$ when $(((t_{1},t_{2}),t_{3}),O)$ is the true
topology. In addition, HyDe relies on the existence of concurrent parental
lineages to test for hybrid speciation. Where only one parental lineage is
sampled, as, for example, in network n10h2 in Figure 5, we investigate the
influence of introgression events from these ’ghost’ lineages on HyDe’s
output. Last, TICR requires an input population tree to be used for the
expected qcCFs. Deviations from the expected qcCFs could indicate deviations
from the ILS-only model, but also, it could indicate that the wrong population
tree was used for comparison. In our simulation studies, we use the known
major tree as the input population tree for TICR so that any significant TICR
results are interpreted with the possibility of hybridization.

(a) n4h1

(b) n4h1${}_{\textrm{introgression}}$

(c) n5h2

(d) n8h3

Figure 4: Networks obtained from [14]. All branches are labeled in coalescent
units. (a) n4h1, a four taxon network where taxon 2 arises from ancestral
lineages to taxon 1 and taxon 3 (denoted ’parents’ in HyDe). This network
depicts a hybrid speciation scenario and was simulated with various values for
$\gamma\in\\{0.0,0.1,0.2,0.3,0.4,0.5\\}$ (inheritance probability on the minor
hybrid edge in blue). (b) n4h1${}_{\textrm{introgression}}$ depicts an
introgression event where $\gamma=0.5$ with a single hybrid labeled 2 that has
introgressed with taxon 1. (c) n5h2 depicts overlapping hybridizations where
one reticulation (ancestral to the (2,3) clade) forms the parental lineage of
another reticulation event (taxon 2). (d) n8h3 depicts three singular
hybridization events, each affecting single taxa, two of identical depths and
one deeper event.

(a) n102h

(b) n10h1${}_{\textrm{shallow}}$

(c) n10h1${}_{\textrm{deep}}$

(d) n25h5

Figure 5: (a) Network n10h2 replicated from [29] with ten taxa and two
hybridization events, its singular hybridizations, (b)
n10h1${}_{\textrm{shallow}}$ and (c) n10h1${}_{\textrm{deep}}$, and the
largest network tested, (d) n25h5. All branches are labeled in coalescent
units. Network n10h1${}_{\textrm{shallow}}$ (b) is the n10h2 network
containing only the shallowest hybridization ($\gamma=0.2$), which impacts one
taxon, labeled 9. Network n10h1${}_{\textrm{deep}}$ (c) is a derivation the
n10h2 network containing only the deepest hybridization event ($\gamma=0.3$),
which impacts taxa 3-9. Network n25h5 (d) contains 25 taxa and five
hybridization events, three of which impact only a singular downstream taxon
($\gamma=0.449,0.395,0.369$), with two that are deeper in the tree
($\gamma=0.024,0.334$), creating overlapping hybridization events.

(a) n15h3

(b) n15h1${}_{\textrm{shallow}}$

(c) n15h1${}_{\textrm{intermediate}}$

(d) n15h1${}_{\textrm{deep}}$

Figure 6: (a) Network n15h3 replicated from [29] with 15 taxa and 3
hybridizations events, and its singular hybridizations, (b)
n15h1${}_{\textrm{shallow}}$, (c) n15h1${}_{\textrm{intermediate}}$, and (d)
n15h1${}_{\textrm{deep}}$. All branches are labeled in coalescent units.
Network n15h1${}_{\textrm{shallow}}$ (b) is the n15h3 network containing only
the shallowest hybridization ($\gamma=0.2$), which impacts only taxon 10.
Network n15h1${}_{\textrm{intermediate}}$ (c) is the n15h3 network with only
the intermediate-depth hybridization ($\gamma=0.2$) impacting the clade (2,3).
Network n15h1${}_{\textrm{deep}}$ (d) is the n15h3 network with only the
deepest hybridization event ($\gamma=0.3$). impacting the clade (2,3,4,5,6).
Figure 7: Simulation Pipeline. The software ms is used to simulate $g$ gene
trees from a known network where $g\in\\{30,100,300,1000,3000\\}$. These gene
trees are used as input for the hybrid detection methods TICR and MSCquartets.
Additionally, sequences of length $L$ are generated using seqgen from $g$ gene
trees where $g\in\\{30,100,300,1000,3000,10000\\}$. Each gene tree is used to
generate sequences with 100 base pairs generated per gene tree as in [14].
These sequences are used as input for HyDe, D3, Patterson’s D-Statistic and
Dp. This process was repeated for each network structure in Figures 4, 5, and
6 with 30 replicates each.

Each method is also evaluated in terms of computing time, as measured in CPU
time in seconds, given their gene tree or sequence data inputs, for the
purpose of predicting how well each summary method accommodates the addition
of sampled taxa. As genetic sequence information has become more widely
available, so too have the datasets that biologists use to construct
phylogenies and infer these reticulation events. In practice, often tens or
hundreds of taxa are compared [6, 33]. Since the Patterson’s D-Statistic, D3,
and Dp were computed from the output of HyDe, which describes all possible
site pattern frequencies, timing was omitted for D-statistic related tests.

###  2.4 Hybridizations in the bee subfamily Nomiinae

To demonstrate the use of these hybrid detection methods on real data, we
compare method performance on a dataset of ultraconserved elements (UCEs) from
the bee subfamily Nomiinae. This data originates from a paper investigating
the impacts of gene tree estimation error on species tree reconstruction [6],
and was used to demonstrate improved tree reconstruction with weighted [38], a
new version of ASTRAL that weights quartets based on their uncertainty (branch
support) and terminal branch lengths in input gene trees. The dataset is
available for download on
https://datadryad.org/stash/dataset/doi:10.5061/dryad.z08kprrb6.

This dataset contains sequences and gene trees of up to 852 UCEs, for a total
concatenated sequence length of 576,041 base pairs for each of 32 taxa. In the
original study [6], gene trees were estimated using six different methods, (1)
IQ-Tree2 with the GTR-G substitution model, (2) IQ-Tree2 with the substitution
model chosen by ModelFinder, (3) MrBayes with the GTR-G substitution model,
(4) MrBayes with reversible jump MCMC, (5) PhyloBayes, and (6) RAxML. The
original investigation found a consensus tree using PhyloBayes on concatenated
UCEs.

Here, we use each of the proposed sets of gene trees created using the six
different methods, as input for MSCquartets, and apply a Bonferroni correction
to evaluate significant quartets which may contain hybridization. We use the
proposed species tree and gene trees in combination for TICR, for which we
interpret a poor fit of gene trees’ quartet concordance factors to either
possibility of incorrect species tree, presence of hybridization, or a
combination of the two. Next, we use the original UCE sequences and
concatenate them to run HyDe, using Lasioglossum albipes as the outgroup, as
indicated by the consensus tree constructed with PhyloBayes [6], and a
Bonferroni correction for significance. As the original study included two
outgroups, Lasioglossum albipes and Dufourea novaeangilae, we removed Dufourea
novaeangilae from all gene trees and sequences prior to running hybrid
detection methods because these methods require only one outgroup.

##  3 Results

###  3.1 Simulations

Figure 8 shows the proportion of times that TICR correctly rejects the major
tree, and thus, detects the presence of hybridizations under the different
networks under study. We highlight that TICR accurately detects hybridizations
for the case of single shallow hybridizations, n10h1${}_{\textrm{shallow}}$
and n15h1${}_{\textrm{shallow}}$. However, TICR struggles to detect deeper
hybridizations as in n10h1${}_{\textrm{deep}}$,
n15h1${}_{\textrm{intermediate}}$, and n15h1${}_{\textrm{deep}}$ or multiple
hybridizations in the same network as in n5h2, n8h3, n10h2, n15h3, and n25h5.
Cases not shown in the figure correspond to all bars equal to zero. Being a
quartet-based method, TICR has no power to detect hybridizations on networks
with four taxa (as n4h1 or n4h1${}_{\textrm{introgression}}$) and those
results are thus not included in the figure. We can conclude that TICR is a
powerful method to detect single shallow hybridization events, even if they
involve ghost lineages (as n15h1${}_{\textrm{shallow}}$). This method,
however, is not powerful enough to detect deep or multiple hybridizations.

Figure 8: Proportion of replicates in which TICR correctly rejects the major
tree at a level of $0.05$. X axis corresponds to the different types of
networks and the Y axis corresponds to the proportion. TICR rejects the null
hypothesis of the major tree (equivalent to the input network sans
hybridization events) for networks with single, shallow hybridizations. TICR
fails to detect deeper hybridizations such as n10h1${}_{\textrm{deep}}$,
n15h1${}_{\textrm{intermediate}}$, and n15h1${}_{\textrm{deep}}$ which are
omitted in the figure as all bars are collapsed to zero. For networks with
multiple hybridization events (n5h2, n8h3, n10h2, n15h3, and n25h5), TICR is
also unable to reliably detect non-treelike patterns (some omitted again
because all bars are collapsed to zero).

Figures 10, 11, 12, 13 and 14 show the false positive rate (orange), precision
(black) and false negative rate (red) for MSCquartets, HyDe, Patterson’s
D-Statistic, and D3 on all the networks under study (Figures 4, 5, and 6). The
results for Dp are omitted as they are similar to those for Patterson’s
D-Statistic. We describe each figure below and here we describe the
computation of the false positive/negative rates and precision (see also
Figure 9). Every triple (or quartet for MSCquartets) could have a hybrid or
not. For example, in the n4h1 network in Figure 4, the triple $\\{1,2,3\\}$
contains a hybrid (taxon $2$), but the triple $\\{1,3,4\\}$ does not contain
any hybrid. If the triple (quartet) contains a hybrid, and the method detects
it (pvalue$<\alpha$ for significance level $\alpha$), we consider this a true
positive (TP). If the triple (quartet) contains a hybrid, but the method does
not detect it (pvalue$>\alpha$), we consider this a false negative (FN). If
the triple (quartet) does not contain a hybrid, and the method finds no hybrid
(pvalue$>\alpha$), we consider this a true negative (TN). If the triple
(quartet) does not contain a hybrid, but the method detects a hybrid
(pvalue$<\alpha$), we consider this a false positive (FP). The False Positive
Rate (FNR) is computed as $FP/(FP+TN)$. The False Negative Rate (FNR) is
computed as $FN/(FN+TP)$. The precision is computed as $TP/(TP+FP)$. For HyDe,
an additional metric, Wrong Hybrid Rate (WHR) describes the rate at which
hybridization is detected but is falsely attributed to the incorrect hybrid
taxon. That is, if the triple contains a hybrid, and HyDe detects it, but
identifies the wrong taxon as the hybrid taxon, we consider this a wrong
hybrid (WH). For example, in the n4h1 network in Figure 4, the triple
$\\{1,2,3\\}$ contains a hybrid (taxon $2$). HyDe could test whether $1$ and
$3$ are the parents of hybrid taxon $2$ (correct hybrid), or whether $2$ and
$3$ are parents of hybrid taxon $1$ (wrong hybrid). If the latter test is
significant, then HyDe correctly identified that there is a hybrid
relationship among these taxa, but wrongly identified the hybrid taxon. We
define the Wrong Hybrid Rate (WHR) as $WHR=WH/(WH+TP+FP)$. We use
PhyloNetworks [30], a Julia package that allows for efficient manipulation of
phylogenetic networks to easily identify triples or quartets with hybrid
relationships in all networks under study.

Figure 9: Visual description of false positives, false negatives, wrong
hybrids (for HyDe only), true positives and true negatives.

Figure 10 shows the false positive rate (orange), precision (black) and false
negative rate (red) for MSCquartets, HyDe, Patterson’s D-Statistic, and D3 on
the n4h1 network with mixing parameters $\gamma\in 0,0.1,0.2,0.3,0.4,0.5$. For
HyDe, an additional metric, wrong hybrid rate (purple) describes the rate at
which hybridization is detected but is falsely attributed to the incorrect
hybrid taxon. This network is among the simplest cases with a single shallow
hybridization, and thus, all methods have a high precision and low false
positive/negative rates for as little as 100,000 sites or 300 gene trees. In
this figure, all test are Bonferroni-corrected at a level of significance
$\alpha=0.05/\textrm{number of tests}$, but we also show the uncorrected
version ($\alpha=0.05$) in Figure 18 in the Appendix.

Figure 10: False positive rate (orange), precision (black) and false negative
rate (red) for MSCquartets, HyDe, Patterson’s D-Statistic, and D3 on the n4h1
network with mixing parameters $\gamma\in 0,0.1,0.2,0.3,0.4,0.5$, where 0
depicts no gene flow, and 0.5 depicts equal contribution from both parental
lineages. All tests are Bonferroni-corrected at a level of significance
$\alpha=0.05/\textrm{number of tests}$ (see Figure 18 in the Appendix for
results on $\alpha=0.05$). For HyDe, an additional metric, Wrong Hybrid Rate
(purple) describes the rate at which hybridization is detected but is falsely
attributed to the incorrect hybrid taxon. In the absence of gene flow, HyDe
and D3 have the lowest proportion of false positives, but the difference is
not considerable when compared to that of Patterson’s D-statistic and
MSCquartets. With all other levels of admixture, as both admixture and sample
size increases, all tests display a decrease in false negative and false
positive rates. In cases where gene flow is present, D3 requires slightly
higher sequence lengths to eliminate false negatives when compared to
Patterson’s D-Statistic. However, given a sufficient sequence length, all
tests are able to reliably detect the presence of hybridization.

Figure 11 shows the false positive rate (orange), precision (black) and false
negative rate (red) for MSCquartets, HyDe, Patterson’s D-Statistic, and D3 on
the networks: n4h1${}_{\textrm{introgression}}$ (network with single shallow
introgression event), n5h2 (network with two overlapping hybridization events)
and n8h3 (network with three overlapping hybridization events). As in [14], an
overlapping hybridization event is defined as a hybridization where one taxa
is the parent of multiple hybridization events. For HyDe, an additional
metric, wrong hybrid rate (purple) describes the rate at which hybridization
is detected but is falsely attributed to the incorrect hybrid taxon. The
network n4h1${}_{\textrm{introgression}}$ displays an introgression event
which is easily detected by all methods (high precision and low false negative
rates). All methods also display no false positive rates on this network, as
all triples or quartets tested contain a hybrid relationship. For the case of
two hybridizations (n5h2), all methods display a high precision and low false
negative rate, except for HyDe which has a higher false negative rate than
others. False positive rate is low and comparable for all methods in this
network. Last, for three hybridizations (n8h3), all methods have high
precision and higher false negative rates. As more taxa become part of the
network, certain combinations contain hybrids that arise from ghost lineages,
which may not have strong signal to detect the hybridization events. In this
figure, all test are Bonferroni-corrected at a level of significance
$\alpha=0.05/\textrm{number of tests}$, but we also show the uncorrected
version ($\alpha=0.05$) in Figure 19 in the Appendix.

Figure 11: False positive rate (orange), precision (black) and false negative
rate (red) for MSCquartets, HyDe, Patterson’s D-Statistic, and D3 on the
networks: n4h1${}_{\textrm{introgression}}$ (network with single shallow
hybridization event), n5h2 (network with two overlapping hybridization events)
and n8h3 (network with three overlapping hybridization events). All tests are
Bonferroni-corrected at a level of significance $\alpha=0.05/\textrm{number of
tests}$ (see Figure 19 in the Appendix for results on $\alpha=0.05$). For
HyDe, an additional metric, Wrong Hybrid Rate (purple) describes the rate at
which hybridization is detected but is falsely attributed to the incorrect
hybrid taxon. For the introgression event (n4h1${}_{\textrm{introgression}}$),
MSCquartets, HyDe, Patterson’s D-Statistic, and D3 all behave similarly to the
hybrid speciation scenario in n4h1, with near-perfect recovery and
identification of hybrid scenarios. Low false negative rates are noted for
n5h2 across MSCquartets and D-related methods. For n5h2, HyDe correctly
recovers parent-hybrid relationships, but has a slightly elevated false
negative rate. Network n8h3 displays elevated false negative rates across all
methods. However, precision is high, and false negative rates between all
methods are comparable.

Figure 12 shows the false positive rate (orange), precision (black) and false
negative rate (red) for MSCquartets, HyDe, Patterson’s D-Statistic, and D3 on
the networks: n10h2 (network with two hybridization events),
n10h1${}_{\textrm{shallow}}$ (network with single shallow hybridization event)
and n10h1${}_{\textrm{deep}}$ (network with a single deep hybridization
event). For HyDe, an additional metric, wrong hybrid rate (purple) describes
the rate at which hybridization is detected but is falsely attributed to the
incorrect hybrid taxon. All methods report a higher false negative rate
compared to the simpler networks (n4h1 and n5h2), although precision continues
to be high for all methods, except for HyDe. HyDe’s lower precision is due to
the fact that some hybridizations involve ghost lineages (hybridizations when
one or both parental lineages contributing to the hybrid node are extinct or
unsampled) and HyDe cannot account for this scenario. False positive rate is
controlled in all methods. This combined with the higher negative rate allow
us to conclude that multiple overlapping hybridizations result in loss of
signal for hybridization, rather than contradicting signal pointing at wrong
hybrids. In addition, the false negative rate is high across all methods for
the single deep hybridization case (n10h1${}_{\textrm{deep}}$) which means
that it is not only multiple hybridizations that result in loss of signal,
hybridizations occurring in deeper parts of the tree have also lost signal to
be detected. We also note that unlike previous cases (Figures 10 and 11) where
HyDe’s wrong hybrid rate (purple line) and false positive rate (orange line)
were overlapped, for these cases, the wrong hybrid rate is much higher than
the false positive rate. This implies that HyDe is better able to identify
hybrid relationships for these networks, but not the correct hybrid taxon. In
this figure, all test are Bonferroni-corrected at a level of significance
$\alpha=0.05/\textrm{number of tests}$, but we also show the uncorrected
version ($\alpha=0.05$) in Figure 20 in the Appendix.

Figure 12: False positive rate (orange), precision (black) and false negative
rate (red) for MSCquartets, HyDe, Patterson’s D-Statistic, and D3 on the
networks: n10h2, n10h1${}_{\textrm{shallow}}$ (n10h2 with only the shallow
hybridization event), n10h1${}_{\textrm{deep}}$ (n10h2 with only the deeper
hybridization event). All tests are Bonferroni-corrected at a level of
significance $\alpha=0.05/\textrm{number of tests}$ (see Figure 20 in the
Appendix for results on $\alpha=0.05$). For HyDe, an additional metric, Wrong
Hybrid Rate (purple) describes the rate at which hybridization is detected but
is falsely attributed to the incorrect hybrid taxon. D3, Patterson’s
D-Statistic, and MSCquartets behave comparably across network structures, with
high precision, low false positive rate, but higher false negatives, which is
seen in the deepest hybridization (n10h1${}_{\textrm{deep}}$). A shallower
hybridization that creates a single hybrid taxon
(n10h1${}_{\textrm{shallow}}$) is readily detected by these tests, but still
with high false negative rates. HyDe has the highest false negative rate with
the n10h1${}_{\textrm{deep}}$ hybridization event, and the lowest with
n10h1${}_{\textrm{shallow}}$. While HyDe correctly detects hybridization in
triples, it often misattributes the hybridization to the wrong parent or
child. When both hybridization events are present, as in n10h2, all methods
result in a high number of false negatives with HyDe also having a higher
number of false positives in this case.

Figure 13 shows the false positive rate (orange), precision (black) and false
negative rate (red) for MSCquartets, HyDe, Patterson’s D-Statistic, and D3 on
the networks: n15h3 (network with three hybridization events),
n15h1${}_{\textrm{shallow}}$ (network with a single shallow hybridization
event), n15h1${}_{\textrm{intermediate}}$ (network with a single intermediate
hybridization event), and n15h1${}_{\textrm{deep}}$ (network with a single
deep hybridization event). For HyDe, an additional metric, wrong hybrid rate
(purple) describes the rate at which hybridization is detected but is falsely
attributed to the incorrect hybrid taxon. As already shown in the case of
n10h2 (Figure 12), all methods have a higher false negative rate, but
controlled false positive rate except for Patterson’s D-statistic with a high
false positive rate for the case of three hybridizations (n15h3). HyDe shows
lower precision compared to other methods which is due to the fact that the
hybridizations in n15h3 involve ghost lineages which HyDe cannot account for.
The fact that there are high false negative rates even for the single shallow
hybridization (n15h1${}_{\textrm{shallow}}$) across of methods could provide
some evidence that ghost lineages create challenges, not just for HyDe. In
this figure, all test are Bonferroni-corrected at a level of significance
$\alpha=0.05/\textrm{number of tests}$, but we also show the uncorrected
version ($\alpha=0.05$) in Figure 21 in the Appendix.

Figure 13: False positive rate (orange), precision (black) and false negative
rate (red) for MSCquartets, HyDe, Patterson’s D-Statistic, and D3 on the
networks: n15H3, n15h1${}_{\textrm{deep}}$ (n15h3 with only the deepest
hybridization event), n15h1${}_{\textrm{shallow}}$ (n15h3 with only the
shallowest hybridization event), and n15h1${}_{\textrm{intermediate}}$ (n15h3
with a hybridization event of intermediate depth). All tests are Bonferroni-
corrected at a level of significance $\alpha=0.05/\textrm{number of tests}$
(see Figure 21 in the Appendix for results on $\alpha=0.05$). For HyDe, an
additional metric, Wrong Hybrid Rate (purple) describes the rate at which
hybridization is detected but is falsely attributed to the incorrect hybrid
taxon. MSCquartets results in higher false negative rates given deeper
hybridization events, with a stable, high precision rate given more than 30
gene trees as input. All methods have a controlled false positive rate, except
Patterson’s D-Statistic with n15h3. HyDe has a similar or higher false
negative rate when compared to other tests, and again misidentifies the exact
hybrid-parent relationship in which there is a hybridization present. However,
in the cases of single hybridizations, its false positive rate remains low.
For all single hybridization events, Patterson’s D-Statistic and D3 perform
comparably, with the same pattern in false negative rates given hybridization
depth. However, in the case of n15h3, D3 outperforms Patterson’s D-Statistic
in terms of its higher precision and lower false positive rate. In this case,
Patterson’s D-Statistic has a very high false positive rate, of approximately
25%.

Figure 14 shows the false positive rate (orange), precision (black) and false
negative rate (red) for MSCquartets, HyDe, Patterson’s D-Statistic, and D3 on
the largest network under study (n25h5). For HyDe, an additional metric, wrong
hybrid rate (purple) describes the rate at which hybridization is detected but
is falsely attributed to the incorrect hybrid taxon. Again, all methods show a
high false negative rate and low false positive rate both of which could be
explain by a weakening of the hybridization signal when multiple
hybridizations are affecting the same taxa. All methods have a high precision
which means that when a hybrid is detected, it is very likely a true hybrid.
HyDe has slightly lower precision compared to other methods, but this is due
to the fact that HyDe (unlike other methods) test for a very specific parent-
hybrid relationship. When HyDe is tested in the setup of clear parent-hybrid
relationships (Figures 10, 11, and 17), HyDe indeed displays high precision.
It is notable that HyDe’s precision is better for n25h5 compared to n15h3 or
n10h2. This is due to the fact that the hybridizations in n10h2 and n15h3
involve ghost lineages which is not accounted for in HyDe. In this figure, all
test are Bonferroni-corrected at a level of significance
$\alpha=0.05/\textrm{number of tests}$, but we also show the uncorrected
version ($\alpha=0.05$) in Figure 22 in the Appendix.

Figure 14: False positive rate (orange), precision (black) and false negative
rate (red) for MSCquartets, HyDe, Patterson’s D-Statistic, and D3 on the
largest network tested with the highest number of hybridization events
(n25h5). All tests are Bonferroni-corrected at a level of significance
$\alpha=0.05/\textrm{number of tests}$ (see Figure 22 in the Appendix for
results on $\alpha=0.05$). For HyDe, an additional metric, Wrong Hybrid Rate
(purple) describes the rate at which hybridization is detected but is falsely
attributed to the incorrect hybrid taxon. High precision is noted for all
tests that evaluate for hybridization within a triple or quartet. All methods
also show a high false negative rate, and low false positive rate. Given that
HyDe tests specific parent-hybrid relationships, it has lower precision scores
than other methods.

###  3.2 Computational time

Each method is predicted to increase linearly with respect to the number of
gene trees or sequence length used as input. Due to the nature of summary
methods triple- or quartet-wise analysis, an increase in network size
corresponds to a cubic or quartic increase in time. These trends are displayed
in Figure 15.

Figure 15: CPU time of HyDe (left), MSCquartets (center), and TICR (right) as
network size and sequence length or number of gene trees change. Note the
different limits on the Y axis for HyDe, which is the fastest of all three
methods. For HyDe, there is a cubic increase in time is observed with respect
to network size. For MSCquartets, there is a quartic increase in time is
observed with respect to network size.

###  3.3 Nomiinae bee subfamily

Figure 16 displays the heatmap of the proportion of times that each taxon is
involved in a significant hybridization event as identified by MSCquartets,
HyDe, Patterson’s D-statistics and D3. This figure was created using ggtree
[36].

We selected the estimated gene trees from the IQtree2-GTRG model to use as
input in MSCquartets and display the proportion of times that each taxa is
involved in a significant hybridization (Figure 16 for the Bonferroni-
corrected significance level of $0.05/14950$, or $3.3\times 10^{6}$ and Figure
23 in the Appendix for significance level of $\alpha=0.05$) Stictonomia spp.
is implicated in all significant quartets, with Stictonomia schubotzi
appearing with the highest frequency. In addition, HyDe detects Stictonomia
spp. implicated 923 times over 2851 significant hybrid speciation events. We
also show the proportion that each taxon is identified as a parent (ancestral
lineage contributing genetic material to the hybrid) in the HyDe tests which
provides a broader picture of hybridization than any of the other methods. The
results of Patterson’s D-statistics or D3 align with those of HyDe and
MSCquartets in the identification of Stictonomia spp. (especially Stictonomia
schubotzi) as involved in hybridization events. The D-related tests, however,
cannot separate the hybrid taxon from the parents as HyDe. Though MSCquartets
primarily implicates Stictonomia spp. in hybridization, other methods find
widespread hybridization across the tree. These differences may be due to
methods’ sensitivity to the depth of hybridization events.

These results suggest that these closely related species may not be
reproductively isolated, which can lead to gene tree estimation error, and
difficulty in reconstructing the phylogenetic tree. In the original study [6],
gene tree estimation error was identified as a source of the discordance and
conflict. However, here we identify hybridization events as a plausible
explanation for the gene tree discordance.

Figure 16: Left: Phylogenetic tree of selected species from the bee subfamily
Nomiinae, as provided by [6]. Letters A-C identify nodes that were identified
as conflicting relationships by the analyses conducted in [6]. Right: Heatmap
of frequency of taxa identified as part of hybridization events (scaled
between 0 and 1 for comparison across methods). MSCquartets was conducted with
gene trees estimated using the IQTree2-GTRG model. Tests are Bonferroni-
corrected (see Figure 23 in the Appendix for significance level
$\alpha=0.05$). The clade with the highest proportion of proposed
hybridizations contains Stictonomia spp.

##  4 Discussion

Here, we present a deep investigation of the performance of genome-wide hybrid
detection methods. We found that all five methods compared (TICR [32],
MSCquartets [18, 19], HyDe [15], Patterson’s D-Statistic [24], Dp [10] and D3
[9]) have similar good performance (high precision and low false
positive/negative rates) on single shallow hybridizations involving few taxa
(n4h1 or n5h2). As more hybridizations are added involving similar groups of
taxa (n8h3, n10h2 and n15h3), all methods have a higher false negative rate
which suggests that combinations of gene flow events weaken the signal to
detect such hybridizations as opposed to creating discordant signal to
identify wrong hybridizations (which would have been evidenced by an increased
false positive rate). This is also confirmed by the results of TICR which is
unable to reject the major tree in most cases, except for those involving
single shallow hybridizations.

HyDe had a lower precision that other methods when ghost lineages were
involved (n10h2 and n15h3) which aligns with previous studies on the subject
[34, 23]. It also showed a higher rate of wrong hybrid identified within the
hybrid triple. However, HyDe is the only method able to detect the hybrid
taxon and the parent taxa involved in the hybridization as shown in the bee
dataset. The results of methods using site pattern frequencies or pairwise
differences is highly influenced by the topology of the underlying species
tree from which taxa arise. Longer coalescent times introduce noise to
sequence data such that comparison to a distant outgroup or comparison between
distant species is no longer advantageous, as the infinite-sites mutation
model on which Patterson’s D-Statistic is based expects a single mutation per
site [12]. With increased branch lengths, or increased distance between taxa,
convergent substitutions can cause ABBA and BABA (and other) site patterns can
accumulate [12].

Finally, we re-analyzed the dataset of the bee subfamily Nomiinae. While the
original study [6] concludes that gene estimation error could be the source of
discordance in the clade, here we show that hybridization is another plausible
explanation for the discordant patterns with all methods identifying
Stictonomia spp. (especially Stictonomia schubotzi) as involved in
hybridization events. As tools for discovering reticulations across the tree
of life are being developed and improved, so are tools to analyze and
manipulate complex networks. Method developers should seek to assess both
scalability and accuracy given growing, complex datasets.

##### Practical advice for evolutionary biologists.

From our investigation, we can conclude that MSCquartets [3] is an accurate
alternative to detect hybridization events under a variety of different
scenarios. HyDe is the only method that can distinguish which taxon is a
hybrid and not a parent of hybrid. However, HyDe cannot perform well when the
parents of hybridization are unsampled or extinct. All methods are unable to
detect hybridizations when multiple events are affecting the same set of taxa
or when hybridizations are deep. When a given parent-hybrid relationship wants
to be tested, HyDe outperforms D-Statistics-like tests by allowing the
identification of hybrid taxa vs parent taxa. Finally, TICR is a powerful
method to detect single shallow hybridization events, even if they involve
ghost lineages.

##### Differences with [14].

While working on this manuscript, another research group published a very
similar simulation study (see [14]). This work already addresses several of
the questions we try to address here. However, our work contains the following
notable extensions: 1) [14] did not investigate the performance of TICR [32]
or MSCquartets [3] which we do include here and we show that MSCquartets
outperforms most methods under varied hybridization scenarios; 2) the effect
of ghost lineages were not investigated in [14] which we study here and which
have a high impact in the performance of HyDe; 3) we study here more scenarios
involving shallow vs deep hybridizations as well as the effect of multiple
hybridizations, and 4) the networks under study in our work are larger and
have more hybridization events compared to those in [14].

##### Limitations and future work.

A limitation our comparison of the methods that take gene tree as input (TICR
and MSCquartets) to the methods that take site pattern frequencies as input
(HyDe, Patternson’s D, D3) is that we used true gene trees, which in practice
would be estimated from molecular sequences. Error in the estimated gene trees
could impact the estimation of hybridization events by MSCquartets or testing
for non-tree-like evolution by TICR, and thus these methods had an advantage
in our study. Along these lines, multiple sequence alignment errors occur in
phylogenomics data sets [39] could impact the performance of all methods
tested. We simulated data under a substitution-only model and thus all methods
were given true alignments as input. Lastly, the substitution model is well
behaved (homogeneous, stationary, and reversible), which is assumed by
HyDe;however, these assumptions can be violated in practice [21]. It is not
clear to what extent these practicalities will impact methods, and future work
should explore them in simulations and in real data sets.

##### Acknowledgements.

This work was supported by the Department of Energy [DE-SC0021016 to CSL] and
by the National Science Foundation [DEB-2144367 to CSL]. We thank Laura
Kubatko and Sungsik Kong for meaningful discussions about HyDe.

## References

  * [1] Elizabeth S. Allman, James H. Degnan, and John A. Rhodes.  Identifying the rooted species tree from the distribution of unrooted gene trees under the coalescent.  Journal of Mathematical Biology, 62(6):833–862, 2011.
  * [2] Elizabeth S. Allman, Jonathan D. Mitchell, and John A. Rhodes.  Gene tree discord, simplex plots, and statistical tests under the coalescent.  bioRxiv, 2020.
  * [3] Elizabeth S. Allman and John A. Rhodes.  Phlogenetic invariants.  In Reconstructing Evolution: New Mathematical and Computational Advances, pages 108–146. Oxford University Press, 2007.
  * [4] N. H. Barton and G. M. Hewitt.  Analysis of hybrid zones.  Annual Review of Ecology and Systematics, 16:113–148, 1985.
  * [5] Paul D Blischak, Julia Chifman, Andrea D Wolfe, and Laura S Kubatko.  HyDe: A Python Package for Genome-Scale Hybridization Detection.  Systematic Biology, 67(5):821–829, 03 2018.
  * [6] Silas Bossert, Elizabeth A Murray, Alain Pauly, Kyrylo Chernyshov, Seán G Brady, and Bryan N Danforth.  Gene Tree Estimation Error with Ultraconserved Elements: An Empirical Study on Pseudapis Bees.  Systematic Biology, 70(4):803–821, 12 2020.
  * [7] Remco Bouckaert, Timothy G. Vaughan, Joëlle Barido-Sottani, Sebastián Duchêne, Mathieu Fourment, Alexandra Gavryushkina, Joseph Heled, Graham Jones, Denise Kühnert, Nicola De Maio, Michael Matschiner, Fábio K. Mendes, Nicola F. Müller, Huw A. Ogilvie, Louis du Plessis, Alex Popinga, Andrew Rambaut, David Rasmussen, Igor Siveroni, Marc A. Suchard, Chieh-Hsi Wu, Dong Xie, Chi Zhang, Tanja Stadler, and Alexei J. Drummond.  Beast 2.5: An advanced software platform for bayesian evolutionary analysis.  PLOS Computational Biology, 15(4):1–28, 04 2019.
  * [8] Ruoyi Cai and Cécile Ané.  Assessing the fit of the multi-species network coalescent to multi-locus data.  Bioinformatics, 37(5):634–641, 12 2020.
  * [9] Matthew W Hahn and Mark S Hibbins.  A three-sample test for introgression.  Mol Biol Evol, 36(12):2878–2882, Dec 2019.
  * [10] Jennafer A P Hamlin, Mark S Hibbins, and Leonie C Moyle.  Assessing biological factors affecting postspeciation introgression.  Evol Lett, 4(2):137–154, Apr 2020.
  * [11] Mark S. Hibbins and Matthew W. Hahn.  The effects of introgression across thousands of quantitative traits revealed by gene expression in wild tomatoes.  PLOS Genetics, 17(11):1–20, 11 2021.
  * [12] Mark S Hibbins and Matthew W Hahn.  Phylogenomic approaches to detecting and characterizing introgression.  Genetics, 220(2), 11 2021.  iyab173.
  * [13] Richard R. Hudson.  Generating samples under a Wright–Fisher neutral model of genetic variation .  Bioinformatics, 18(2):337–338, 02 2002.
  * [14] Sungsik Kong and Laura S Kubatko.  Comparative Performance of Popular Methods for Hybrid Detection using Genomic Data.  Systematic Biology, 70(5):891–907, 01 2021.
  * [15] Laura S. Kubatko and Julia Chifman.  An invariants-based method for efficient identification of hybrid species from large-scale genomic data.  BMC Evolutionary Biology, 19(1):112, 2019.
  * [16] Wayne P. Maddison.  Gene Trees in Species Trees.  Systematic Biology, 46(3):523–536, 09 1997.
  * [17] Alexey Markin, Sanket Wagle, Tavis K Anderson, and Oliver Eulenstein.  RF-Net 2: fast inference of virus reassortment and hybridization networks.  Bioinformatics, 38(8):2144–2152, 02 2022.
  * [18] Jonathan D Mitchell, Elizabeth S Allman, and John A Rhodes.  Hypothesis testing near singularities and boundaries.  Electron J Stat, 13(1):2150–2193, 2019.
  * [19] Jonathan D. Mitchell, Elizabeth S. Allman, and John A. Rhodes.  Hypothesis testing near singularities and boundaries.  Electronic Journal of Statistics, 13(1):2150 – 2193, 2019.
  * [20] B.M.E. Moret, L. Nakhleh, T. Warnow, C.R. Linder, A. Tholse, A. Padolina, J. Sun, and R. Timme.  Phylogenetic networks: modeling, reconstructibility, and accuracy.  IEEE/ACM Transactions on Computational Biology and Bioinformatics, 1(1):13–23, 2004.
  * [21] Suha Naser-Khdour, Bui Quang Minh, Wenqi Zhang, Eric A Stone, and Robert Lanfear.  The Prevalence and Impact of Model Violations in Phylogenetic Analysis.  Genome Biology and Evolution, 11(12):3341–3352, 09 2019.
  * [22] Lam-Tung Nguyen, Heiko A Schmidt, Arndt von Haeseler, and Bui Quang Minh.  Iq-tree: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies.  Mol Biol Evol, 32(1):268–274, Jan 2015.
  * [23] Xiao-Xu Pang and Da-Yong Zhang.  Impact of Ghost Introgression on Coalescent-based Species Tree Inference and Estimation of Divergence Time.  Systematic Biology, 07 2022.  syac047.
  * [24] Nick Patterson, Priya Moorjani, Yontao Luo, Swapan Mallick, Nadin Rohland, Yiping Zhan, Teri Genschoreck, Teresa Webster, and David Reich.  Ancient admixture in human history.  Genetics, 192(3):1065–1093, 11 2012.
  * [25] Fernando Racimo, Sriram Sankararaman, Rasmus Nielsen, and Emilia Huerta-Sánchez.  Evidence for archaic adaptive introgression in humans.  Nat Rev Genet, 16(6):359–371, Jun 2015.
  * [26] Andrew Rambaut and Nicholas C. Grass.  Seq-Gen: an application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees.  Bioinformatics, 13(3):235–238, 06 1997.
  * [27] John A Rhodes, Hector Baños, Jonathan D Mitchell, and Elizabeth S Allman.  Mscquartets 1.0: quartet methods for species trees and networks under the multispecies coalescent model in r.  Bioinformatics, 37(12):1766–1768, Jul 2021.
  * [28] Sagi Snir and Satish Rao.  Quartet maxcut: A fast algorithm for amalgamating quartet trees.  Molecular Phylogenetics and Evolution, 62(1):1–8, 2012.
  * [29] Claudia Solís-Lemus and Cécile Ané.  Inferring phylogenetic networks with maximum pseudolikelihood under incomplete lineage sorting.  PLOS Genetics, 12(3):e1005896–, 03 2016.
  * [30] Claudia Solís-Lemus, Paul Bastide, and Cécile Ané.  PhyloNetworks: A Package for Phylogenetic Networks.  Molecular Biology and Evolution, 34(12):3292–3298, 09 2017.
  * [31] Alexandros Stamatakis.  Raxml version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies.  Bioinformatics, 30(9):1312–1313, May 2014.
  * [32] Noah W. M. Stenz, Bret Larget, David A. Baum, and Cécile Ané.  Exploring Tree-Like and Non-Tree-Like Patterns Using Genome Sequences: An Example Using the Inbreeding Plant Species Arabidopsis thaliana (L.) Heynh.  Systematic Biology, 64(5):809–823, 06 2015.
  * [33] Anton Suvorov, Bernard Y. Kim, Jeremy Wang, Ellie E. Armstrong, David Peede, Emmanuel R.R. D’Agostino, Donald K. Price, Peter J. Waddell, Michael Lang, Virginie Courtier-Orgogozo, Jean R. David, Dmitri Petrov, Daniel R. Matute, Daniel R. Schrider, and Aaron A. Comeault.  Widespread introgression across a phylogeny of 155 drosophila genomes.  Current Biology, 32(1):111–123.e5, 2022.
  * [34] Théo Tricou, Eric Tannier, and Damien M de Vienne.  Ghost Lineages Highly Influence the Interpretation of Introgression Tests.  Systematic Biology, 02 2022.  syac011.
  * [35] Dingqiao Wen, Yun Yu, Jiafan Zhu, and Luay Nakhleh.  Inferring phylogenetic networks using phylonet.  Systematic biology, 67(4):735–740, 2018.
  * [36] Guangchuang Yu, David K. Smith, Huachen Zhu, Yi Guan, and Tommy Tsan-Yuk Lam.  ggtree: an r package for visualization and annotation of phylogenetic trees with their covariates and other associated data.  Methods in Ecology and Evolution, 8(1):28–36, 2017.
  * [37] Yun Yu, James H. Degnan, and Luay Nakhleh.  The probability of a gene tree topology within a phylogenetic network with applications to hybridization detection.  PLOS Genetics, 8(4):1–10, 04 2012.
  * [38] Chao Zhang and Siavash Mirarab.  Weighting by gene tree uncertainty improves accuracy of quartet-based species trees.  bioRxiv, 2022.
  * [39] Chao Zhang, Yiming Zhao, Edward L. Braun, and Siavash Mirarab.  Taper: Pinpointing errors in multiple sequence alignments despite varying rates of evolution.  Methods in Ecology and Evolution, 12(11):2145–2158, 2021.
  * [40] Chi Zhang, Huw A Ogilvie, Alexei J Drummond, and Tanja Stadler.  Bayesian inference of species networks from multilocus sequence data.  Molecular biology and evolution, 35(2):504–517, 2018.

## Supplementary figures

### Reproducing the power calculations in [14]

Our original simulations involved only the networks n10h2, n15h3 and n25h5. We
noticed HyDe’s low precision and high false negative rates which contradicted
the positive results found in [14]. We then decided to reproduce the power
calculations in [14] as shown in Figure 17 and decided to include the networks
n4h1, n4h1${}_{\textrm{introgression}}$, n5h2 and n8h3 in the simulation
study. We discovered that HyDe has good performance when the parent-hybrid
relationship is clearly specified, and it is unable to detect hybrid
relationships involving ghost lineages.

Figure 17: Power of HyDe test on detecting a single hybridization event by
mixing parameter, using hybridlambda. This testing is performed on a gene
trees generated from a species topology
$(O:2.0,((P1:0.50,\\#H1:0::\gamma):0.50,(P2:0.5,(HYB:0.5)\\#H1:0::1-\gamma):0.5):1.0)$.
This topology has relatively shorter branch lengths than those tested in other
networks in this paper. It has four taxa and one hybridization event.

### Simulation figures

Figure 18: False positive rate (orange), precision (black) and false negative
rate (red) for MSCquartets, HyDe, Patterson’s D-Statistic, and D3 on the n4h1
network with mixing parameters $\gamma\in 0,0.1,0.2,0.3,0.4,0.5$, where 0
depicts no gene flow, and 0.5 depicts equal contribution from both parental
lineages. All tests are Bonferroni-corrected at a level of significance
$\alpha=0.05$. At a significance level of 0.05, non-zero false positive rates
can be seen in the case of zero gene flow for all but D3. As gene flow and
amount of genetic information increases, each method displays a decrease in
false negative rate. Figure 19: False positive rate (orange), precision
(black) and false negative rate (red) for MSCquartets, HyDe, Patterson’s
D-Statistic, and D3 on the networks: n4h1${}_{\textrm{introgression}}$
(network with single shallow hybridization event), n5h2 (network with two
overlapping hybridization events) and n8h3 (network with three overlapping
hybridization events). All tests are at a level of significance $\alpha=0.05$.
Each method performs well on the introgression event of n4h1, and n5h2 with
some false negatives seen with HyDe. Each method has a similar false negative
rate for h8h3, with precision near 100% for all but HyDe, which still had 75%
precision at minimum for all sequence lengths.  Figure 20: False positive rate
(orange), precision (black) and false negative rate (red) for MSCquartets,
HyDe, Patterson’s D-Statistic, and D3 on the networks: n10h2,
n10h1${}_{\textrm{shallow}}$ (n10h2 with only the shallow hybridization
event), n10h1${}_{\textrm{deep}}$ (n10h2 with only the deeper hybridization
event). All tests are at a level of significance $\alpha=0.05$. D3 performs
well on n10h2 and its singular hybridizations n10h1${}_{\textrm{shallow}}$ and
n10h1${}_{\textrm{deep}}$ relative to other methods given its high precision.
In contrast, Patterson’s D-Statistic has slightly higher false positives, at
the cost of precision. MSCquartets also has increased false positive and false
negative rates. HyDe has higher false negative rates for n10h2 and
n10h1${}_{\textrm{shallow}}$. There are also high false positives ( 20%) in
the case of n10h2, with a pattern for hybrid clade misidentification more
frequently than precise identification of hybrid clades. Figure 21: False
positive rate (orange), precision (black) and false negative rate (red) for
MSCquartets, HyDe, Patterson’s D-Statistic, and D3 on the networks: n15h3,
n15h1${}_{\textrm{deep}}$ (n15h3 with only the deepest hybridization event),
n15h1${}_{\textrm{shallow}}$ (n15h3 with only the shallowest hybridization
event), and n15h1${}_{\textrm{intermediate}}$ (n15h3 with a hybridization
event of intermediate depth). All tests are at a level of significance
$\alpha=0.05$. D3 has the highest and most stable precision across all n15h3
and single hybrid derivatives. False negative rates remain similar across all
tests, with the exception of Patterson’s D-Statistic on n15h3, which is
slightly lower, paired with a higher false positive rate. Patterson’s
D-Statistic also shows decreased precision scores. MSCquartets’ precision
remains high ($\>80\%$) in all cases except n15h1${}_{\textrm{intermediate}}$.
HyDe also has a lower precision, around $10\%-25\%$ across all network
formats. Figure 22: False positive rate (orange), precision (black) and false
negative rate (red) for MSCquartets, HyDe, Patterson’s D-Statistic, and D3 on
the largest network tested with the highest number of hybridization events
(n25h5). All tests are at a level of significance $\alpha=0.05$. MSCquartets,
Patterson’s D-Statistic, and D3 all behave comparably on n25h5, the largest
tested network with three shallow hybridization events. HyDe also shows higher
relative precision and lower rates of hybrid clade misattribution. Figure 23:
Left: Phylogenetic tree of selected species from the bee subfamily Nomiinae,
as provided by [6]. Letters A-C identify nodes that were identified as
conflicting relationships by the analyses conducted in [6]. Right: Heatmap of
frequency of taxa identified as part of hybridization events (scaled between 0
and 1 for comparison across methods). MSCquartets was conducted with gene
trees estimated using the IQTree2-GTRG model. Tests are performed at
significance level $\alpha=0.05$.