Fast Computations for Measures of Phylogenetic Beta Diversity

For many applications in ecology, it is important to examine the phylogenetic relations between two communities of species. More formally, let T be a phylogenetic tree and let A and B be two samples of its tips, representing the examined communities. We want to compute a value that expresses the phylogenetic diversity between A and B in T. There exist several measures that can do this; these are the so-called phylogenetic beta diversity (β-diversity) measures. Two popular measures of this kind are the Community Distance (CD) and the Common Branch Length (CBL). In most applications, it is not sufficient to compute the value of a beta diversity measure for two communities A and B; we also want to know if this value is relatively large or small compared to all possible pairs of communities in T that have the same size. To decide this, the ideal approach is to compute a standardised index that involves the mean and the standard deviation of this measure among all pairs of species samples that have the same number of elements as A and B. However, no method exists for computing exactly and efficiently this index for CD and CBL. We present analytical expressions for computing the expectation and the standard deviation of CD and CBL. Based on these expressions, we describe efficient algorithms for computing the standardised indices of the two measures. Using standard algorithmic analysis, we provide guarantees on the theoretical efficiency of our algorithms. We implemented our algorithms and measured their efficiency in practice. Our implementations compute the standardised indices of CD and CBL in less than twenty seconds for a hundred pairs of samples on trees with 7 ⋅ 104 tips. Our implementations are available through the R package PhyloMeasures.


Introduction
Ecologists often distinguish three kinds of diversity. Alpha diversity describes the diversity of one sample (such as the number of plant species in a vegetation plot), beta diversity describes the dissimilarity between a pair of samples, and gamma diversity describes the diversity of a large set of samples [1,2]. These concepts can be applied to a number of diversity measures, including species richness, functional diversity and phylogenetic diversity [3,4]. Phylogenetic beta diversity describes the phylogenetic distance among pairs of communities [5]. It is an increasingly widely-T , we show how to derive analytical expressions for computing exactly the mean and the variance of the CD and the CBL among pairs of species samples in T that have given sizes. Based on these expressions, for each of these two measures we describe efficient algorithms for computing the measure's standardized index for a given pair of samples (and also the measure's original value for this pair). These algorithms are developed based on standard techniques in Algorithms Design, and we provide theoretical guarantees for their performance. We implemented all of the algorithms the we describe, and we conducted experiments that exhibit their efficiency in practice. We measured the performance of our algorithms on trees of several sizes, extracted from a phylogeny that has 71,181 tips and 83,751 nodes in total. As we show later in detail, our implementations ran very fast even for the largest trees that we considered; for each of the two measures, our programs managed to compute in less than a minute the standardized indices for one hundred pairs of species samples on the complete tree with the 71,181 tips. We have made the programs that we developed publicly available through the open source R package PhyloMeasures [24]. To exhibit the strength of fast beta diversity computations, we also present an application example; we compute heat maps that illustrate the beta diversity values between a focal point assemblage and assemblages distributed on a world grid. Given our efficient algorithms, it is now possible to produce several such maps in high resolution, something that was infeasible with previously existing software.

Terminology and Notation
Let T be a phylogenetic tree. We use V to indicate the set of nodes (representing the species/ taxa) in T , and we use E to denote the set of edges (links between nodes) in this tree. For an edge e 2 E, we denote the (always positive) weight of this edge by w e . Depending on the context, the edge weights in T may represent time intervals, molecular distance, or some other notion of difference between taxa. The analysis in this paper does not depend on this notion of difference. We denote the set of leaf nodes in T (taxa that do not have any child species) by S. We refer to these nodes as the tips of T . We indicate the number of the tips in T by sðT Þ or simply s, and we indicate the number of all the nodes in T by n. We consider that T is a rooted tree.
For any edge e be in T , we use Ch(e) to denote the edges that are adjacent to the child node of e. We call these edges the children of e. Let u be a node in T , and let e be the edge that connects u with its parent node. We use interchangeably T ðuÞ and T ðeÞ to denote the subtree of T whose root is u. We use S(u) and S(e) to indicate the set of tips that appear in T ðuÞ. We denote the number of these tips by s(u) and s(e).
Let u, v be two nodes in T . We call a simple path between these nodes the cycle-free sequence of edges that we have to traverse in order to reach u from v. We call the cost of this path the sum of the weights of all the edges in the path. We denote this cost by cost(u, v). Since T is a tree, there exists a unique simple path between any pair of nodes in T . We call the height of T the maximum number of edges that appear on a simple path between the root of T and any leaf. We represent the height of T by hðT Þ. Let R S be any sample (subset) of the tips in T . We denote the number of tips in this sample by |R|. We indicate the set of all paths that connect two elements in R by Paths(R), that is: We denote the set whose elements are all subsets of S that have cardinality exactly r by Sub(S, r). For an edge e 2 E and a subset R of the tips of T , we denote the elements of S(e) that are also elements of R by S R (e), that is S R (e) = S(e)\R. We indicate the the number of these tips as s R (e).
For a given phylogenetic tree T we call the total path cost of T the sum of the costs of all distinct simple paths that connect tips of T . We denote this quantity by TCðT Þ. Thus, the toal path cost of a tree T is equal to: costðu; vÞ: Let e be an edge of T . We call the total path cost of e the sum of the costs of all simple paths in Paths(S) that contain e. We denote this quantity by TC(e), thus: For a node u that is a tip of T , we call the total path cost of u the sum of the costs of all simple paths between u and any other tip of T . We indicate this quantity by TC(u). Note that TC(u) = TC(e), where e is the edge adjacent to u.

