Population size differences can lead to biases in phylogenetic inference and introgression detection in the presence of purifying selection

Supplementary Materials Population Size Differences Can Lead to Biases in Phylogenetic Inference and Introgression Detection in the Presence of Purifying Selection Chong He1*, Meng-Yun Chen2, Hao Zhu134* 1 Bioinformatics Section, School of Basic Medical Sciences, Southern Medical University, Guangzhou, 510515, China 2 Guangzhou Key Laboratory of Subtropical Biodiversity and Biomonitoring, Guangdong Provincial Key Laboratory of Biotechnology for Plant Development, School of Life Sciences, South China Normal University, Guangzhou 510631, China 3 Guangdong-Hong Kong-Macao Greater Bay Area Center for Brain Science and Brain-Inspired Intelligence, Southern Medical University, Guangzhou, 510515, China 4 Guangdong Provincial Key Lab of Single Cell Technology and Application, Southern Medical University, Guangzhou, 510515, China *Corresponding author: [email protected]; [email protected] Supplementary Methods Examining the Values of P(V1 | V) under Neutral Evolution and under Purifying Selection P(V1 | V) is the probability that a genealogical tree is a type V1 genealogical tree given that this genealogical tree is a type V genealogical tree (Fig. 7). Under neutral evolution, it is theoretically expected that P(V1 | V) = 2/3 (Fig. 9B, Supplementary Appendix 2). The SLiM simulation described here aimed to empirically examine whether P(V1 | V) is indeed 2/3 under evolution. We also examined whether the value of P(V1 | V) under purifying selection is greater or smaller than that under neutral evolution. We set the simplest biallelic context in our simulations because our aim was to determine the direction (greater or smaller) of the difference between the values of P(V1 | V) under purifying selection and under neutral evolution. If more than two alleles are involved, the magnitude of the difference can change, nevertheless, but the direction of the difference will not change. The “Mutation.time” and “Mutation.node” attributes in the Python package tskit were used to extract information about the time each mutation occurs and the branch each mutation falls on. We analyzed gene trees that experienced deep coalescence and at the same time contained one and only one mutation occurring at a time earlier than time t123 (this mutation corresponds to the mutation M in our theoretical analysis). Among these gene trees, we calculated the ratio between gene trees where the mutation occurring in S123 falls on the first diverging branch (type V1) and gene trees where the mutation occurring in S123 falls on any of three terminal branches (type V). Supplementary Appendix S1 Here, we provide a more rigorous explanation of why under patterns A3, A4, and A5, P(V1 | Ax) = P(V1 | V) (x = 3, 4, 5). As patterns A3, A4, and A5 are similar, here, we take pattern A3 as an example (g1 and g2 are connected to “noncarriers” and g3 is connected to a “carrier” at t123, Fig. 7). If pattern A3 occurs, G123 is certainly a type V genealogical tree. Hence, P(V1 | A3) may be written as P(V1 | A3V). Then, according to Bayes’ theorem, P(V1 | A3V) may be written as follows: 𝑃(𝑉1|𝐴3)=𝑃(𝑉1|𝐴3𝑉)=𝑃(𝑉1𝐴3|𝑉) 𝑃(𝐴3|𝑉)=𝑃(𝑉1𝐴3|𝑉) 𝑃(𝑉1𝐴3|𝑉)+𝑃(𝑉1 𝐴3|𝑉) =𝑃(𝑉1|𝑉)𝑃(𝐴3|𝑉1𝑉) 𝑃(𝑉1|𝑉)𝑃(𝐴3|𝑉1𝑉)+𝑃(𝑉1 |𝑉)𝑃(𝐴3|𝑉1 𝑉) Biologically speaking, regardless of mediators’ genealogical positions in G123 (depending on evolutionary events earlier than t123), “carriers” are still “carriers”, and “noncarriers” are still “noncarriers”. How g1, g2, and g3 are connected to mediators (depending on evolutionary events latter than t123) only depends on whether mediators are “carriers” or “noncarriers”. This is irrelevant to mediators’ genealogical positions of in G123. This essentially is a case of the Markov property: Given the “current state” (the states of mediators at t123), the “earlier state” is irrelevant to the “future state”. Accordingly, 𝑃(𝐴3|𝑉1𝑉)=𝑃(𝐴3|𝑉1 𝑉)=𝑃(𝐴3|𝑉), and the above expression becomes: 𝑃(𝑉1|𝐴3)=𝑃(𝑉1|𝑉) 𝑃(𝑉1|𝑉)+𝑃(𝑉1 |𝑉)=𝑃(𝑉1|𝑉) Following a similar logic, it can also be proven that under patterns A6, A7, and A8, P(W1 | Ax) = P(W1 | W) (x = 6, 7, 8). Supplementary Appendix S2 The equivalence relationship between pi and E(fi) can be more rigorously understood as follows: 𝑝𝑖=𝑃{𝑎𝑖 is a “carrier”} =𝑃{𝑔𝑖 is a descendant of a “carrier”} =∫ 𝑃(𝑔𝑖 is a descendant of a “carrier” | 𝑓𝑖)𝑝(𝑓𝑖) +∞ −∞ 𝑑𝑓𝑖 =∫ 𝑃(a descendant of a “carrier” is sampled from 𝑆𝑖 | 𝑓𝑖)𝑝(𝑓𝑖)𝑑𝑓𝑖 +∞ −∞ =∫ 𝑓𝑖𝑝(𝑓𝑖) +∞ −∞ 𝑑𝑓𝑖 =𝐸(𝑓𝑖) Where, fi denotes the frequency of descendants of “carriers” in species Si, p(fi) denotes the probability density function of fi. Additionally, it may be worth noting that the equivalence relationship between pi and E(fi) is unaffected by the presence of mutations occurring later than t123 (shallow mutations). This is because shallow mutations should be considered as “random noise” relative to mutation M itself. More specifically, 𝑃(𝑔𝑖 is a descendant of a “carrier” | 𝑓𝑖) is a probability over all possible carrying states of shallow mutations: 𝑃(𝑔𝑖 is a descendant of a “carrier” | 𝑓𝑖) =∑𝑃(𝑔𝑖 is a descendant of a “carrier” |𝑢,𝑓𝑖)𝑃(𝑢) ∞ Where, u denotes each possible carrying state of shallow mutations. Therefore, pi/E(fi) is a parameter that reflects the intrinsic property of mutation M. Supplementary Appendix S3 Here, we prove that under a neutral evolution assumption with a constant mutation rate and a constant population size, P(V1 | V) = 2/3. We use X1 and X2 to denote the waiting time for the first coalescent event and the waiting time for the second coalescent event, respectively (Fig. 9B). The waiting time for the first coalescent event X1 and the waiting time for the second coalescent event X2 follow an exponential distribution with rate parameter λ1 and an exponential distribution with rate parameter λ2, respectively. The expected waiting time for the first coalescent event and that for the second coalescent event are 1/λ1 and 1/λ2, respectively. Let R denote the ratio distribution of X2 to X1, i.e., R = X2/X1. As the first and the second coalescent events are mutually independent, the probability density function of R is as follows (Blom 1989): 𝑓(𝑟)={∫ 𝑥1∙𝜆1𝑒−𝜆1𝑥1∙𝜆2𝑒−𝜆2(𝑟𝑥1)𝑑𝑥1 ∞ 0, 𝑟>0 0, 𝑜𝑡ℎ𝑒𝑟𝑠 When r > 0, f(r) may be simplified as follows: 𝑓(𝑟)=∫ 𝑥1∙𝜆1𝑒−𝜆1𝑥1∙𝜆2𝑒−𝜆2(𝑟𝑥1)𝑑𝑥1 ∞ 0 =𝜆1𝜆2∫ 𝑥1∙𝑒−(𝜆1+𝜆2𝑟)𝑥1𝑑𝑥1 ∞ 0 =𝜆1𝜆2 (𝜆1+𝜆2𝑟)2Γ(2) Above, Γ(2) is a gamma function, which is equal to 1. Hence: 𝑓(𝑟)=𝜆1𝜆2 (𝜆1+𝜆2𝑟)2 =𝜆1𝜆2 ⁄ (𝜆1𝜆2 ⁄+𝑟)2 In a three-tip genealogical tree, the length of the first diverging branch is X1 + X2. The total length of the three terminal branches is 3X1 + X2. If mutation M falls on one of the three terminal branches, the three-tip genealogical tree is a type V genealogical tree. Among the three terminal branches, if mutation M falls on the first diverging branch, the three-tip genealogical tree is a type V1 genealogical tree (Fig. 9B). Therefore, let φ denote the conditional value of P(V1 | V) given that the values of X1 and X2 are X1 = x1 and X2 = x2, respectively. Assuming a constant mutation rate, φ is as follows: 𝜑= 𝑥1+𝑥2 3𝑥1+𝑥2=1+𝑟 3+𝑟 Then, the overall value of P(V1 | V) across all possible values of X1 and X2 may be written as follows: 𝑃(𝑉1|𝑉)=∫ 𝜑(𝑟)𝑓(𝑟)𝑑𝑟 ∞ 0 =∫ (1+𝑟 3+𝑟)[𝜆1𝜆2 ⁄ (𝜆1𝜆2 ⁄+𝑟)2]𝑑𝑟 ∞ 0 =(𝜆1𝜆2 ⁄ )∫(1− 2 3+𝑟)[1 (𝜆1𝜆2 ⁄+𝑟)2]𝑑𝑟 ∞ 0 =(𝜆1𝜆2 ⁄ )[∫ 1 (𝜆1𝜆2 ⁄+𝑟)2𝑑𝑟−∫ 2 (3+𝑟)(𝜆1𝜆2 ⁄+𝑟)2𝑑𝑟 ∞ 0 ∞ 0] Assuming a constant population size, the expected waiting time for the second coalescent event is 3 times longer than that for the first coalescent event (Kingman 1982; Tajima 1983), i.e., λ1/λ2 = 3. Hence: 𝑃(𝑉1|𝑉)=(𝜆1𝜆2 ⁄ )[∫ 1 (𝜆1𝜆2 ⁄+𝑟)2𝑑𝑟−∫ 2 (3+𝑟)(𝜆1𝜆2 ⁄+𝑟)2𝑑𝑟 ∞ 0 ∞ 0] =3×[∫ 1 (3+𝑟)2𝑑𝑟−∫ 2 (3+𝑟)3𝑑𝑟 ∞ 0 ∞ 0] =3×[−1 3+𝑟|∞ 0+1 (3+𝑟)2|∞ 0] =3×(1 3−1 9)=2 3 Supplementary Section S1: The Conclusion Regarding Gene Tree Distributions under Purifying Selection Remains Valid When Multiple Ancestral Mutations Are Assumed For demonstration, in the main text, we focus on the simplest one-mutation case. Here, we show that no matter how many mutations are assumed in ancestral species S123, the conclusions regarding mediators’ input and output remain valid. Let us assume n mutations occurring S123 (n can be any integer). Then, there are (𝑛 𝑛)+(𝑛 𝑛−1)⋯+(𝑛 0) possible allelic states. We consider an arbitrary pair of these allelic states, typex and typex+k, where typex denotes an allelic state with a combination of x mutations, and typex+k denotes an allelic state where k mutations are added to typex. If the mediators’ output is completely random, regardless of what allelic states are typex and typex+k, a typex gene copy should be equally likely to occupy each of the three tips of G123, and a typex+k gene copy should also be equally likely to occupy each tip of G123. The k mutations mentioned above can be viewed as a “super mutation”, and typex and typex+k gene copies can be viewed as “noncarriers” and “carriers”, respectively. Therefore, the effect of this “super mutation” can be understood in the same way as the effect of mutation M (Fig. 9). Above, there are two different allelic states among the three tips of G123. When more than one mutation is assumed in S123, it is also possible that there are three different allelic states among the three tips of G123. We now add another allelic state—typex+l. Accordingly, typex+k and typex+l gene copies are two types of “carriers”, and typex gene copies are “noncarriers”. If there are three different allelic states among the three tips of G123, there are two “super mutations” in G123—one is constituted by k mutations, the other is constituted by l mutations (see the figure below). If mediators’ output is completely random, the “noncarrier” in G123 should be equally likely to occupy each of the three tips. That is to say, the probability that the two “super mutations” fall on the first diverging branch should be equal to the probability that they fall on one of the two nested branches (see the figure below). This conflicts with the fact that the first diverging branch is longer than the two nested branches. Therefore, gene copies at t123 are not equally distributed among the tips of G123. Instead, it is expected that when an extant gene copy is connected to a “carrier”, its connection tends to be established with the first diverging branch of G123 (see the figure below). This effect acts in the same direction as the effect described in the main text. Therefore, we exclude the possibility that there is a pair of effects that cancel each other out. Collectively, when multiple mutations are assumed in S123, it remains biologically untenable to claim that changes in mediators’ input have no effect on gene tree distributions. Furthermore, when multiple mutations are assumed, how likely gi is connected to each type of gene copies at t123 (typex, typex-k, or typex-k-l) is still determined by the expected fates of these gene copies in Si. Therefore, regardless of how many mutations are assumed in S123, it can still be concluded that, under purifying selection, population size differences can lead to changes in mediators’ input, and consequently lead to changes in gene tree distributions. Illustration of the situations where the tips of G123 contain three types of gene copies. A. If mediators’ output is completely random, the type x gene copy (blue circle) should be equally likely to occupy each of the three possible positions. B. If mediators’ output is completely random, the two situations shown here should be equally likely to occur, which conflicts with the fact that the first diverging branch is longer than the two nested branches. The left situation is more likely to occur than the right one. Therefore, a “noncarrier” tends to bring the gene copy connected to it to a nested branch. 𝑃(𝑇1|𝐶3    )= ∑ 𝑃(𝑇1|𝐴𝑥)𝑃(𝐴𝑥|𝐶3    ) 𝑥=1,4,5,6 =1 3𝑞1𝑞2𝑞3+[1−𝑃(𝑉1|𝑉) 2]𝑞1𝑝2𝑞3+[1−𝑃(𝑉1|𝑉) 2]𝑝1𝑞2𝑞3+𝑃(𝑊1|𝑊)𝑝1𝑝2𝑞3 𝑞1𝑞2𝑞3+𝑞1𝑝2𝑞3+𝑝1𝑞2𝑞3+𝑝1𝑝2𝑞3 =1 3𝑞1𝑞2+[1−𝑃(𝑉1|𝑉) 2]𝑞1𝑝2+[1−𝑃(𝑉1|𝑉) 2]𝑝1𝑞2+𝑃(𝑊1|𝑊)𝑝1𝑝2 Again, we have shown that P(V1 | V) ≥ 2/3 and P(W1 | W) = 1. Hence, the above expression becomes: 𝑃(𝑇1|𝐶3    )≤1 3𝑞1𝑞2+1 6𝑞1𝑝2+1 6𝑝1𝑞2+𝑝1𝑝2 =1 3−1 6(𝑝1+𝑝2)+𝑝1𝑝2 =1 3+1 2𝑝1(𝑝2−1 3)+1 2𝑝2(𝑝1−1 3)(8) According to expression (8), if p1 < 1/3 and p2 < 1/3 (as shown above, almost always satisfied), 𝑃(𝑇1|𝐶3    )<1 3. Because 𝑃(𝑇1|𝐶3)>1 3 and 𝑃(𝑇1|𝐶3    )<1 3, it is expected that as p3 increases (i.e., g3 becomes more likely to be connected to a “carrier” at t123), T1 becomes more likely to occur. This conclusion is the same as that drawn by analyzing 𝜕𝑃(𝑇1) 𝜕𝑝3. In fact, the expression for 𝜕𝑃(𝑇1) 𝜕𝑝3 may be rearranged as follows: 𝜕𝑃(𝑇1) 𝜕𝑝3=−1 3𝑞1𝑞2+1 3𝑝1𝑝2 +𝑃(𝑉1|𝑉)𝑞1𝑞2−[1−𝑃(𝑉1|𝑉) 2]𝑞1𝑝2−[1−𝑃(𝑉1|𝑉) 2]𝑝1𝑞2 −𝑃(𝑊1|𝑊)𝑝1𝑝2+[1−𝑃(𝑊1|𝑊) 2]𝑝1𝑞2+[1−𝑃(𝑊1|𝑊) 2]𝑞1𝑝2 =1 3𝑝1𝑝2+𝑃(𝑉1|𝑉)𝑞1𝑞2+[1−𝑃(𝑊1|𝑊) 2]𝑝1𝑞2+[1−𝑃(𝑊1|𝑊) 2]𝑞1𝑝2 −1 3𝑞1𝑞2−[1−𝑃(𝑉1|𝑉) 2]𝑞1𝑝2−[1−𝑃(𝑉1|𝑉) 2]𝑝1𝑞2−𝑃(𝑊1|𝑊)𝑝1𝑝2 =𝑃(𝑇1|𝐶3)−𝑃(𝑇1|𝐶3    ) An increase in p3 can be understood as a process where a fraction of gene trees with a “noncarrier” a3 are replaced by a fraction of gene trees with a “carrier” a3. In this fraction of gene trees, as gene trees with a "noncarrier" a3 are replaced by gene trees with a “carrier” a3, the probability of occurrence of gene tree topology T1 will change from 𝑃(𝑇1|𝐶3    ) to 𝑃(𝑇1|𝐶3). When this gene tree replacement process is expressed in mathematical language, it is the expression presented above. Therefore, the analysis of 𝜕𝑃(𝑇1) 𝜕𝑝3 and the analysis of 𝑃(𝑇1|𝐶3    ) and 𝑃(𝑇1|𝐶3) are essentially equivalent. Clarifying the equivalence between the analysis of 𝜕𝑃(𝑇1) 𝜕𝑝3 and the analysis of 𝑃(𝑇1|𝐶3    ) and 𝑃(𝑇1|𝐶3) help better understand the biological meaning of the expression for 𝜕𝑃(𝑇1) 𝜕𝑝3. 𝜕𝑃(𝑇2) 𝜕𝑝3 and 𝜕𝑃(𝑇3) 𝜕𝑝3 can be similarly understood by decomposing them into 𝑃(𝑇2|𝐶3    ) and 𝑃(𝑇2|𝐶3), and into 𝑃(𝑇3|𝐶3    ) and 𝑃(𝑇3|𝐶3), respectively. Reference Blom G. 1989. Ratio of random variables. In: Fienberg S, Olkin I, editors. Probability and Statistics, Theory and Applications. New York, NY: Springer-Verlag. p. 91–92. Kingman JFC. 1982. The coalescent. Stoch. Process. their Appl. 13:235–248. Tajima F. 1983. Evolutionary relationship of DNA sequences in finite populations. Genetics 105:437–460. Supplementary Figure S1. A. A comparison of simulation results of different levels of scaling. The result in the left side was generated with parameters identical to that shown in Figure 5F, the result in the right side was generated using 5-fold larger population sizes, a 5-fold lower selection coefficient, and a 5-fold lower mutation rate (i.e. using a scaling factor that is 1/5 of that of Fig. 5F). These results are similar. Each simulation result shown here only contains 3000 × 20 = 60,000 gene trees (in Fig. 5, each contains 6000 × 20 = 120,000). Therefore, the variation among repliates is greater than that in Figure 5. Figure S1 Identical to Fig. 5F (N123 = 200, μ = 1.5 × 10-5, and s = -0.0075) 5-fold larger population sizes (N123 = 1000, μ = 3 × 10-6, and s = -0.0015) N1=(1/25)×N12 N2=N12 N3=N12 Population sizes used by Vanderpool et al. (2020) N12: N123: N1: N2: N3 = 25:25:1:25:25 (σ123 = -7.5) Population sizes used by He et al. (2020) N12: N123: N1: N2: N3 = 5:5:1:25:25 (σ123 = -1.5) Supplementary Figure S2. A. The species tree set in the simulations of Vanderpool et al. (2020). The left side shows the SLiM codes used by Vanderpool et al. (2020). These codes correspond to the species tree shown on the right side. This species tree is not comparable to any of those used by He et al. (2020). Please compare the species tree shown here and that shown in Figure 5. Additionally, each locus in Vanderpool et al.’s simulations has 41000 possible allelic states, whereas each locus in the simulations of He et al. (2020) is biallelic. Because of these differences, Vanderpool et al.’s simulations cannot be used to examine the replicability of the simulation results of He et al. (2020). Among the differences mentioned above, the most crucial difference lies in N123. The population sizes used by Vanderpool et al. (2020) are in the ratio N12:N123:N1:N2:N3 = 25:25:1:25:25, whereas the population sizes used by He et al. (2020) are in the ratio N12: N123: N1: N2: N3 = 5:5:1:25:25. That is to say, the N123 used by Vanderpool et al. (2020) is effectively 5 times greater than that used by He et al. (2020). Consequently, the scaled selection coefficient in ancestral species S123 (σ123 ) in Vanderpool et al.’s simulations has an absolute value 5 times greater than that in the simulations of He et al. (2020) (7.5 versus 1.5). Because of this difference, in Vanderpool et al.’s simulations, deleterious mutations occurring in ancestral species S123 have a much lower probability of being present in extant species S1, S2, and S3. In such a situation, the ancestries of g1, g2, and g3 are often traced back to three identical ancestors. It is unsurprising that asymmetry in gene tree distribution is not pronounced in such a situation. In reality, sligtly deleterious mutations have a chance of being present in species S1, S2, and S3. Therefore, the simulation settings used by Vanderpool et al. (2020) are not realistic. B. We changed the population sizes used by Vanderpool et al. (2020) to those used by He et al. (2020). Other settings used by Vanderpool et al. (2020) remained unchanged. This change in population sizes casues a change in the distribution of gene tree topologies. Such a phenomenon contradicts Vanderpool et al.’s claim “there should be no effect of negative selection on the distribution of tree topologies”. N1N2N3 0.005×(2N12) 4×(2N12) 10×(2N12) N12 N123 Species tree used by Vanderpool et al. (2020) Population sizes used by Vanderpool et al. (2020) N12: N123: N1: N2: N3 = 25:25:1:25:25 (σ123 = 0 ) Neutral Purifying selection Purifying selection A B N123=N12 Figure S2 Supplementary Figure S3. The values of P(V1|V) under neutral evolution and under purifying selection. Under neutral evolution, the observed values of P(V1|V) agree with the theoretical expectation described in the main text, P(V1|V) = 2/3. Under purifying selection, the observed values of P(V1|V) are slightly greater than 2/3. The result of the Mann–Whiteney U test suggests that the difference between the values of P(V1|V) under neutral evolution and those under purifying selection is statistically significant (p = 5.2 × 10-5). This phenomenon can be understood as follows: under purifying selection, “noncarriers” tend to produce more offspring in comparison with under neutral evolution. The necessary condition for the occurrence of V1 genealogical trees is that a “noncarrier” produce more than one offspring. Therefore, under purifying selection, it is expected that V1 genealogical trees have a greater probability of occurrence in comparison with under neutral evolution. In addition, this phenomenon can be understood by the ancestral selection graph (ASG) model. In the ASG model, the influence of selection is represented by “incoming branches”. The presence of incoming branches can cause a proportion of type W1 genealogical trees to be transformed into type V1 genealogical trees. Therefore, the value of P(V1|V) is expected to increase under purifying selection. Figure S3 Neutral (σ123 = 0) Purifying selection (σ123 = -1.5) P(V1|V)