A global test of hybrid ancestry from genome-scale data

1 Introduction

Advances in modern sequencing technologies have led to the availability of a substantial amount of DNA sequence data from hundreds, if not thousands, of species. These data are often used to understand the evolutionary relationships between species as well as the process of evolution more generally. One way of representing the evolutionary relationships between a collection of species is with a species-level phylogeny, also known as a phylogenetic tree. A species-level phylogeny is an acyclic graph in which all edges are considered to be directed from the root node to the leaves. An internal node at which a lineage splits into two lineages that each subsequently evolve independently represents a speciation event. In a phylogeny, all internal nodes except the root node are connected to three edges, making the degree of all internal nodes three. The root node, which represents the common ancestor of all the species in the data, is connected with two edges leading the root node to have degree two, while the leaves (also known as taxa) have degree one and represent the present-day populations from which the DNA sequence data were collected (see right panel of Figure 1).

Figure 1:

Model for the species-level relationships among four taxa under the coalescent model with hybridization. Here taxon H is a hybrid of taxa P ₁ and P ₂, τ ₁, τ ₂, and τ ₃ are the speciation times for the internal nodes from the present. The circles on the leftmost figure represent internal nodes. The red circle (n ₂) represents the hybrid node as there are three descendent from that node. The network on the left can be decomposed into two species-level phylogenies; one with H and P ₂ as sister taxa with probability γ, and one with H and P ₁ as sister taxa with probability 1 − γ.

One major drawback of using a phylogenetic tree to represent evolutionary history is that it assumes speciation occurs through bifurcation, meaning that an ancestor can have precisely two descendants. However, this might not always be the case since there are other evolutionary mechanisms that generate novel species. One such mechanism is hybridization, which occurs when two distinct populations interbreed and produce a new species that shares genetic material with both parents. However, hybridization does not always result in creation of a new species, in which case the mechanism is known as introgression (Dowling and DeMarais 1993; Good et al. 2003; Grant et al. 2004; Mallet 2005, 2007; Roques et al. 2001; Salzburger et al. 2002; Thórsson et al. 2001; Weigel et al. 2002). Recent studies have shown numerous instances of hybridization in plants and animals, despite the belief that hybridization is rare (Mallet 2007; Mavárez et al. 2006; Rieseberg 1997). If a collection of species is known to include species that have arisen via hybridization, one should infer a species network rather than a phylogenetic tree. Since the internal nodes in a network are allowed to have degree more than three, networks can capture evolutionary processes such as hybridization (see left panel of Figure 1). However, inferring a network from data can be computationally intensive, especially when the number of taxa and/or the length of the DNA sequence alignment are large (Rabier et al. 2021). In addition, network analysis in the absence of hybridization can lead to incorrect conclusions concerning the ancestry of the taxa identified to be of hybrid origin. Thus, an important first step before deciding whether to perform a network or a tree analysis is to gather concrete evidence about the possibility of hybridization among the species on which inference is to be made.

In the literature, methods for detecting hybridization are ubiquitous. However, several of these methods ignore a crucial evolutionary process called incomplete lineage sorting (ILS). In practice, the DNA sequence data come from genes, and each of these genes has its own evolutionary history, which may or may not match the evolutionary history of the species. When two gene copies from the current generation fail to coalesce in the immediately previous generation, incomplete lineage sorting (ILS) results, leading to gene-tree species-tree discordance (Maddison 1997). Coalescent theory links the gene trees to a species tree by providing the probability distribution of all possible rooted gene trees arising within a rooted species tree and thus can be used to model ILS (Pamilo and Nei 1988; Tajima 1983; Takahata and Nei 1985). To distinguish between hybridization and incomplete lineage sorting, Joly et al. (2009) proposed a statistic based on the observation that the genetic distance between two DNA sequences should be smaller in the presence of hybridization than the distance that would be observed if only ILS is acting on the genes. One can simulate this distance under the null hypothesis that ILS alone explains the discordance using coalescent theory. Once the distribution of the distance under the null hypothesis is estimated, one can compare the observed distance to the null distribution to assess the evidence for hybridization (Joly et al. 2009).