The Community Distance
Let T be a phylogenetic tree, and let A, B S be two samples of its tips with |A| = a, and |B| = b. The Community Distance (CD) between A and B is equal to the sum of the costs of all paths that connect a tip in A with a tip in B, divided by the total number of these paths. Therefore, the community distance between A and B is equal to: costðu; vÞ: Samples A and B may not necessarily be of equal size. The CD is analogous to the Mean Pairwise Distance measure for computing the average distance between two samples of species.
The β Net Relatedness Index. Given a phylogenetic tree T and two samples of its tips A, B such that |A| = a and |B| = b, the standardized index of the CD for these samples is equal to: where E CD ðT ; a; bÞ and sd CD ðT ; a; bÞ are respectively the expectation and the standard deviation of the CD for all pairs of tip samples such that one sample contains a tips and the other sample contains b tips. The following theorem provides an analytical expression for the expectation of the CD. Theorem 1. Let T be a phylogenetic tree that has s tips, and let a, b be positive integers with a, b s. The expected value of the CD among all pairs of tip samples in T , such that one sample consists of a tips, and the other consists of b tips, is equal to: which yields the analytical expression for the expectation of the CD measure.
Next we present an analytical expression for calculating the standard deviation of the CD measure.
Theorem 2. Let T be a phylogenetic tree that has s tips, and let a, b be positive integers with a, b s. The standard deviation of the CD among all pairs of tip samples in T , such that one sample consists of a tips, and the other consists of b tips, is equal to: w e Á TCðeÞ À E 2 CD ðT ; a; bÞ r ; where: k 1 ¼ 4ðaÀ1ÞðbÀ1Þ abs 2 ðsÀ1Þ 2 , k 2 ¼ 2ðaÀ1ÞðbÀ1Þ abs 2 ðsÀ1Þ 2 þ bÀ1 abs 2 ðsÀ1Þ þ aÀ1 abs 2 ðsÀ1Þ , k 3 ¼ 2ðaÀ1ÞðbÀ1Þ abs 2 ðsÀ1Þ 2 þ 2 abs 2 . Proof. The standard deviation of the CD is equal to: ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi E A 2 SubðS; aÞ B 2 SubðS; bÞ ½CD 2 ðT ; A; BÞ À E 2 CD ðT ; a; bÞ r : We already provided an analytical expression for the expectation of the CD; hence, we now focus on deriving an expression for the expectation of the squared value of this measure. We get that: costðu; vÞ Á costðx; zÞ Á Pðu; v; A; BÞ Á Pðx; z; A; BÞ if jfu; vg \ fx; zgj ¼ 1: 2aðaÀ1ÞbðbÀ1Þ Thus, we can rewrite Eq (2) as follows: where k 1 ¼ 4aðaÀ1ÞbðbÀ1Þ þ aðaÀ1Þb s 2 ðsÀ1Þ þ abðbÀ1Þ s 2 ðsÀ1Þ , and k 3 ¼ 2aðaÀ1ÞbðbÀ1Þ Next, we simplify the expression that appears in Eq (3). The last of the three sums in this expression can be rewritten as: X fu;vg2S The second double sum in Eq (3) can be simplified as follows: The first double sum in Eq (3) can be rewritten as: costðu; vÞ½TCðT Þ À TCðuÞ À TCðvÞ þ costðu; vÞ ¼ Combining Eqs (5), (6) and (7) with Eq (3) we get the expression that appears in the definition of the theorem.

