==== Front Bioinformatics Bioinformatics bioinformatics Bioinformatics 1367-4803 1367-4811 Oxford University Press 37387175 10.1093/bioinformatics/btad263 btad263 Evolutionary, Comparative and Population Genomics AcademicSubjects/SCI01060 Phylogenetic diversity statistics for all clades in a phylogeny Grover Siddhant Department of Computer Science, Iowa State University, Ames, IA 50010, United States Markin Alexey Virus and Prion Research Unit, National Animal Disease Center, USDA-ARS, Ames, IA 50010, United States https://orcid.org/0000-0002-3138-5535 Anderson Tavis K Virus and Prion Research Unit, National Animal Disease Center, USDA-ARS, Ames, IA 50010, United States Eulenstein Oliver Department of Computer Science, Iowa State University, Ames, IA 50010, United States Corresponding authors. Department of Computer Science, Iowa State University, Ames, IA 50010, United States. E-mails: grover16@iastate.edu (S.G.) and oeulenst@iastate.edu (O.E.) 6 2023 30 6 2023 30 6 2023 39 Suppl 1 ISMB/ECCB 2023 Proceedings i177i184 © The Author(s) 2023. Published by Oxford University Press. 2023 https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract The classic quantitative measure of phylogenetic diversity (PD) has been used to address problems in conservation biology, microbial ecology, and evolutionary biology. PD is the minimum total length of the branches in a phylogeny required to cover a specified set of taxa on the phylogeny. A general goal in the application of PD has been identifying a set of taxa of size k that maximize PD on a given phylogeny; this has been mirrored in active research to develop efficient algorithms for the problem. Other descriptive statistics, such as the minimum PD, average PD, and standard deviation of PD, can provide invaluable insight into the distribution of PD across a phylogeny (relative to a fixed value of k). However, there has been limited or no research on computing these statistics, especially when required for each clade in a phylogeny, enabling direct comparisons of PD between clades. We introduce efficient algorithms for computing PD and the associated descriptive statistics for a given phylogeny and each of its clades. In simulation studies, we demonstrate the ability of our algorithms to analyze large-scale phylogenies with applications in ecology and evolutionary biology. The software is available at https://github.com/flu-crew/PD_stats. Department of Agriculture 10.13039/100000199 Agricultural Research Service 10.13039/100007917 5030-32000-231-000-D 5030-32000-231-095-S National Institute of Allergy and Infectious Diseases 10.13039/100000060 National Institutes of Health 10.13039/100000002 Department of Health and Human Services 10.13039/100000016 75N93021C00015 USDA Agricultural Research Service 10.13039/100007917 0500-00093-001-00-D ==== Body pmc1 Introduction Habitat loss, biological invasions, emerging infectious diseases, and climate change are altering global ecological systems. Extinctions have vastly outpaced speciation events over the past century (Barnosky et al. 2011, Ceballos et al. 2015), and novel pathogens are emerging at an increasing rate with two pandemics within the past 20 years (Gibb et al. 2020, Keusch et al. 2022). Quantifying biodiversity is consequential as it permits assessment of when and where significant biological changes are occurring. Traditional approaches to measuring changes in diversity patterns have taken the form of generating phylogenetic trees from genetic sequence data and then extracting measures of phylogenetic diversity (PD) from the inferred trees (Cadotte et al. 2010). A now classic measure to quantify diversity is the PD index introduced by Faith (1992). PD measures the phylogenetic history among taxa occurring in a given sample and is calculated on a rooted phylogenetic tree for a set of taxa in the tree (i.e. represented as leaves). It is defined as the sum of the edge lengths that span these taxa along with the root of the tree (Fig. 1). This index directly interprets the evolutionary history of the taxon set and has been applied to generate solutions to practical problems in ecology and evolutionary biology (Li et al. 2020). Consequently, there has been extensive computational work considering the mathematical and combinatorial properties of the index (Steel 2005, Hartmann and Steel 2006, Minh et al. 2006, Moulton et al. 2007, Bordewich and Semple 2012). Specifically, there are polynomial time algorithms to solve the maximum PD problem on phylogenetic trees (Halldórsson et al. 1999, Steel, 2005, Minh et al. 2006, Hartmann and Steel 2006). However, PD applications (e.g. Faith and Richards 2012, Faith 2018, Faith and Simon 2018) could benefit from understanding descriptive statistics associated with PD across the entire phylogeny in addition to the maximum PD and the maximal PD set. These statistics allow for the identification of critical nodes within the phylogeny and for comparison of diversity statistics between nodes. Figure 1. An example of a phylogenetic tree. The internal nodes below the root are labeled by roman numerals. The phylogenetic diversity of a set of taxa {A, B, D} is 3+4+1+1+2+1=12 as can be seen from the highlighted (bold) edges. There has been progress toward calculating the minimum, average, and variance of PD on a phylogeny with n leaves and a subset of those leaves of size k. Hartmann et al. proposed using the minimum PD value to compute the gain in PD for the given set of taxa when compared with the worst-case scenario (i.e. a set of taxa with the smallest PD) (Hartmann and André 2013). Similarly, Manson et al. suggested using minimum PD to find the worst-case conservation scenario — i.e. find k taxa whose extinction will decrease the overall PD the most (Manson et al. 2022). Both Hartmann and André (2013) and Manson et al. (2022) provided algorithms for finding a value and a set for minimum PD. However, the complexity of the algorithm was not analyzed by Hartmann and André (2013), and the algorithm described by Manson et al. (2022) is prohibitive for larger phylogenies (and larger k) due to its O(nk3) run time. The average PD was also explored by Hartmann et al. to compute how “beneficial” a taxon set is when compared with a randomly chosen set of k taxa (Hartmann and André 2013), i.e. to identify sets of taxa with higher than average PD. Tsirogiannis et al. (2012, 2014) proposed algorithms to compute both average PD and variance PD. Though these approaches computed the PD statistics approximately in linear time, they do not have an approximation guarantee as they approximate ratios of binomial coefficients using the hypergeometric distribution. In this work, we introduce dynamic programming algorithms designed to efficiently and accurately compute the distribution of PD on a phylogeny by calculating the associated descriptive statistics for each clade on a given tree. Our solution for the Minimum PD problem has O(kn) time-complexity, significantly improving on the previous algorithm by Manson et al. (2022) that runs in O(nk3) time, where n is the number of taxa of a phylogeny and k is the user-specified size of a PD set (query size). We achieved this complexity bound by applying the technique from Halldórsson et al. to our dynamic programming formulation (Tamir 1996, Halldórsson et al. 1999). We extended our dynamic programming framework to the Average PD and Variance PD problems, obtaining precise algorithms that run in O(k2n log2n) time. We designed exact algorithms using big-number arithmetics with binomial coefficients. Crucially, our algorithms compute PD statistics (minimum, maximum, average, and variance) not only for the overall tree but also for every clade in the tree. Moreover, given a fixed input value k, the algorithms compute PD statistics for all 2≤ki≤k. Using a simulation study, we demonstrated that our dynamic programming framework can process large-scale trees with 10,000 taxa in under an hour on a standard PC, computing all of the PD statistics for each clade and every taxon size k, where 1|v|, the corresponding values are undefined). 3 Solving PD Problems We introduce efficient dynamic programing algorithms to solve the PD problems for all clades in a given tree with n leaves and an integer k. First, we prove a O(kn) recurrence for the minimum and maximum PD problems. Then, for the first time, we describe efficient algorithms for exactly solving the average PD and variance PD problems in O(k2n log2n) time. 3.1 Minimum PD problem We define δ(v,k)-subsolutions at subtrees and describe an O(kn) time algorithm to compute δ(k) for each clade in T with taxon set X. In a tree T, the subproblem at v∈V(T), for the MinPD problem, is defined as the problem to compute δ(k) for Tv. Thus, for 1≤k≤|v| in a tree rooted at a node v, we define δ(v,k) to be the min PDTv(S). That is, and δ(ρ,k)=δ(k).Proposition 1For every node v∈V(T)with children x and y (if they exist), and 1≤k≤|v|:δ(v,k)=0, if v∈X, elsewhere w(x↓)is the weight of the edge from v to x, and w(y↓)is the weight of the edge from v to y. Proof. For v being a leaf and any k (the base case), δ(v,k)=0, since there is no diversity to be selected from the subtree rooted at a leaf. For an internal node v with children x and y, there are multiple possibilities for splitting k taxa between the subtrees rooted at x and y. δ(v,k):=minS∈vkPDTv(S), δ(v,k)=min{w(x↓)+δ(x,k),w(y↓)+δ(y,k),w(x↓)+w(y↓)+minr+l=k1≤r≤|x|1≤l≤|y|(δ(x,r)+δ(y,l)) Case I: all the k taxa in a minimal PD set S are from the subtree rooted at x. Note that for any S∈xk, we have PDTv(S)=PDTx(S)+w(x↓). Therefore, set S that minimizes PDTx(S) also minimizes PDTv(S). That is, δ(v,k)=δ(x,k)+w(x↓). Case II: all the k taxa in a minimal PD set S are from the subtree rooted at y. Similarly to Case I, in this case δ(v,k)=δ(y,k)+w(y↓). Case III: k taxa are split between the subtrees rooted at x and y. Let us fix some k1 and k2, such that 1≤k1≤|x|, 1≤k2≤|y|, k1+k2=k. Then for S1∈xk1 and S2∈yk2, we have PDTv(S1∪S2)=PDTx(S1)+PDTy(S2)+w(x↓)+w(y↓). Then minimizing over all such sets S1 and S2 reduces to independent minimization in the left and right subtrees. It is then not difficult to see that minS∈vkPD(S)=w(x↓)+w(y↓)+mink1+k2=k1≤k1≤|x|1≤k2≤|y|(δ(x,k1)+δ(y,k2)). Combining all three above cases, we obtain the formula from Proposition 1. □ Our pseudocode for the algorithm using the recurrence above is presented in Algorithm 1. Algorithm 1 MinPD Problem 1: procedureComputeMin(T, k) 2: δ←[1;2n−1][0;k] 3:    forv∈Xdo 4:      for i = 0 to k do 5:        δ[v,i]=0 6:      end for 7:    end for 8:    forv∈V∖X(in post order) do 9:     //x, yare children ofv 10:      for i = 1 to min(k,|v|)do 11:        minval=∞ 12:        forr=max(0,i−|y|) to min(|x|,i)do 13:          l←i−r 14:          val = δ[x,r]+w(x↓)·min(r,1) 15:           + δ[y,l]+w(y↓)·min(l,1) 16:          ifval < minval 17:           minval = val 18:      end for 19:      δ[v,k]=val 20:     end for 21:    end for 22: end procedure Theorem 2 Algorithm 1 solves the MinPD problem in O(kn) time. Proof. Algorithm 1 is built on the recurrence relation proved in Proposition 1. The complexity of Algorithm 1 is as follows. Line 2 creates an empty 2D array of O(kn) space. The outer loop in line 3 runs in O(n) time, while the inner loop in line 4 runs in O(k) time; thus the runtime for lines 3−7 is O(kn). At first, it is clear that lines 8−20 take O(k2n) time; however, Tamir (1996) and Halldórsson et al. (1999) proved that dynamic programming of this type in fact runs in O(kn) time. We explicitly summarize their results in Lemma 3. Lemma 3 (Tamir 1996, Halldórsson et al. 1999) For a non-leaf node v with children x and y, let Iv denote the number of pairs (l, r) such that 0≤l≤|x|, 0≤r≤|y|, and l+r≤min(k,|v|). Then, I(T):=∑v∈V(T)Iv∈O(kn) It is not difficult to see that I(T) bounds the runtime of Algorithm 1. Thus, MinPD for every clade can be computed in O(kn) time and space. Once δ(k) is obtained, the Minimal PD set can be found via backtracking in O(k+n) time. □ 3.2 Maximum PD problem The original MaxPD problem can be solved in O(n) time (Spillner et al. 2008). Then, solving this problem for every clade in a tree would require O(n2) time. However, adapting the dynamic programming approach from the previous section, we can solve the same problem in O(kn) time, which is more desirable, especially for smaller k. The adaptation would only require changing the minimization criterion in Algorithm 1 to maximization. Additionally, it might be possible to extend the algorithm from Spillner et al. (2008) to achieve a similar runtime. 3.3 Average PD problem We begin by defining AvgPD subproblems and then prove the recurrence relation for computing the β(v,k)-subsolution at Tv (defined below). We then present an O(k2n log2n) time and O(k2n log n) space algorithm to compute α(k) for each clade in T. We define β(k):=∑S∈ρkPD(S). That is, β(k)=α(k)·(|X|k). We then define the subproblem β(v,k) at a vertex v∈V(T), for the AvgPD problem, as the problem to compute β(k) in Tv subtree. That is, for 1≤k≤|v|, β(v,k):=∑S∈vkPDTv(S). Then, β(k)=β(ρ,k) and α(k)=β(k)/(|X|k). Proposition 4 For every node v∈V(T)with children x and y (if they exist), and 1≤k≤|v|: β(v,k)=0, if v∈X, else β(v,k)=[∑r+l=k1≤r≤|x|1≤l≤|y|(|y|l)β(x,r)+(|x|r)β(y,l)+(|x|r)(|y|l)(w(x↓)+w(y↓))]+[β(x,k)+(|x|k)w(x↓)]+[β(y,k)+(|y|k)w(y↓)] Proof. For v∈X, β(v,1)=PDTv({v})=0. We now prove the recurrence for internal v with children x and y. We split β(v,k)=∑S∈vkPDTv(S) into three categories. Case I : all the k taxa in S are chosen from the subtree rooted at x, i.e. S∈xk. Then, βI(v,k)=∑S∈xk(PD(S)+w(x↓))=β(x,k)+(|x|k)w(x↓). Case II : all the k taxa in S are chosen from the subtree rooted at y. Similarly to Case I, we have βII(v,k)=β(y,k)+(|y|k)w(y↓). Case III : k taxa are split between subtrees rooted at x and y. Let k1 and k2 be such that 1≤k1≤|x|, 1≤k2≤|y|, k1+k2=k. Then we define β(v,k1,k2):=∑S1∈xk1S2∈yk2PDTv(S1∪S2). It is then not difficult to see that βIII(v,k)=∑k1+k2=k1≤k1≤|x|1≤k2≤|y|β(v,k1,k2). Next, observe the following β(v,k1,k2)=∑S1∈xk1S2∈yk2(PDTx(S1)+PDTy(S2)+w(x↓)+w(y↓))=∑S1∈xk1((|y|k2)(PDTx(S1)+w(x↓)+w(y↓))+β(y,k2))=(|x|k1)(|y|k2)(w(x↓)+w(y↓)))+(|y|k2)β(x,k1)+(|x|k1)β(y,k2) Then, observing gives us the equation in Proposition 4. □ β(v,k)=βI(v,k)+βII(v,k)+βIII(v,k) Finally, we have where α(v,k) is the solution for the the problem to compute α(k) in Tv. Our pseudo code for the algorithm using the recurrence above is presented above. α(v,k)=((|v|k)−1)·β(v,k) Algorithm 2 AvgPD Algorithm 1: procedureComputeAvg(T, k) 2: Precompute (ij) values for i≤n,j≤k via a Pascal triangle. 3:   α←[1;2n−1][1;k] 4:   β←[1;2n−1][1;k] 5:   forv∈Xdo 6:     α[v,1]=0 7:     β[v,1]=0 8:   end for 9:   forv∈V∖X(in post order) do 10:    //x, yare children ofv 11:     for i = 1 to min(k,|v|)do 12:       sum_x = β[x,i]+(|x|i)w(x↓) 13:       sum_y = β[y,i]+(|y|i)w(y↓) 14:       sum = sum_x + sum_y 15:       forr=max(1,i−|y|) to min(|x|,i−1)do 16:         l←i−r 17:         sum+=(|y|l)β[x,r]+(|x|r)β[y,l] 18:            +(|x|r)(|y|l)(w(x↓)+w(y↓)) 19:       end for 20:       β[v,i]=sum 21:       α[v,i]=sum/(|v|i) 22:     end for 23:    end for 24: end procedure Theorem 5 Algorithm 2 solves the AvgPD in O(k2n log 2n)time and O(k2n log n)space. Proof. The analysis of the complexity of Algorithm 2 is similar to the analysis of Algorithm 1. The main difference is the fact that (nk) values can be very large, hence we need to account for the arithmetic on large numbers. Note that (nk)≤nk; consequently, the number of bits required to encode such numbers is O(log(nk))=O(k log n). Note that precomputing all values (ij) for i≤n and j≤k can be done by computing a part of a Pascal triangle in O(kn·k log n) time. The remainder of the algorithm requires big-number arithmetic. Note that multiplication of two O(k log n)-bit numbers can be done in O(k logn log k)=O(k log 2n) time (Harvey and Van Der Hoeven 2021). Therefore, the remainder of the algorithm requires O(kn·k log 2n) time, where kn follows from the proof of Theorem 2. □ 3.4 Variance PD problem Recall that by definition, ψ(k):=∑S∈ρk(PD(S)−α(k))2(|X|k)=((|X|k)−1)·(∑S∈ρkPD(S)2+∑S∈ρkα(k)2−2∑S∈ρkPD(S)·α(k))=∑S∈ρkPD(S)2+(|X|k)α(k)2−2(|X|k)α(k)2(|X|k)=∑S∈ρkPD(S)2(|X|k)−α(k)2 Using Algorithm 2, we can compute α(k) in O(k2n log n) time. Therefore, we are only left to compute the sum of squares ∑S∈ρkPD(S)2, which we denote by γ(k). For convenience, we also define for any v∈V(T) and 1≤k≤|v|. γ(v,k):=∑S∈vkPDTv(S)2, Proposition 6 For every node v∈V(T)with children x and y (if exist), and 1≤k≤|v|, γ(v,k)=0, if v∈X, elseγ(v,k)=∑1≤r,l≤k−1r+l=k[(|y|l)γ(x,r)+(|x|r)γ(y,l)+(|x|r)(|y|l)λ2+2(β(x,r)β(y,l)+λ((|y|l)β(x,r)+(|x|r)β(y,l)))]+[γ(x,k)+(|x|k)w(x↓)2+2w(x↓)β(x,k)]+[γ(y,k)+(|y|k)w(y↓)2+2w(y↓)β(y,k)] where λ represents w(x↓)+w(y↓)and β(v,k)=∑S∈vkPDTv(S). Proof. For a leaf v∈X, γ(v,1)=PDTv({v})2=0. We now prove the recurrence for an internal vertex v with children x and y. We split γ(v,k) into three categories. Case I: all the k taxa in S are chosen from the subtree rooted at x, i.e. S∈xk. Then, γI(v,k)=∑S∈xk(PDTx(S)+w(x↓))2=γ(x,k)+(|x|k)w(x↓)2+2w(x↓)β(x,k). Case II: all the k taxa in S are chosen from the subtree rooted at y. Similarly to Case I, we have γII(v,k)=γ(y,k)+(|y|k)w(y↓)2+2w(y↓)β(y,k), Case III : k taxa are split between the subtrees rooted at x and y. Let k1 and k2 be such that 1≤k1≤|x|, 1≤k2≤|y|, k1+k2=k. For convenience, we assign λ:=w(x↓)+w(y↓). Then, we define γ(v,k1,k2):=∑S1∈xk1S2∈yk2PDTv(S1∪S2)2 Observe now that γ(v,k1,k2)=∑S1∈xk1S2∈yk2(PDTx(S1)+PDTy(S2)+λ)2=∑S1∈xk1S2∈yk2(PDTx(S1)2+PDTy(S2)2+λ2++2λ(PDTx(S1)+PDTy(S2))++2PDTx(S1)PDTy(S2))=(|y|k2)γ(x,k1)+(|x|k1)γ(y,k2)+(|x|k1)(|y|k2)λ2+2λ((|y|k2)β(x,k1)+(|x|k1)β(y,k2))+2β(x,k1)β(y,k2) and γIII(v,k)=∑k1+k2=k1≤k1≤|x|1≤k2≤|y|γ(v,k1,k2). Proposition 7 then follows from the observation that γ(v,k)=γI(v,k)+γII(v,k)+γIII(v,k). □ Theorem 7 Algorithm 3 solves the VarPD in in O(k2n log 2n)time and O(k2n log n)space. Proof. The proof is similar to the proof of Theorem 5. □ Algorithm 3 VarPD Algorithm 1: procedureComputeVarPD(T, k) 2:   Compute (ij) values for i≤n,j≤k via a Pascal triangle. 3:   Compute α and β arrays using Algorithm 2. 4:   ψ←[1;2n−1][1;k]; γ←[1;2n−1][1;k] 5:   forv∈Xdo 6:     for i = 1 to k do 7:       ψ[v,i]=γ[v,i]=0 8:     end for 9:   end for 10:   forv∈V∖X (in post order) do 11:    //x, yare children ofv 12:     for i = 1 to min(k,|v|)do 13:   sq_x= γ(x,k)+(|x|k)w(x↓)2+2w(x↓)β[x,k] 14:   sq_y= γ(y,k)+(|y|k)w(y↓)2+2w(y↓)β[y,k] 15:   sq = sq_x + sq_y 16:   forr=max(0,i−|y|) to min(|x|,i)do 17:         sq += (|y|l)γ[x,r]+(|x|r)γ[y,l] 18:         +(|x|r)(|y|l)λ2+2[β[x,r]β[y,l] 19:         +λ((|y|l)β[x,r]+(|x|r)β[y,l])] 20:   end for 21:       γ[v,i]=sq; ψ[v,i]=γ[v,i](|v|i)-α[v,i]2 22:     end for 23:   end for 24: end procedure 4 Experimental evaluation We demonstrate the applicability of the PD statistics algorithms by evaluating their scalability. Additionally, we compare our exact algorithms for calculating AvgPD and VarPD to calculating these measures through randomly sampling taxon sets. The algorithms were implemented in Python 3 and executed on a PC with an Intel Core i7 3.19 Ghz processor using 8.0 GB of RAM under the Windows 10 operating system. Dendropy (Sukamaran and Holder 2010) and Biopython (Cock et al. 2009) were used in our implementation. 4.1 Scalability analysis To demonstrate the applicability of our algorithms, we generated large-scale phylogenetic trees using the birth/death model (1.0 birth rate and 0.5 death rate) with the number n of taxa ranging from 1000 to 10,000 and a step of 1000. For each fixed number of taxa, we generated 10 trees. The PD set size k was set to match the total number of taxa n to measure the worst-case complexity of our algorithms. Thus, we compute the PD statistics for all clades in a tree and all k ranging from 2 to n−1. Figure 2 demonstrates that for a 10,000 taxon tree, all of the PD statistics on all clades (and all k) can be computed within an hour. Figure 2. Runtime of the algorithms to compute PD statistics on trees with 1000 – 10,000 taxa for k = n. Since our algorithms for AvgPD and VarPD problems depend on the preprocessing step for computing the binomial coefficients, we report the runtime of those algorithms separately from the preprocessing runtime. 4.2 AvgPD and VarPD algorithms are significantly faster than random sampling When exact algorithms for computing the mean/variance of a distribution are not known, random sampling is typically used to approximate those moments. Here we compare the runtime of our methods for AvgPD and VarPD problems against a random sampling strategy. We simulated a birth/death tree (1.0 birth rate and 0.5 death rate) with n = 500 taxa and tested different k between 50 and 450 with a step of 50. After computing the exact AvgPD and VarPD values using our algorithms, we sampled subsets of taxa of size k uniformly at random and computed the PD of the chosen samples. After every 200 samples, we computed the overall mean/variance PD of the samples until that value converged to the truth within 0.1%. Below, we report the average runtimes of the algorithms over various k. For the AvgPD problem, random sampling required 14.2 s to converge, while our algorithm required 2 s to compute, on average. The difference was more significant for the VarPD problem, where random sampling required 28 min to converge (on average), while our algorithm required only 2.8 s. Further, for AvgPD, the gap between our algorithm and random sampling increases with larger n. For example, for n = 2000 and k = 20, random sampling required >1 h to converge, while our algorithm required only 6.5 s to complete. 5 Conclusion In this work, we devised effective, efficient algorithms for computing the PD statistics. Notably, our algorithms are not restricted to calculating these values for the entire tree; but for every subtree/clade present in T. Furthermore, our algorithms provide descriptive statistics of PD, including the minimum PD, maximum PD, average PD, and variance PD for all values ranging from 1 up to k. A major merit of our algorithms for the minimum and maximum PD is that these values can be computed simultaneously by building on a similar dynamic programming approach. All the algorithms described in this work are efficient and scale well in practice. Consequently, these algorithms can be applied in the field without limitation based on the number of taxa or k. Quantifying statistics such as the average PD and variance PD across all nodes in a phylogeny can provide deep insights into evolutionary histories. The average PD is the mean PD value over all sets of k taxa on a tree. In other words, it is the expected PD value if taxa were chosen uniformly at random. Similarly, the variance PD is the variance of the PD values given an equiprobable distribution of selecting any k taxa. Standardizing these diversity measures is crucial, making values comparable for different values of k and across different phylogenies. Further, when computed for all taxa in a phylogeny, these statistics can be used to assess the role of evolutionary history in shaping ecological communities, to determine how this evolutionary history influences niche and resource use, and can be used to study trait evolution and biogeography (Webb et al. 2002). Acknowledgements Mention of trade names or commercial products in this article is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the USDA. USDA is an equal opportunity provider and employer. Conflict of interest: The authors declare no conflicting interests. Funding This work was supported in part by the U.S. Department of Agriculture (USDA) Agricultural Research Service [ARS project number 5030-32000-231-000-D; ARS project number 5030-32000-231-095-S]; the National Institute of Allergy and Infectious Diseases, National Institutes of Health, Department of Health and Human Services [contract number 75N93021C00015]; and the SCINet project of the USDA Agricultural Research Service [ARS project number 0500-00093-001-00-D]. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. Data availability The data underlying this article are available in P D_stats, at https://github.com/flucrew/PD_stats/. ==== Refs References Barnosky AD , MatzkeN, TomiyaS et al Has the Earth’s sixth mass extinction already arrived? Nature 2011;471 :51–7.21368823 Bordewich M , SempleC. Budgeted nature reserve selection with diversity feature loss and arbitrary split systems. J Math Biol 2012;64 :69–85.21301849 Cadotte MW , Jonathan DaviesT, RegetzJ et al Phylogenetic diversity metrics for ecological communities: integrating species richness, abundance and evolutionary history. Ecol Lett 2010;13 :96–105.19903196 Ceballos G , EhrlichPR, BarnoskyAD et al Accelerated modern human–induced species losses: entering the sixth mass extinction. Sci Adv 2015;1 :e1400253.26601195 Cock PJ , AntaoT, ChangJT et al Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics 2009;25 :1422–3.19304878 Faith DP. Conservation evaluation and phylogenetic diversity. Biol Conserv 1992;61 :1–10. Faith DP. Phylogenetic diversity and conservation evaluation: perspectives on multiple values, indices, and scales of application. Phylogenet Divers 2018;1–26. Faith DP , RichardsZT. Climate change impacts on the tree of life: changes in phylogenetic diversity illustrated for Acropora corals. Biology (Basel) 2012;1 :906–32.24832524 Faith DP , VeronS, PavoineS, et al Indicators for the expected loss of phylogenetic diversity. Phylogenet Divers 2018;:73–91. Gibb R , ReddingDW, ChinKQ et al Zoonotic host diversity increases in human-dominated ecosystems. Nature 2020;584 :398–402.32759999 Halldórsson MM , IwanoK, KatohN et al Finding subsets maximizing minimum structures. SIAM J Discrete Math 1999;12 :342–59. Hartmann K , AndréJ. Should evolutionary history guide conservation? Biodivers Conserv 2013;22 :449–58. Hartmann K , SteelM. Maximizing phylogenetic diversity in biodiversity conservation: greedy solutions to the Noah’s Ark problem. Syst Biol 2006;55 :644–51.16969940 Harvey D , Van Der HoevenJ. Integer multiplication in time O (n log 2). Ann. of Math 2021;193 :563–617. Keusch GT , AmuasiJH, AndersonDE et al Pandemic origins and a one health approach to preparedness and prevention: solutions based on SARS-CoV-2 and other RNA viruses. Proc Natl Acad Sci U S A 2022;119 :e2202871119.36215506 Li D , OldenJD, LockwoodJL et al Changes in taxonomic and phylogenetic diversity in the Anthropocene. Proc Biol Sci 2020;287 :20200777.32546087 Manson K , SempleC, SteelM. Counting and optimising maximum phylogenetic diversity sets. J Math Biol 2022;85 :1–23.35767083 Minh BQ , KlaereS, von HaeselerA. Phylogenetic diversity within seconds. Syst Biol 2006;55 :769–73.17060198 Moulton V , SempleC, SteelM. Optimizing phylogenetic diversity under constraints. J Theor Biol 2007;246 :186–94.17275037 Spillner A , NguyenBT, MoultonV. Computing phylogenetic diversity for split systems. IEEE/ACM Trans Comput Biol Bioinform 2008;5 :235–44.18451432 Steel M. Phylogenetic diversity and the greedy algorithm. Syst Biol 2005;54 :527–9.16051588 Sukumaran J , HolderMT. DendroPy: a Python library for phylogenetic computing. Bioinformatics 2010;26 :1569–71.20421198 Tamir A. An O(pn2) algorithm for the p-median and related problems on tree graphs. Operat Res Lett 1996;19 :59–64. Tsirogiannis C , BrodyS, DimitrisC. Efficient computation of popular phylogenetic tree measures. In: International Workshop on Algorithms in Bioinformatics, Wroclaw, Poland. 2012, pp. 30–43. Tsirogiannis C , BrodyS, AdrijaK. New algorithms for computing phylogenetic biodiversity. In: International Workshop on Algorithms in Bioinformatics, Ljubljana, Slovenia. 2014, pp. 187–203. Webb CO , AckerlyDD, McPeekMA et al Phylogenies and community ecology. Annu Rev Ecol Syst 2002;33 :475–505.