Tests based on the specific site pattern frequencies (a particular position at an alignment is called a site, and the combination of possible nucleotides at a particular site for the sampled individuals is called a site pattern) have also been developed to test ancient admixture utilizing Patterson’s D-Statistic (Durand et al. 2011; Green et al. 2010; Patterson et al. 2012). In addition, several likelihood-based methods have been proposed, including Meng and Kubatko (2009), Kubatko (2009), and Yu et al. (2014). In 2019, Kubatko and Chifman proposed a method to detect hybridization based on “phylogenetic invariants” (polynomials in the site pattern probabilities that evaluate to zero on any probability distribution that is consistent with the tree topology and associated model). Their method considers models that can generate data that allows the possibility of hybridization as well as incomplete lineage sorting. Based on a comparison of popular hybridization detection methods, Kong and Kubatko (2021) concluded that Kubatko and Chifman’s (2019) invariant-based method detects hybrid species quickly and accurately, generally outperforming the other methods tested in terms of power.

One common characteristic of several of these available methods is that they test whether a suspected species is a hybrid given knowledge of its two parents. However, to decide whether or not a network analysis is necessary, we need to perform a global test to assess evidence for the presence of one or more hybrid taxa in the entire data set. Kubatko and Chifman (2019) extended their method to deal with phylogenies containing more than four taxa by testing every species as a hybrid and applying a Bonferroni correction to adjust for multiple testing. However, the number of tests one needs to consider can be massive, especially when the number of taxa is large, making the Bonferroni correction less powerful in detecting hybridization events. For instance, from a data set containing m species, one can pick three species and test three hypotheses considering one of the three species as a hybrid and the other two as parents. Hence, one can test a total of 3 × m 3 hypothesis using a data set containing m species. If, for instance, a data set contains eight species, we need to perform 168 tests, and the number of required tests grows significantly as the number of taxa increases. In addition, the individual tests are correlated because the test statistic for each individual test is calculated using sequences from closely related species, making the Bonferroni method conservative. Thus, to perform a global test for detecting hybridization, we need a method that can handle arbitrary correlation structures. One possibility is to use a combination test.

A combination test combines the individual p-values to aggregate multiple small effects rather than adjusting individual p-values to maintain the false rejection rate. Fisher’s combination test is the oldest method of combining individual p-values to test an overall null hypothesis (Fisher 1932). Consider a scenario where we carry out l individual tests. Fisher’s combination test statistic to perform this test is as follows:

where p _i , i = 1, 2, …, l, is the p-value from the ith individual test. The test statistic above converges to a Chi-squared distribution with 2l degrees of freedom under the assumption that all individual nulls are true and the p _i are independent. When the individual p-values are small, the value of the test statistic in Equation (1) tends to be larger, leading to evidence suggesting that not all of the null hypotheses are true. The test has substantial power to reject the global null hypothesis when a large fraction of the individual p-values are small. This will not be the case, however, when only a small fraction of the individual tests are expected to be significant (Arias-Castro et al. 2011; Koziol and Perlman 1978), a situation that we might often encounter in many applications of bioinformatics. For instance, while assessing genome-scale data to detect hybridization events for 20 taxa, one needs to perform a total of 2280 tests, only a few of which are expected to be significant. The application of Fisher’s combination test, in this case, may lead to failure to detect any hybrid species.

Multiple attempts have been made to improve the power of the test in situations like this. The most popular combination tests that attempt to improve power are Tippett’s minimum p-value (Tippett et al. 1931) and higher criticism tests (Donoho and Jin 2004). In Tippett’s minimum p-value approach, the test statistic is determined by the minimum of the individual p-values. That is, for the l individual p-values (p ₁, p ₂, …, p _l ), the test statistic is defined as: S _T = min(p ₁, p ₂, …, p _l ). Assuming all the individual p-values are independent and all the individual nulls are true, Tippett’s minimum p-value statistic (S _T ) converges to a Beta(1, l) distribution. For l individual tests with α = 0.05, the higher criticism test statistic proposed by Tukey is

In Tukey’s interpretation, the null hypothesis should be rejected when HC is greater than 2, assuming that the individual p-values are independent. When the assumption of independence of individual p-values holds, Tippett’s minimum p-value and the higher criticism tests both have high power. However, both these tests, as well as Fisher’s combination test, fail to perform well when there is dependence among the individual tests.

In reality, it is often the case that the individual tests have a dependence structure and the p-values are correlated. For example, in a phylogenetic tree, species are related by their common ancestor. There are methods available based on permutation or numerical simulation to incorporate the correlation structure (Liu and Xie 2019). However, these approaches are computationally involved and time-consuming when performing a large number of individual tests involving massive data sets. Recently, an analytical approximation method has been proposed using the higher criticism test to incorporate the dependence structure (Barnett et al. 2017). However, this approximation is computationally involved and does not perform well when the p-values are extremely small.