The Common Branch Length
Let T be a phylogenetic tree and let A, B S be two samples of its tips. The Common Branch Length (CBL) between A and B is the sum of the weights of all edges that appear both in T ðAÞ and in T ðBÞ. Recall that T ðRÞ denotes the smallest subtree in T that contains all the tips of a sample R. Therefore, the CBL between A and B is equal to: Samples A and B may not have the same number of elements. The CBL is analogous to the Phylogenetic Diversity measure for computing a diversity value between two samples of species.
The Common Length Index. For a phylogenetic tree T and two samples of its tips A, B such that |A| = a, and |B| = b, we denote the standardized index of the CBL by: where E CBL ðT ; a; bÞ and sd CBL ðT ; a; bÞ are respectively the expected value and the standard deviation of the CBL for all possible pairs of tip samples where one sample contains exactly a tips and the other contains exactly b tips. We call this index the Common Length Index of A and B. The following two theorems provide analytical expressions for the expectation and standard deviation of the CBL. Theorem 3. Let T be a phylogenetic tree that contains s tips, and let a, b be two positive integers such that a, b s. The expected value of the CBL among all pairs of tip samples in T such that one sample consists of a tips and the other sample consists of b tips, is equal to: Proof. The expectation of the CBL for two independtly selected samples A and B is equal to: and the expression for the expectation follows. Theorem 4. Let T be a phylogenetic tree that contains s tips, and let a, b be two positive integers such that a, b s. The standard deviation of the CBL among all pairs of tip samples in T such that one sample consists of a tips and the other sample consists of b tips, is equal to: In Theorem 3 we provided an expression for the expectation of the CD. It remains to derive an expression for the expectation of the squared value of this measure. We get that: Events fe; l 2 T ðAÞg and fe; l 2 T ðBÞg are independent, so the probability value in the last sum can be rewritten as: Pr½fe; l 2 T ðAÞg \ fe; l 2 T ðBÞg ¼ Pr½e; l 2 T ðAÞ Á Pr½e; l 2 T ðBÞ: Value Pr½e; l 2 T ðAÞ is the probability that both edges e, l appear in T ðAÞ. Recall that T ðAÞ is the smallest subtree of T that contains all tips in A. An edge e does not appear in T ðAÞ if either all, or none of the elements in A are tips in the subtree of e. Given that, we get: For the probability value Pr½ðe = 2 T ðAÞÞ \ ðl = 2 T ðAÞÞ we distinguish three cases: (I). Edge l 2 T ðeÞ: (II). Edge e 2 T(l). This case is symmetric to case (I): (III). Subtrees T ðeÞ and T ðlÞ do not contain each other: The probability value Pr½ðe; l 2 T ðBÞÞ can be expressed in a similar manner. The analytical expression for the standard deviation of the CBL follows by combining Eqs (11)- (17).

Design of Algorithms
Based on the analytical expressions that we presented in the previous sections, we designed efficient algorithms for calculating the value and the standardized indices of the CD and the CBL. Next, we present in short how we designed these algorithms, and we also give a theoretical measure for their efficiency. Before continuing with the description of the algorithms, we explain standard concepts and notation from the field of Algorithms Design that we use to describe the efficiency of our algorithms.
In many applications, it is important to define a simple bound that describes the order of growth for a given function. For example, let G be a set of n real numbers. Suppose that we want to count the number of possible subsets that can be created by picking four elements of G. The exact number of these subsets is equal to f ðnÞ ¼ 1 24 ðn 4 À 6n 3 þ 11n 2 À 6nÞ. But, this expression is quite complicated. Suppose that we want to express how fast f(n) grows as the value of n increases. In that case, for large values of n, the term that influences the value of f(n) the most is n 4 . The standard way to express this in short is to say that f(n) = O(n 4 ). The notation O(n 4 ) is used to indicate that the value of f(n) grows roughly as fast as n 4 when n becomes large. This is known as the Big-Oh notation. More formally, let g(n) and f(n) be two functions.
We say that f(n) = O(g(n)) if there exist two positive constants c and n 0 such that f(n) c Á g(n) for every n ! n 0 .
In the fields of Algorithms Design, the O(Á) notation is used extensively for measuring several aspects of an algorithm's efficiency. The standard way of measuring the efficiency of an algorithm is to count the number of basic operations that would take place during its execution. When we refer to basic operations, we mean all simple mathematical operations (such as comparisons, additions, divisions), but also standard operations like initialising a variable, or assigning a value to this variable. Let A be an algorithm, and suppose we have an input dataset for this algorithm that has size n; for instance a phylogenetic tree T that consists of n elements. Instead of counting exactly the number of basic operations that take place when A processes T , we can bound this number using the Big-Oh notation. Ideally, we would like to prove that, for some function g(n), algorithm A takes O(g(n)) operations to process any input of size n. In that case, we say that the worst case running time complexity of A, or simply the worst case time complexity of this algorithm is O(g(n)).
Next we describe in short the algorithms that we designed for computing the standardized indices of the CD and the CBL. For each of the algorithms that we describe, we also provide a bound for its worst case time complexity. Let T be a phylogenetic tree, and let A and B be two samples of tips in T such that A has a tips and B has b tips. As described in Eqs (1) and (8), to compute the standardized index of either CD or CBL we need to calculate three values; we need to calculate the value of one of these measures for samples A and B, and we need to compute the expectation and the standard deviation of this measure among all pairs of samples in T where one sample has a tips and the other sample has b tips. Next, for each measure, we describe an algorithm for computing each of these three values.

Algorithms for computing the CD and NRI β
Computing the value of the CD for a given pair of samples. To compute the value of the CD efficiently, we rewrite the expression of this measure so that it can be evaluated with a small number of basic operations. Let e be an edge in T , and let Num(A, B, e) denote the number of simple paths that connect a tip in sample A with a tip in sample B, and contain edge e.
Therefore, computing CDðT ; A; BÞ boils down to computing values s A (e) and s B (e) for every edge e 2 T . We can compute all these values in O(n) operations using a simple recursive algorithm; for each edge e that we process, we first calculate values s A (e 0 ) and s B (e 0 ), for every edge e 0 2 Ch(e); recall that Ch(e) is the set of the edges that are adjacent to the child node of e. Then, we can calculate values s A (e) and s B (e) for e by simply adding the corresponding values of the edges in Ch(e). We can do this with only O(n) operations, by traversing the edges in T appropriately. We start from the root of T , and traverse the tree top to bottom. For each edge e that we encounter for the first time, we consider two cases; if e is adjacent to a tip, we check if this tip belongs to A or B, and we set values s A (e) and s B (e) to one or zero accordingly. Then, we move upwards in the tree and we use s A (e) and s B (e) to calculate the corresponding values for the parent edge of e. If e is not adjacent to a tip, then for every edge l 2 Ch(e) we first visit the subtree of l and we calculate recursively the values s A (Á) and s B (Á) for every edge in this subtree. After calculating these values for every l 2 Ch(e), we can compute s A (e) and s B (e) in O(Ch(e)) time by evaluating sums s A (e) = ∑ l 2 Ch(e) s A (l) and s B (e) = ∑ l 2 Ch(e) s B (l). Such a traversal is known as a post-order traversal, and requires visiting each edge in T at most two times. Since T has O(n) edges, this traversal requires O(n) operations. We also perform O(Ch(e)) arithmetic operations to compute s A (e) and s B (e) for every edge in T , which sums up to O(n) time operations in total. After calculating values s A (e) and s B (e), we can compute CDðT ; A; BÞ with O(n) arithmetic operations using Eqs (18) and (19). Therefore, the total number of operations that are needed to compute the value of the CD for two tip samples in T is O(n).
Computing the expectation of the CD. In Theorem 1 we provide an analytical expression for the expectation of the CD. Evaluating this expression boils down to evaluating TCðT Þ. Recall that TCðT Þ is equal to the sum of the costs of all simple paths that connect pairs of tips in T . We can rewrite this values as: TCðT Þ ¼ X e2T sðeÞ Á ðs À sðeÞÞ: Therefore, to compute the expectation of the CD, it remains to compute s(e) for every edge e in T . This can be done in O(n) operations with a post-order traversal of T , as described also for the algorithm that calculates the value of the CD for two samples. Hence, we can compute the expectation for this measure with only O(n) operations.
Computing the standard deviation of the CD. For the standard deviation of the CD we presented an analytical expression in Theorem 2. This expression contains sums that have O (n) terms in total. Among other terms, these sums contain values TCðT Þ and TC(e) for every edge e in T . In the algorithm that we described for computing the expectation of the CD, we showed already how we can compute TCðT Þ with O(n) operations. We can also compute TC (e) for every edge e in the tree, with O(n) operations in total. This can be done as follows. For an edge e, we use w T ðeÞ to denote the value: For any edge e in the tree, value TC(e) is equal to: TCðlÞ þ ðsðeÞ À sðlÞÞ Á w T ðlÞ þ sðlÞ Á ðw T ðeÞ À w T ðlÞÞ: ð20Þ To compute the standard deviation of the CD, we perform two post-order traversals of the tree T . In the first traversal, we compute value w T ðeÞ for each edge e based on the corresponding values of the edges in Ch(e). In the second traversal, we use the values w T ðÁÞ that we just calculated to compute TC(e) for every edge in the tree. Then, using these values we can evaluate the expression in Eq (20). It takes O(n) operations to traverse the tree twice and, given values TCðT Þ and T ðeÞ for every edge in T , it takes O(n) arithmetic operations to calculate the expression in Eq (20). Hence, we can compute the standard deviation of the CD with O(n) operations in total.
Using the algorithms that we describe above, we can compute the value of the CD for two samples A and B, but also the expectation and the deviation of this measure in O(n) time in the worst case. Therefore, by combining these algorithms, we can derive an algorithm that computes NRI β in O(n) time in the worst case.