To incorporate the arbitrary dependence structure of individual tests and efficiently handle a large number of tests, Liu and Xie (2020) proposed a combination test called the Cauchy combination test (CCT). Following Fisher’s combination test, the test statistic of the Cauchy combination test is calculated using the weighted sum of transformed p-values. Let p _i be the individual p-value, for i = 1, 2, …, l. Then, the Cauchy combination test statistic is defined as

where the weights ω _i ’s are nonnegative and ∑ i = 1 l ω i = 1 . The authors showed that the statistic in Equation (2) converges to a standard Cauchy distribution under the global null hypothesis even when there is an arbitrary correlation structure between the individual p-values. The only assumption made by the authors is that any two of the individual test statistics should have a bivariate normal distribution. That is, considering X _m to be the mth test statistic, m = 1, 2, …, l, for any 1 ≤ i < j ≤ l , ( X i , X j ) T follows a bivariate normal distribution. However, Liu and Xie (2020) argued that the bivariate normality assumption is mild, and the test is robust even when the assumption is not satisfied.

In 2022, Zhongxue Chen showed situations where the Cauchy combination test failed to reject the global null hypothesis even if some of the individual tests were rejected correctly. The author proposed a robust test combining the MinP (Tippett et al. 1931) and the Cauchy combination (Liu and Xie 2020) test that performs well under a wide range of situations. Consider combining l individual tests with p-values p ₁, p ₂, …, p _l . Let the p-value for the global null hypothesis using the Cauchy combination test be p _c , and the p-value using the MinP test be p _m with p _m = l × min(p ₁, p ₂, …, p _l ). Using the resulting p-values from the Cauchy combination and MinP test, Chen (2022) proposed a test statistic called MinP-CCT-MinP (MCM) with which the p-value for the global test is calculated as follows:

Chen (2022) showed that the proposed test is robust and performs well in situations where the Cauchy combination test might fail while maintaining the correct type I error rate. In addition, the author proved that, unlike the Cauchy combination test, the MCM test need not assume that the joint distribution of pairs of individual tests is bivariate normal.

Here, we consider the test statistic proposed by Kubatko and Chifman (2019) to test whether a particular species is a hybrid of two specified parents using both the Cauchy combination and MCM test frameworks proposed by Liu and Xie (2020) and Chen (2022), respectively, to suggest a global test statistic for detecting the presence of hybridization within a set of species. We perform an extensive simulation study to investigate the performance of the proposed test. We also apply our method to empirical data sets and compare results with existing studies of those data sets. We focus our work on the test of Kubatko and Chifman (2019), rather than on other available methods e.g., Patterson’s D-Statistics (Durand et al. 2011; Green et al. 2010; Patterson et al. 2012), because previous work (Kong and Kubatko 2021) has shown that the method of Kubatko and Chifman (2019) is more powerful than methods based on D-Statistics. Recalling Equations (2) and (3), it is clear that if we compare the power of a global test using CCT or MCM for two different individual tests, the individual test that is more powerful will always result in a more powerful global test, and it is thus unnecessary to further consider tests based on D-statistics here. Nonetheless, we note that the method could equally well be applied to the case where D-Statistics are used for the individual tests. Finally, we note that a key feature of our approach is that dependence in the individual tests is allowed; i.e., there is no assumption that the p-values from the individual tests are uniformly distributed when the null hypotheses are true. Indeed, an empirical data set for which no hybrids are believed to exist displays a non-uniform distribution of p-values.

2 Methods

We begin this section by providing the mathematical details of both the individual test of hybridization on quartets of species and the global test of hybridization for the entire set of species. Readers who are most interested in the empirical aspects of the proposed method of hybrid detection can skip the Methods and move straight to the Results.

As discussed in the previous section, ILS is a crucial evolutionary mechanism that can lead to gene-tree species-tree discordance and thus should be incorporated into phylogenetic applications. Hence, we will consider a data generation mechanism that allows both hybridization and ILS, following the model proposed by Meng and Kubatko (2009). As already discussed, a bifurcating species-level phylogeny is inadequate to represent a hybridization event. Instead, a species network that connects two edges of a phylogenetic tree using a horizontal line is commonly used to depict a hybrid node. Consider the phylogenetic network in the leftmost panel of Figure 1. Here H is a hybrid species that inherits a portion (γ) of its genetic material from P ₂ and the rest (1 − γ) from parent P ₁; τ _i is the speciation time for node i measured from the present.