Algorithms for computing the CBL and CLI
Computing the value of CBL for a given pair of samples. Recall that, for two tip samples A and B, the value of the CBL is equal to the sum of the weights of all edges e such that e belongs to both subtrees T ðAÞ and T ðBÞ. Let R be a sample that consists of |R| = r tips in T . If edge e belongs to subtree T ðRÞ, then it holds that 0 < s R (e)<r. Hence, the weight of an edge e is counted in the value CBLðT ; A; BÞ if 0 < s A (e)<a and 0 < s B (e)<b. Therefore, to decide for every edge e 2 T whether to include w e in the CBL value, we simply have to compute values s A (e) and s B (e). In the algorithm that computes the value of CD, we showed how we can compute values s A (e) and s B (e) with O(n) operations. Thus, the total number of required operations for computing the value of CBL is also O(n).
Computing the expectation of the CBL. Theorem 3 provides an analytical expression for the expectation of the CBL. This expression consists of a sum with O(n) terms, and can be evaluated with O(n) arithmetic operations, given that we have already computed value s(e) for each edge e in the tree. When describing an algorithm that calculates the expectation of the CD, we showed how we can compute these values with O(n) operations. Therefore, we can compute the expectation of the CBL in O(n) time.
Computing the standard deviation of the CBL. The standard deviation of the CBL can be calculated using the expression that we provide in Theorem 4. This expression contains a sum of O(n 2 ) terms. Evaluating this sum directly takes O(n 2 ) operations, which can be inefficient in practice.
In a previous paper, we presented a technique for computing similar expressions [13]. We can use this technique for evaluating the sum in Eq (9). As we presented in our previous work, this technique is quite efficient when T is relatively balanced. More precisely, the worst case time complexity of this algorithm is OðSIðT ÞÞ, where SIðT Þ denotes the Sackin's Index of a tree T . The Sackin's Index of T is equal to the sum of the depths of all tips in T . The depth of a tip v is equal to the number of edges that appear on the simple path between v and the root of the tree. For a perfectly balanced tree that consists of n nodes, the depth of each tip is O(log (n)). In this case, the value of of the Sackin's Index is equal to O(n log n). However, if the tree is skewed, there may exist many tips with depth close to n and therefore the value of the Sackin's Index is O(n 2 ). Hence, in the worst case, the algorithm that we consider for computing the standard deviation of the CBL takes quadratic time with respect to the size of the input tree. In practice, phylogenetic trees are relatively balanced, and therefore the proposed algorithm is quite efficient. In the next section, we provide evidence for this argument; there we present experiments that we conducted using our algorithms on large tress.