We assume that the aligned DNA sequences arise along the network. Let Y _L be the observed nucleotide for a specific site in the DNA sequence for taxon L. Then a site pattern arising from the species network in Figure 1 can be expressed as Y O Y P 1 Y H Y P 2 . For a four-taxon network, the number of possible site patterns is 4⁴ = 256 since Y _L ∈ {A, G, C, T}. Let p ijkl = Pr ( Y O = i , Y P 1 = j , Y H = k , Y P 2 = l | ( S γ , τ ) ) be the probability of site pattern ijkl with i, j, k, l ∈ {A, C, G, T}, and let P = (P _AAAA |(S _γ , τ), …, P _TTTT |(S _γ , τ)) be the set of all possible site pattern probabilities on network (S _γ , τ). If we let N _ijkl be the observed number of sites with site pattern ijkl and N be the total number of sites in the aligned DNA sequence data, one can estimate p _ijkl as p ̂ ijkl = N ijkl N , and the set of observed site pattern probabilities can be written as P ̂ = ( P ̂ AAAA | ( S γ , τ ) , … , P ̂ TTTT | ( S γ , τ ) ) .

As the vector N P ̂ represents the observed count of each of the 256 possible site patterns, it is easy to see that

With i ≠ j ∈ {A, G, C, T} consider the following linear relationships:

Kubatko and Chifman (2019) showed that f ₂ and f ₄ evaluate to zero assuming the site patterns arise from tree S ₁ (middle panel of Figure 1), while f ₁ and f ₃ are non-zero, and the opposite is true considering f ₁ and f ₃ with species tree S ₂ (rightmost panel of Figure 1). The authors also showed that f 1 f 2 = γ 1 − γ , and using this relationship, they defined a test statistic called the Hils statistic to test H ₀: γ = 0:

Here, μ f 1 and μ f 2 are the mean of f ₁ and f ₂, respectively, σ f 1 and σ f 2 are the standard deviation of f ₁ and f ₂, respectively, and σ f 1 f 2 is the covariance between f ₁ and f ₂. It has been shown that under H ₀, the statistic H follows a standard normal distribution and the null is rejected for large values of the observed statistic (Kubatko and Chifman 2019).

Now, consider testing the following hypothesis given aligned DNA sequence data for m + 1 taxa (m ingroup species and one outgroup):

• H ₀: None of the species are hybrid versus H _A : At least one of the species is hybrid

Consider choosing three ingroup taxa and the outgroup making a four-taxon subset from the set of taxa. Using the selected subset, we can consider a quartet network assuming that one of three ingroup taxa is a hybrid and the other two are the parental species. Note that we can build three such networks from one four-taxon subset by designating each ingroup taxon as the hypothesized hybrid species. For instance, if the selected four-taxon subset contains outgroup O and ingroup taxa P ₁, P ₂, and H, we can build three networks assuming H as a hybrid of parents P ₁ and P ₂, P ₁ as a hybrid of parents H and P ₂, and P ₂ as a hybrid of parents H and P ₁. Using the whole collection of species, one can build 3 × m 3 quartet networks with the outgroup and three taxa (triplet) selected from the m ingroup taxa. For instance, one can build a total of 3 × 5 3 = 30 quartet networks from a collection of species containing one outgroup and five taxa.

Note that for a network with m taxa, 4 ^m site patterns are possible. To get the site patterns for a quartet network built by selecting a four-taxon subset from the collection of taxa, one can marginalize over the missing species. For instance, consider creating quartet networks from the species network (S _γ , τ) in Figure 2a containing six species. Let the kth quartet network ( S γ k ) contains species O as outgroup, h as hybrid, and p ₁, p ₂ as the parental species. Then the site pattern probability P AAGT S γ k | ( S γ , τ ) can be found by:

here, i, j ∈ {A, G, C, T} represents the possible site patterns of species a and b from network (S _γ , τ) in Figure 2a that were not selected in the kth quartet network ( S γ k ) , respectively. It is clear that after marginalization, the possible number of site patterns for a quartet network is 4⁴ = 256.

Figure 2:

Five hybridization scenarios were considered in this study: (a) one recent hybridization event within six species, (b) one ancient hybridization event within six species, (c) one recent and one ancient hybridization event within seven species, (d) one recent hybridization event within twenty species, and (e) one ancient hybridization event within twenty species. The outgroup species is denoted by o. For panels a, b, and d, the hybrid species is denoted by h, and parents of the hybrid species are denoted by p ₁ and p ₂. For panel c, the two hybrid species are denoted by h ₁ and h ₂; the parents of h ₂ are p ₁ and p ₂, and the parents of h ₁ are the ancestor of a and b, and the ancestor of p ₁, h ₂, and p ₂. For panel e, the hybrid species is denoted by h, and its parents are the ancestor of a, b, …, l, and the ancestor of m, n, …, u. γ is the proportion of genetic material shared by the leftmost parent to the hybrid. The time between speciation event i and i + 1 is denoted by τ _i . Note that only the first sub-figure (panel a) is fully labeled.

Once the set of four taxa for the kth (k = 1, 2, …, l = 3 × m 3 ) quartet network is selected, one can decompose the network ( S γ k ) into two species trees (S _k1 and S _k2) following Figure 1. Also, following the idea of Kubatko and Chifman (2019), one can define f 1 S γ k , f 2 S γ k , f 3 S γ k , and f 4 S γ k using the quartet network S γ k and show that f 1 S γ k and f 3 S γ k are zero when evaluated on the site pattern probabilities arising from one of the species trees (e.g., S _k1) and not zero when evaluated on the site pattern probabilities arising from the other species tree (e.g., S _k2). On the other hand, for f 2 S γ k and f 4 S γ k , the opposite is true. However, none of the linear relationships are zero when evaluated on the site pattern probabilities from the hybrid species network S γ k with γ _k ∈ (0, 1). One can obtain a function of γ _k using these linear relationships as follows:

Now, one can easily estimate these functions by substituting the observed site pattern probabilities and then use the properties of the multinomial distribution and easily estimate the means μ f 1 S k γ , μ f 2 S k γ and variances σ f 1 S k γ 2 , σ f 2 S k γ 2 of these estimated functions. It is also possible to estimate the covariance σ f 1 S k γ , f 2 S k γ between these functions. Once the means, variances, and covariances are computed, one can use the H statistic in Equation (4) to test H ₀:γ _k = 0 for the kth quartet network ( S γ k ) :

Under the null hypothesis that γ _k = 0, the term μ f 1 S γ k μ f 2 S γ k is 0. Also, note that under the null hypothesis, H ^k follows the standard normal distribution. Thus, the p-value of this test can be computed easily as follows:

where h ^k is the observed value of the test statistic and Φ(h ^k ) = Pr(Z < h ^k ) with Z ∼ normal(0, 1). Note that one can repeat this process and compute the p-values for all possible tests. For instance, for a species network with five ingroup and one outgroup species, one can compute a total of l = 3 × 5 3 = 30 p-values. These tests will not be independent as species in a tree are related through their common ancestors. Let p ¹, p ², …, p ^l be the corresponding p-values of the l tests. To test the global null hypothesis: H ₀: None of the species are hybrid, we propose the following statistic:

here, the ω _i ’s are non-negative weights such that ∑ k = 1 n ω i = 1 . Under the global null hypothesis H ₀, the test statistic G follows a standard Cauchy distribution (Liu and Xie 2020). Note that a large value of the test statistic G provides evidence against the global null hypothesis. Thus the c.d.f. of the standard Cauchy distribution can be used to compute the p-value for the global null as follows:

Here, g is the realized value of G. The global null hypothesis is rejected if the p-value is smaller than α given that the chosen significance level is α. Note that in applying the results of the Cauchy combination test (Liu and Xie 2020), we have assumed that for any 1 ≤ i < j ≤ l, ( H i , H j ) T follows a bivariate normal distribution. The bivariate normality assumption is reasonable as the site pattern probabilities are well approximated by a normal distribution when the number of sites is large. Also, the summation of any two test statistics (e.g., H ¹ + H ³) is just the sum of two approximate normal random variables, and thus the sum is also normal. Hence, it is safe to assume bivariate normality here.

Finally, to apply the MCM proposed by Chen (2022), we do the following:

1. Step 1: Using the individual p-values (p ₁, p ₂, …, p _l ), calculate the p-value (p _m ) for the global test using the MinP test with p _m = l × min{p ₁, p ₂, …, p _l }.