Experimental Evaluation
We implemented all the algorithms that we present in this paper, and we measured their performance in practice. Our implementations are developed in C++, and are publicly available through the software package PhyloMeasures [24]. This package provides functions for computing the value and the standardized indices of several phylogenetic biodiversity measures, and it is available both as an R package and a C++ library.
For the experiments that we conducted, we used a large phylogenetic tree from which we extracted several subtrees of several sizes. The tree that we used was constructed by Goloboff et. al [25]. This is the largest evolutionary tree of eukaryotic organisms that has been so far constructed from molecular and morphological data. It consists of 71,181 tips and 83,751 nodes in total. This tree is unrooted; for the needs of our experiments we picked arbitrarily an internal node and used this as the root. We call this dataset the eukaryotes dataset.
From the eukaryotes we extracted fourteen trees, each tree having 5,000k + 1181 tips with k 2 {1, 2, . . ., 14}. These subtrees were produced by successively pruning chunks of 5,000 leaves from the eukaryotes tree. We represent the set of these trees by EK. Therefore, for any two trees T ; T 0 2 EK we have that either T is a subtree of T 0 , or vice versa. For each tree T 2 EK we produced one hundred samples of tips with sizes sðT Þ=k with k ranging from one to a hundred. We denote this set of samples by samples ðT Þ. For every T 2 EK we executed the algorithms we implented for computing the values and the standardized indices of the CD and the CBL. For each algorithm and for each tree T 2 EK, we measured the execution time for processing all samples in samples ðT Þ. The results of these experiments are illustrated in Fig 1. The experiments were performed using the R version of the PhyloMeasures package on a computer with an Intel core i5-2430M processor. This is a four-core CPU with 2.40GHz per core. The main memory of this computer is 7.8 Gigabytes. Our implementations run on a Linux Ubuntu operating system, release 12.04. The experiments were executed using R version 3.1.2 (Pumpkin Helmet).
All of the examined implementations run very fast even for the largest trees in EK. For the complete eukaryotes tree that consists of 71,181 tips, the algorithm that computes the value of the CD takes 3.01 seconds to process one hundred samples, the algorithm that computes the CBL value takes 1.61 seconds, and the algorithms that compute the standardized indices of the CD and the CBL take 3.12 and 54.52 seconds respectively. Note that these are the running times that each program takes for computing the results for one hundred samples of species. The program that computes the index of the CBL runs slower than the other programs, yet we see that it is still quite efficient. It seems that, for the datasets that we used, its execution time does not scale as a quadratic function of the tree size.

Applications
To illustrate an application of the algorithms, we used their PhyloMeasures implementation to calculate CBL and standardized CBL for mammals globally. We used range maps from Efficient Algorithms for the Phylogenetic Common Branch Length and Community Distance the IUCN and rasterized them on a Behrmann equal area grid with a resolution of 193km (at 30N or S), corresponding roughly to two degrees. Each of the resulting grids consists of 71 × 180 cells, which means 12780 cells in total for each grid.
We combined these grids with the phylogeny of Bininda-Emonds et al [26]; this is a tree that portrays the phylogenetic relations between all mammal species. The number of leaf nodes in this tree is 4,510. Using this tree, for each grid that we created we then calculated CBL and standardized CBL between a focal cell and all other cells in the world. The results are spatial maps of the phylogenetic similarity between the focal location and all other sites. The constructed spatial maps are illustrated in Fig 2. In Central Asia, for example, there is evidence for a broad longitudinal band of high similarity, with more rapid turnover along the latitudinal gradient, while phylogenetic similarity is fairly high throughout northern South America, but declines rapidly towards the south of the continent.
Such maps can provide a detailed picture of how phylogenetic similarity changes among species assemblages over geographic space. However, constructing high resolution maps of this kind is practically infeasible without efficient algorithms that compute beta-diversity values, and especially the standardised version of these values. Using our implementations, we computed the exact standardized and non-stanndardized CBL values for these maps in one and a half minutes in total. On the other hand, it would take several days to execute these computations with previously existing software, which would also provide these values only approximately.

Conclusion
Phylogenetic beta diversity metrics are widely and increasingly used in ecology, but are often quite slow to compute. In addition, the relationship between species richness of the samples and the beta diversity metrics is often ignored. These problems are related, as computing the dependence of a metric on species richness is even more computationally intensive using traditional approaches, than is computing a single beta diversity value. Here, we propose solutions to these problems, providing 1) efficient algorithms for computing phylogenetic beta diversity metrics, and 2) algorithms to efficiently and exactly calculate their moments, allowing a simple standardization for species richness. We expect that these algorithms will significantly ease the computational burdens of researchers and lead to wider adoption of phylogenetic beta diversity metrics.