2. Step 2: Calculate the p-value (p _mcm ) for the global test of no hybridization event using the MCM test with p _mcm = 2 × min{p _c , p _m , 0.5}, where p _c is the p-value for the global test from the Cauchy combination test.

3 Results

3.1 Simulation study

To evaluate the performance of the proposed statistic using the Cauchy Combination test, we plot the power curves for each of the scenarios (see Figures 5–8 in the Appendix for results under all of the simulation settings). As a representative example, panel (a) of Figure 3 shows power plots of the proposed test for the twenty-taxon tree with one ancient hybridization event (Figure 2e) with an α = 0.05. When γ = 0, these plots show that the power of the test is less than α = 0.05 across all the scenarios; i.e., under the null hypothesis of no hybridization, the power of the test is less than the selected level of significance. This result shows that the proposed test has the correct type-I error rate. Also, it is noticeable from the results that regardless of other parameters, as the length of the sequence increased, the power of the test also increased (setups with the length of sequence equal to 200k are the most powerful, and the power gradually decreases as the length decreases to 50k).

Figure 3:

Results of the power simulations for the twenty-taxon hybrid species network with one ancient hybridization event in (Figure 2e) using the Cauchy combination test (panel a) and the MCM test (panel b). Four different DNA sequence lengths (50k, 100k, 150k, and, 200k), four different effective population sizes (θ = 0.001, 0.005, 0.01, and, 0.03), and two different branch lengths (τ = 0.50 and 1.0) have been considered. Here, γ represents the proportion of genetic material contributed by each parent to the hybrid taxon.

The effect of branch length is also as expected. Holding everything else constant, tests for data with branch lengths equal to 1.0 are more powerful than those with 0.50. This result is expected because the prevalence of ILS increases in gene trees as branch length decreases, making hybridization harder to detect. We see a similar effect of the parameter θ, which controls the mutation rate (higher values correspond to more mutations). As expected, decreased value of θ results in decreased power to detect hybridization when other parameters are held constant. In addition, the power of the test is also increasing as the value of γ increases, regardless of other parameters. This result is also expected since a bigger value of γ is associated with a greater contribution from each parent, making it easier to detect hybridization. Overall, the test performs well. We can see that in almost all of the cases, the test reaches a high power (0.9–1.0) given a reasonably large sequence length (100 k or more) when the value of γ is 0.25 or greater.

It is also worth noting that the test adheres to all theoretical expectations. The only case when the proposed test does not perform well is when both the branch length and parameter θ are extremely small. As can be seen in the first panel of Figure 3, the test has extremely low power to detect hybridization when the branch lengths are 0.50 and θ = 0.001 even when the sequence length is high and the value of γ is high. However, this is expected because, with branch lengths of 0.50, one should expect a large amount of ILS and an extremely low mutation rate as θ is extremely small. The combination of these two extreme situations makes it extremely hard to detect hybridization.

The performance of the global test using the MCM statistic proposed by Chen (2022) is similar to that of the test using the Cauchy combination test statistic proposed by Liu and Xie (2020). As a representative example, panel (b) of Figure 3 shows power plots of the proposed test for the twenty-taxon tree with one ancient hybridization event (Figure 2e) with an α = 0.05. Figure 3 shows that the global test satisfies the theoretical expectations. We observed that the power of the test increases as the branch length (τ), length of DNA sequences (N), hybridization parameter (γ), and effective population size (θ) increase when the effect of each of these parameters are observed fixing all the other parameters.

We noticed that the power of the global test using the MCM test is slightly higher than the Cauchy combination test. For instance, the power of the test reaches 1.0 for a slightly lower value of γ when compared with the Cauchy combination test. Also, the power of the test using the MCM statistic is slightly higher in some cases when the length of the DNA sequence is small (e.g., N = 50k). However, the type-I error rate under the MCM test framework was greater than the chosen level of significance (α) in a few cases (7 out of 128 cases) (Figure 4), while the global test using the Cauchy combination framework maintained the correct type-I error rate in all the scenarios.

Figure 4:

Type-I error rate under the Cauchy combination test (panel a) and the MCM test (panel b) for all the scenarios (T20_deep: twenty species with one ancient hybridization, T20_recent: twenty species with one recent hybridization, T6_deep: six species with one ancient hybridization, and T6_recent: six species with one recent hybridization). The red dots are situations where the type-I error rates were greater than the chosen level of significance (α = 0.05). For each setup, there are 32 simulated type-I error rates considering lengths of the DNA sequence alignment (50k, 100k, 150k, and, 200k), branch lengths (τ = 0.50 and τ = 1.0), and effective population sizes (θ = 0.001, 0.005, 0.01, and, 0.03). Each of these box plots consists of these 32 type-I error rates.

4 Discussion

We have proposed two methods for assessing whether a data set contains one or more hybrid species that take incomplete lineage sorting into account. Based on DNA sequence data on a set of taxa, we propose our test as the first step before inferring a phylogenetic tree or network. One of the key benefits of the proposed method is its speed. The computation time of the proposed test is almost negligible since we used observed site pattern frequencies to calculate the test statistic. The speed makes the test feasible for genome-scale data sets with a large number of taxa. On one hand, individual tests are based on observed site patterns derived from a multinomial distribution, which allows a normal distribution to be derived from the asymptotic distribution of the site patterns for the test statistic. On the other hand, the null distribution of the global test was standard Cauchy, leading to a straightforward method for calculating the p-value.

Our simulation studies show that the methods are powerful to detect any presence of hybridization events even when the hybridization event is close to the root. Results show that deep or ancient hybridization events are harder to detect compared to recent hybridization events. However, if the branch length and length of the DNA sequence are large enough, the proposed method can still detect the presence of hybridization with high power. The proposed methods are extremely efficient: for six taxon networks with 250,000 sites, the whole process can be completed within 25 s, while for a data set with 20 taxa and one million sites, the whole process can be completed within 10 min. This shows that the proposed methods are extremely scalable for large data sets with a large number of sites. It is also possible to use parallel computing for the individual tests and then combine the results at the end if more time efficiency is needed.

The methods are based on the Hils statistic proposed by Kubatko and Chifman (2019) in combination with two combination tests (the Cauchy combination test and the MCM test). We notice that the type I error rate under the MCM test is sometimes slightly larger than the chosen level of significance, and the Cauchy combination test is more conservative. Thus, there is a possibility of improvement by making the type I error rate smaller for the MCM test while maintaining the same power. There are also opportunities to try to reduce the number of individual tests needed to perform the global test. As we are only interested to know whether one or more hybrid taxa are present in a data set, it is reasonable to perform one test per quartet tree to determine if there are any hybrid species in that set and finally combine the results.

We note that the proposed framework can be considered as an example of solving problems where it is difficult to perform a global test but easy to perform individual tests that have an unknown dependence structure. The framework can be followed in any field of application where a global test is desired, but only p-values from individual tests are available. Example applications can include meta-analysis, where p-values from different studies are combined to test a global hypothesis, finding rare genes, and testing more general phylogenetic hypotheses. Also, as we have mentioned in the Introduction, any individual test (e.g., Patterson’s D-Statistics) can be used instead of the test proposed by Kubatko and Chifman (2019) to test the global null of the absence of hybrid species in a dataset. However, because of the dependent nature of a phylogenetic tree, these individual tests will always have a non-negligible dependence structure. In this paper, we solve the problem of using a single hypothesis test for the entire collection of species, removing the need to test subsets separately and aggregate the results in an ad-hoc manner.

5 Conclusions

In many areas of biological research, it is central to understand the relationships between species. There are two primary ways to represent such relationships: a species-level phylogeny and a species network. Therefore, depending on the data, it is necessary to decide which analysis method to use. In doing so, we face several challenges, such as tests hypothesis in large genomic data sets in a computational efficient manner, modeling the data at the gene level as well as species level, and developing methods to account for potential evolutionary mechanisms (e.g., incomplete lineage sorting, gene duplication, hybridization) that generate incongruence between individual gene histories. The multispecies coalescent model is commonly used to model incomplete lineage sorting and has also been used to develop a framework for testing for hybridization. We use existing individual tests for hybridization together with combination tests to develop a novel approach for testing for the presence of hybridization within a given collection of species. We demonstrate that the proposed method is powerful in detecting hybridization events in several possible situations via simulation. We also demonstrate the method is capable of detecting hybridization events in empirical data. It is also notable that our method is extremely efficient and thus can handle genome-scale data with a large number of taxa. Thus, the proposed method provides researchers with an efficient way of detecting hybridization events in a data set and formally deciding whether a tree or a network analysis is more appropriate for a given data set.