Emmanuel Paradis
email: emmanuel.paradis@ird.fr

The distribution of branch lengths in phylogenetic trees

Keywords: diversity, extinction, phylogeny, speciation

A lot of effort has been devoted to analyse the distribution of branching times observed in a phylogenetic tree. On the other hand, the distribution of branch lengths has not received similar attention. In this paper, the distribution of branch lengths is studied. It is shown that different types of branches within a tree have distinct distributions. Some equations to predict these distributions are derived with respect to diversification parameters and whether the size of the tree is known or not. A simulation study validated these predictions. The inferred distributions are used to develop graphical and statistical tools to assess the goodness-of-fit of diversification models. An application is presented on a recently published dated phylogeny of Carnivora. Some future developments are discussed.

Introduction

Phylogenetic trees have become a foundamental tool to study macroevolutionary processes.

Several information can be extracted from a phylogeny: its general tree shape or various summaries of branch lengths [START_REF] Mooers | Inferring evolutionary process from phylogenetic tree shape[END_REF][START_REF] Pybus | Testing macro-evolutionary models using incomplete molecular phylogenies[END_REF][START_REF] Franc ¸ois | Gaussian approximations for phylogenetic branch length statistics under stochastic models of biodiversity[END_REF]. These information can then be used to infer the speciation and extinction rates of the studied group thus quantifying the tempo of its diversification through time. An important class of models considers the homogeneous models which assume that the rates of speciation and extinction may vary with respect to time but at a given time they are the same for all species present in the clade. From a biological point of view, these models are justified when diversification is driven by global environmental variables such as climate. An important result from [START_REF] Lambert | Birth-death models and coalescent point processes: the shape and probability of reconstructed phylogenies[END_REF] is that under a homogeneous model of diversification, the distribution of branching times depends only on the rates of speciation and extinction. Therefore, the branching times of a phylogenetic tree summarize the information required to estimate these rates under the assumption of homogeneity. However, when this assumption does not hold, one has to consider other quantities such as the distribution of branch lengths. For instance, [START_REF] Venditti | Phylogenies reveal new interpretation of speciation and the Red Queen[END_REF] analysed the distribution of internal branch lengths for 101 phylogenies in order to assess the constancy of speciation rates accross a wide range of taxonomic groups. Besides, the distribution of branch lengths had retained some attention, for instance for the study of phylogenertic diversity [START_REF] Faith | Conservation evaluation and phylogenetic diversity[END_REF][START_REF] Mooers | Branch lengths on birth-death trees and the expected loss of phylogenetic diversity[END_REF][START_REF] Stadler | Distribution of branch lengths and phylogenetic diversity under homogeneous speciation models[END_REF]. These previous studies have mostly considered models of constant diversification, i.e., with constant rates of speciation and extinction.

The main objective of this paper is to develop an approach to infer the distribution of the branch lengths of a phylogenetic tree considered as the result of a diversification process with known rates of speciation and extinction. Contrary to previous contributions, the present approach handles time-dependent variation in speciation and extinction rates. Two cases are considered: when the size of the tree is unknown, and when it is known so that the predictions of the mean branch length can be conditioned on tree size. A simulation study was performed to assess the accuracy of the predictions of the equations derived in this work. Some applications of these equations are considered with a focus on a method to assess the goodness of the fit of a diversification model using graphical tools and statistical tests based on the cumulative distribution function. An analysis of a phylogeny of Carnivora (Mammalia) is presented to illustrate their practical use.

Models

THE GENERAL CASE UNCONDITIONAL ON TREE SIZE

We consider a dated phylogenetic tree with N species. We assume that this tree is rooted and fully dichotomous; consequently, it has N -1 (= m) internal nodes. Several quantities can be extracted and analysed from such a tree (Fig. 1). The branching times (also known as node times, node heights, or node depths), measured from the root, are denoted t 1 , . . . ,t m . In this section, the root is the origin of the time scale, so that t 1 = 0, and the present is denoted T . Alternatively, we could consider the branching intervals (known as coalescent intervals in population genetics) which summarise the same information as the branching times and are defined as t i+1t i (i = 1, . . . , m -1). The branch lengths provide more information than the branching times because different sets of branch lengths can result in the same branching times. However, one needs the topology of the tree to fully determine the branching times from the branch lengths.

Birth-death models are a class of continuous-time models where each species (or individual) is continuously exposed to speciation (birth) and extinction (death). The time-dependent birth-death model assumes that the probability of a speciation event varies with time and follows a function λ(t), and similarly for the probability of extinction µ(t). This formulation implies no lineage-specific variation: at a given time, all species have the same probabilities of speciation and extinction. We first remind some general mathematical background.

Under a time-dependent birth-death model, we can write the probability that a lineage, originating from a single species at time t, has exactly one descendant at time T [START_REF] Kendall | On the generalized "birth-and-death" process[END_REF]:

Pr(t, n T = 1) = e -ρ(t,T ) W (t) 2 , with ρ(t, T ) = T t µ(u) -λ(u)du, and
W (t) = e -ρ(t,T ) 1 + T t e ρ(t,u) µ(u)du .
The probability of a lineage being extinct during the same time interval is:

Pr(t, n T = 0) = 1 - e -ρ(t,T ) W (t) ,
from which we can write the probability that a lineage is not extinct:

Pr(t, n T ≥ 1) = 1 -Pr(t, n T = 0) = e -ρ(t,T ) W (t) .
We easily deduce the probability that a lineage has two or more lineages which we will use below:

Pr(t, n T ≥ 2) = 1 -Pr(t, n T = 0) -Pr(t, n T = 1).
The number of species living at time t is a random variable with expectation given by: 0,t) .

E(n t ) = e -ρ(
These equations are the basis to derive the distribution of branching times of an ultrametric phylogenetic tree [START_REF] Nee | The reconstructed evolutionary process[END_REF] [START_REF] Mckenzie | Distributions of cherries for two models of trees[END_REF], and those that are sister-group of a clade of two or more species. These two types of terminal branches are termed cherry and outer branch, respectively (Fig. 2). The cherry branches are by definition duplicated in a tree.

A cherry is the result of a speciation at time t leading to two species surviving at present (time T ) with no other (recorded) speciation event. Under the assumption that these events do not covary, the density ξ 1 (t) of a cherry is given by the product of the probabilities of these three events (speciation at t and survival of both species from t to T ):

ξ 1 (t) = λ(t)Pr(t, n T = 1) 2 . (1)
Note that this is a density function, not a probability function because we cannot assume that ξ 1 (t) sums to one over all possible values of t.

We can use a similar reasoning and apply it to derive the density of outer branches. We consider here three events: a speciation at time t, the survival of a single species from t to T , and the survival of a clade with at least two species from t to T . There is an additional factor 2 because there are two possible combinations (i.e., left and right children from the speciation at t).

Thus, the density ξ 2 (t) of a single terminal branch originating at time t and being sister of a clade with two or more species is:

ξ 2 (t) = 2λ(t) Pr(t, n T = 1) Pr(t, n T ≥ 2). (2)
These two density functions ξ 1 (t) and ξ 2 (t), after being multiplied by E(n t ) and properly normalized, lead to probability density functions (pdf) of the distribution of both types of terminal branch length. The respective expected means can be calculated with:

E(l) = T 0 u ξ(u)du, ( 3
)
where ξ is the normalized function and l is the branch length. The normalizing factor is computed with:

T 0 E(n u )ξ(u)du,
where ξ is replaced by either ξ 1 or ξ 2 depending on the type of branch. The cumulative density function (CDF) of branch lengths, which is by definition the probability that l is smaller than or equal to a given value (t) is:

Pr(l ≤ t) = t 0 ξ(u)du,
The expected mean terminal branch length over the tree (i.e., pooling both types) is computed by the mean of both means weighted by their normalizing factors.

We now turn to the distribution of internal branch lengths. The problem is more difficult as it requires to derive the probabilities of branches connecting two internal nodes. The first (oldest) node connects two sister-clades, one of them includes two or more species and the other one includes one or more species. The density function of such a node being observed at time t is denoted as g 1 . This node is the result of three events: a speciation at time t, the survival of a clade with two or more species from t until T , and the survival of a clade with one or more species during the same time. Thus:

g 1 (t) = 2λ(t) Pr(t, n T ≥ 2) Pr(t, n T ≥ 1).
The factor 2 is, again, because there are two possible combinations for the sister-clades. The second (youngest) node connects two clades both made of one or more species, and its density function is:

g 2 (t) = λ(t)Pr(t, n T ≥ 1) 2 .
The density of internal branches is found by multiplying the density functions of these two nodes and the probability of neither extinction, nor recorded speciation between t and t . The probability of this last event is precisely given by the function g 2 . It is thus necessary to integrate between t and t the sum of the two functions µ(t) and g 2 (t) [START_REF] Cox | Analysis of Survival Data[END_REF]. We finally find that the density of an internal branch starting at time t and ending at time t is:

f (t,t ) = g 1 (t)g 2 (t ) exp - t t g 2 (u) + µ(u)du .
As before, we can use this density, after proper normalization, as a pdf of internal branch lengths. In this case, calculating the expected mean requires a double integration on t and t (with the constraint 0 ≤ t < t < T ):

E(l) = T 0 T t (t -t) f (t,t )dt dt. ( 4
)
Note that tt is the internal branch length (l), and f is, as above, the normalized density.

CONDITIONING ON TREE SIZE

The above development assumes that T (the age of the root of the tree) is known whereas N is unknown. In this section, it is assumed that N is known and fixed whereas T is a random variable.

To derive the distributions of branch lengths, a different approach is needed. The approach used here is inspired from the coalescent approach developed by [START_REF] Stadler | Lineages-through-time plots of neutral models for speciation[END_REF][START_REF] Stadler | On incomplete sampling under birth-death models and connections to the sampling-based coalescent[END_REF]Stadler ( , 2011a)). The coalescent [START_REF] Kingman | The coalescent[END_REF] considers the way a sample of gene lineages from a population coalesce backwards in time. In population genetics, the distribution of the times to coalescence is determined by the genealogical structure within the population [START_REF] Wakeley | Coalescent theory: an introduction[END_REF]. On the other hand, in a phylogenetic tree, the coalescence events follow the speciation and extinction events (see Stadler, 2011a, for a clear exposition of this rationale). This approach leads in a straightforward way to the distribution of branching times in a phylogeny conditioned on N (see below); however, we also need to know the distribution of cherry and outer terminal branches.

Consider a coalescent tree of size N; each node of this tree is the result of a coalescent event between two lineages. These lineages may be either a lineage not yet coalesced, or the result of a previous coalescent event (i.e., lineages already clustered). For simplicity, we may call the first kind of lineage 'singleton' and the second kind 'cluster'. Let α i and ω i denote the respective number of these two kinds of lineages, where the subscript i denote the coalescent event (i = 1, . . . , N -1) with i = 0 designing the initial state (i.e., present time). Clearly, we have α 0 = N and ω 0 = 0 (no coalescent event have yet occurred). The first event (i = 1) is the coalescence of two singleton lineages so that α 1 = N -2 and ω 1 = 1. The second coalescence event may join, either two singleton lineages, so α 2 = N -4 and ω 2 = 2, or a singleton lineage with the cluster lineage built at the previous step, so α 2 = N -3 and ω 2 = 1. We can continue along the successice steps of the coalescent process, and find that, under the assumption that lineages coalesce randomly, α i and ω i follow a random bivariate discrete process with the three following possible transitions with their respective probabilities on the right-hand side:

α i+1 = α i -2 ω i+1 = ω i + 1      α i (α i -1) 2K i , (5)
α i+1 = α i -1 ω i+1 = ω i      α i ω i K i , (6)
α i+1 = α i ω i+1 = ω i -1      ω i (ω i -1) 2K i , (7)
with

K i = (α i + ω i )(α i + ω i -1)/2.
The first transition ( 5) is the coalescence of two singleton lineages; the second transition ( 6) is the coalescence of one singleton lineage with one cluster lineage; and the third transition ( 7) is the coalescence of two cluster lineages. K i denotes the number of combinations among the α i + ω i lineages not yet coalesced. Because α and ω are discrete values, it is straightforward to compute their pdf for each coalecent event (see code provided with this paper). Trivially, the first coalescent event (i = 1) is of the first type (since ω 0 = 0), and at the end of the coalescent process (i = N -1), we have α = 0 and ω = 1 (Fig. 3).

We note that at each step of a coalescent the number of lineages not yet coalesced is given by Ni and is equal to α i + ω i , so the distributions of these two variables are symmetric.

From this random process it is possible to derive the probability that a branching time leads to a cherry or to an outer branch: the first possible transition (5) leads to the creation of two cherry branches, while the second one (6) leads to the creation of an outer branch. The probability that the coalescent event i results in a cherry is thus given by eq. 5 summed over all possible values of α i weighted by their probabilities as derived above:

N ∑ x=0 Pr(α i = x) α i (α i -1) 2K i .
Similarly, the probability the same coalescent event results in an outer branch is given by eq. 6

summed over all possible values of α i and ω i :

N ∑ α i =0 N ∑ ω i =0 α i ω i K i .
However, we do not have yet the distribution of the lengths of these two kinds of branches.

For this we need the distribution of branching times. In population genetics, this distribution is related to the mutation rate and the demographic structure and dynamics of the population [START_REF] Wakeley | Coalescent theory: an introduction[END_REF], but this cannot be applied to phylogenetic trees. Using the description of the coalescent process generating a phylogenetic tree in Stadler (2011a), it is possible to obtain this distribution in a straightforward way. To use this approach, we have to assume that the birth-death process is time-reversible which is a reasonable assumption if the probabilities of speciation and extinction are not affected by the past states of species (in other words, if the diversification process is Markovian). The successive events considered here are speciations and extinctions. When looking at the tree backwards in time (i.e., from present to the past), a speciation leads to the coalescence of two lineages, whereas an extinction leads to the appearance of a lineage (which does not survive until present). Thus, at time t during this process, the clade is constituted of some lineages surviving until present and some lineages that go extinct before present and which are not recorded in the observed phylogenetic tree. Let S t and E t denote the number of each type of lineages at time t. Note that, by contrast with the previous section, time is now measured from present (t = 0). Under the assumption of independence of events (speciation, extinction and coalescence), the changes through time of these two variables can be modelled by the following pair of differential equations:

dS t dt = -λ(t)n t S t n t 2 , dE t dt = n t µ(t) -λ(t) 2S t E t n 2 t -λ(t) E t n t 2 ,
with n t = S t + E t . The initial state is given by S 0 = N and E 0 = 0. These equations are solved until S t = 1. This can be done numerically after setting the value of N and the functions λ(t) and µ(t), for instance with the package deSolve [START_REF] Soetaert | Solving differential equations in R: package deSolve[END_REF]. The solutions give the expected values of S t and E t through time (Fig. 4). Thus, the value t for which S t = Ni is the expected value of the ith coalescent time of the phylogenetic tree. We then just have to weight these expected values with the probabilities obtained above to derive the expected mean cherry and outer branch lengths.

Like we have seen in the previous section, the calculations for internal branches are more complicated because we have to consider the times of two nodes in the tree. Let A j i be the probability that a node created during the coalescent at step j is available for coalescence at step i with j < i. Clearly, a node is created at each step of the coalescent so this node is available for coalescence at the next step, thus we have A i i+1 = 1. For a node created before step i, the probability that it does not coalesce at step i is given by the ratio of the number of combinations excluding this node on the total number of combinations:

N -i -1 2 N -i 2 = (N -i -1)(N -i -2) 2 × 2 (N -i)(N -i -1) = N -i -2 N -i .
We may thus write how A j i changes along the coalescent process:

A j i+1 =      1 if j = i, A j i N -i -2 N -i otherwise.
We now use A j i to derive the frequency distribution of node ages at all steps of the coalescent, and thus compute the expected length of an internal branch. The reasoning is the same than for terminal branches: a recursive calculation is done for all coalescent events, computing the probabilities of different types of events. The difference here is that instead of considering the branching times, we consider their difference weighted by the probabilities that an internal branch is made for each interval.

Simulation Study

METHOD

To validate the above predictions, some trees were simulated using three algorithms. The first algorithm is described in [START_REF] Paradis | Time-dependent speciation and extinction from phylogenies: a least squares approach[END_REF] and simulates a tree in continuous time where λ and µ are allowed to vary with time under any user-specified model; this is a time-forward algorithm where T is fixed and N is a random variable. It is implemented in the R package ape [START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF] with the function rbdtree. The second algorithm is due to Stadler (2011b) and simulates trees in continuous time using a time-backward method with both N and T fixed; it is implemented in the R package TreeSim [START_REF] Stadler | TreeSim: Simulating trees under the birth-death model. R package version 2.1[END_REF] with the function sim.bd.taxa.age.

The third algorithm is similar to the second one except that T is random: the function rphylo in ape was used. All algorithms require to specify the values of the parameters λ and µ. Different values of these parameters were used (see Table 1). For each set of values of λ and µ, 10,000 trees were simulated using each algorithm. For the second and third algorithms, N was fixed to 200 (this value did not influence significantly the results reported here).

The agreement between the model predictions and the simulations were assessed by comparing the predicted means with the observed means of each type of branch length.

Additionally for some simulated trees, the distribution of branch lengths was plotted and compared with the distribution inferred from eqs. 1 and 2.

RESULTS

The distribution of terminal branch lengths differed markedly between the two types of branch:

the cherry branches had an exponential-like distribution whereas the outer branches had a unimodal distribution (Figs. 5 and6). The exact shape of these distributions depended on the values of speciation and extinction rates. The histograms displayed on Figs. 5 and6 show the observed distributions for two trees simulated with the same net diversification rate (λµ = 0.1)

but with different values of the rates λ and µ. The curves superimposed on the histograms are the densities predicted from eqs. 1 and 2.

When the trees were simulated with fixed N and fixed T (package TreeSim), the observed mean branch lengths were well predicted by eqs. 3 and 4 when the value of µ was small relative to the value of λ. However, the difference between the observed and the predicted values increased when the values of µ increased (Table 1). Similarly, the predictions of the mean internal branch length were not very accurate except when µ was small; however, the predicted values varied in the same way as the observed ones.

When the trees were simulated with random N (package ape), the predicted values for terminal branches were close to the observed ones. For both types of terminal branches, the observed means were very close to the predicted values: the difference between both values was less than 0.1 with one exception (the cherry branches with λ = 0.1 and µ = 0.09) where the difference was less than 0.5 (Table 2). The prediction of mean internal branch lengths was not accurate though, as before, the predicted and observed values varied in the same way.

When the trees were simulated with fixed N and random T , the observed mean cherry or outer branch length (Table 3) were slightly different from the predicted means not conditioned on N (Table 1 or 2). However, when the prediction was conditioned on N = 200 (within parentheses in Table 3), the predicted and observed means were in good agreement for all values of λ and µ.

The same observation was made for internal branch lengths.

Application

The ability to predict the distribution of branch lengths may have several practical applications, such as predicting variation in phylogenetic diversity [START_REF] Faith | Conservation evaluation and phylogenetic diversity[END_REF], or deriving likelihood functions for estimation of diversification parameters. Another potential application could be in specifying informative priors for Bayesian phylogenetic analyses: assuming distributions on speciation and extinction rates (for instance, obtained from data on species diversity), it may then be possible to derive prior distributions on the different kinds of branch lengths to use as input for a Bayesian approach on phylogenetic inference.

AN EXAMPLE WITH GOODNESS-OF-FIT TESTS

Here, we detail an application of the above equations to assess the goodness of fit of diversification models. Methods to assess the goodness-of-fit of a model must be distinguished from methods that compare models (e.g., ratios of mean squares, likelihood-ratio tests, information criteria). The latter seek to test the adequacy of two or more models relatively to each others. By contrast, goodness-of-fit methods assess the adequacy of a model in an absolute manner by comparing the predictions of the model with the data. These methods may be graphical (e.g., plot of residuals in a regression analysis) or statistical (e.g., testing the null hypothesis that the model does not fit the data).

The predicted distributions of branch lengths can be used in different ways to assess the adequacy of a specific birth-death model to a phylogeny. It is possible to draw a histogram of the observed distribution of branch lengths and compare it to the expected density. We must be careful that the shape of the histogram can be influenced by the definition of the intervals.

Another method, which avoids this problem, is to compare the observed and expected quantiles of the distribution of branch lengths (a method known as the QQ-plot). Appendix A details a statistical procedure based on this approach.

Computer code written in R (R Core Team, 2014) is provided with this article to perform all calculations described in this paper.

APPLICATION TO A CARNIVORA PHYLOGENY

Nyakatura and Bininda-Emonds ( 2012) published a dated phylogeny with 286 species of Carnivora and 8 other mammal species as outgroup which were removed from the tree before analysis. The "best estimate" tree provided by Nyakatura and Bininda-Emonds was used. The dated phylogeny was analysed with the method from [START_REF] Nee | The reconstructed evolutionary process[END_REF] to estimate λ and µ; the 95% confidence intervals (CIs) were calculated with the method described in [START_REF] Paradis | Analysis of diversification: combining phylogenetic and taxonomic data[END_REF].

From the above estimates, the expected distributions of terminal branch lengths were derived and compared to the observed ones. Giving that the prediction of internal branch lengths is less precise than for terminal branches (see previous section), only the latter were considered.

The estimated speciation and extinction rates from the dated phylogeny were: λ = 0.134 (95% CI: 0.097-0.182) and μ = 0.032 (95% CI: 0.008-0.068). The distribution of terminal branch lengths showed a unimodal shape for both types (Fig. 7). For the cherry branches, there were less short branches than expected as revealed by the histogram and the QQ-plot (Fig. 7A).

For the outer branches, the distribution of branches of length between 0 and 10 Ma seemed to agree with the expectation of the fitted constant-rate model; however, longer branches were more frequent than expected (Fig. 7B).

Discussion

The present paper proposes to infer the distributions of terminal branch lengths by distinguishing two types: the cherry branches which lead to two sister-species, and the outer branches which lead to a single species sister of a clade with at least two species. The simulations confirmed that these two kinds of terminal branches have distinct distributions. Besides, both types of terminal branches are not expected to be observed with the same frequencies. In a phylogenetic tree, cherry branches are, on average, twice more frequent than outer ones; so, the mixture of both is dominated by the exponential nature of the distribution of the first type.

The simulation results were substantially different depending on the algorithm used to simulate the trees. It is interesting to note that the mean length of cherry branches did not vary much with respect to the simulation algorithm used (first column in Tables 123), while a substantial difference was observed for the outer (second column) and the internal branches (third column). It is noteworthy that the predictions matched well the simulation results when the same conditioning (or lack of) was used. The only notable discrepancy relates to the mean internal branch lengths with random N which were not well predicted when the extinction rate was high (Table 2). On the other hand, the predictions conditioned on the value of N were particularly more accurate for all types of branches (Table 3). This shows the advantage of conditioning the predictions when the number of species in the phylogeny is known.

Two approaches were used in this paper. The first one considers time-forward equations where tree size is a random variable. The second approach is based on time-backward equations with known tree size. Both approaches have their respective advantages: the first one makes possible to derive the complete distribution of branch lengths, while the second one gives more accurate predictions of mean branch lengths. It seems possible to use this second approach to derive the complete distribution of branch lengths conditioned on N, but this has not been yet worked out. This is currently under study. Nevertheless, it appeared here that the predicted distributions from the unconditioned case is already accurate for large trees (see Fig. 6). So it will be particularly interesting to assess whether conditioning on N improves these predictions.

The equations presented in this paper may have several applications. The branch lengths of phylogenetic trees are often used to quantify the diversity of a clade because two distantly related species (thus separated by a long path on a phylogeny) likely represent more biological diversity than two closely related species [START_REF] Faith | Conservation evaluation and phylogenetic diversity[END_REF]. This question is still under current research [START_REF] Volkmann | Prioritizing populations for conservation using phylogenetic networks[END_REF]. The probability density function of branch lengths with respect to diversification parameters is thus likely to be helpful to predict the distribution of phylogenetic biodiversity when reconstructed phylogenetic trees are not available. For instance, speciation and extinction rates can be estimated from species diversity of clades (e.g., [START_REF] Magallón | Absolute diversification rates in angiosperm clades[END_REF][START_REF] Paradis | Quantifying variation in speciation and extinction rates with clade data[END_REF][START_REF] Stadler | On age and species richness of higher taxa[END_REF]. These estimates could then be used together with the eqs. 1 and 2 to predict the phylogenetic diversity of such groups even without phylogenies. In a recent paper, [START_REF] Crawford | Diversity, disparity, and evolutionary rate estimation for unresolved yule trees[END_REF] derived expressions for the distributions of phylogenetic diversity and of phenotypic disparity under the assumption of diversification under a Yule model (so µ = 0 and λ constant) so that these quantities may be estimated event without a known phylogeny.

The application with the Carnivora phylogeny had mainly an illustrative purpose. It showed how to use graphical tools to assess the excess or deficit of short and long branches within a tree.

It should be kept in mind that such diagnostics of fit are done with respect to a fitted model. For simplicity, a simple model was used in the present analyses, but it is straightforward to use more complicated models. The results showed that the Carnivora tree has a deficit of short cherry branches compared to what is expected under a null model. It is clear that such deficit may result in a an apparent decline of speciation rates close to present, a pattern that has been frequently observed in real phylogenies [START_REF] Moen | Why does diversification slow down?[END_REF]. The tools introduced here may thus contribute to investigate this widely debated issue.

In a previous work, [START_REF] Venditti | Phylogenies reveal new interpretation of speciation and the Red Queen[END_REF] quantified the distribution of internal branch lengths using standard statistical distributions and found that the exponential distribution was the best fit for the majority of phylogenies. They interpreted this result as evidence that speciations are the results of single rare events leading to reproductive isolation. Interestingly for the present work, they excluded all terminal branches under the motivation that these do not record information on waiting times between speciation events. The developments presented here may bring a new perspective on such analyses by including all information from the phylogenies.

The possibility to derive separate distributions for cherry and outer terminal branches may have some implications for the study of tree shape and its balance (Tarver and Donoghue, 2011).

In particular, it seems possible to derive under which conditions unbalanced trees are generated.

It must be pointed out that the main objective of this paper was to set a framework to derive the distributions of different types of branch lengths within a tree. For simplicity, it was implicitly assumed in the above developments that diversification is homegeneous (λ and µ can vary through time but they have the same values for all species at a given time). However, it is possible to generalise this approach to situations with heterogeneous rates (e.g., when some contemporaneous clades do not diversify at the same rates). The fact that the above equations ( 1and 2) can consider a single branch makes possible to derive distributions which consider heterogeneous rates. A possible application could be, for instance, with the coalescent approach used with fixed tree size, to relate diversification parameters to ancestral character states, and thus analyse a trait-dependent diversification model. This is currently under development.

difficulty with these ECDF-based tests is that their exact distributions are not known, and therefore the critical values must be found by simulations of the null hypothesis [START_REF] Stephens | provide a revisionist perspective of evolutionary history in topological analyses of diversity[END_REF].

Here we consider two tests: the Cramér-von Mises test:

W 2 = n ∑ i=1 z i - 2i -1 2n 2 + 1 12n ,
where n is the number of points considered for the test and z i is the predicted CDF for the ith point. The second test was proposed by [START_REF] Zhang | Powerful goodness-of-fit tests based on the likelihood ratio[END_REF]:

Z C = n ∑ i=1 log 1/z i -1 (n -0.5)/(i -0.75) -1 2 .
For both tests, the null distribution was determined by simulating a large number of trees (typically 10 4 ). The procedure is as follows.

1. Obtain the estimates λ and μ from the observed phylogeny using the maximum likelihood method from [START_REF] Nee | The reconstructed evolutionary process[END_REF].

2. Compute the expected CDF of branch lengths (z i 's) using either eq. 1 (for the cherry branches) or eq. 2 (for the outer branches) with λ and μ as speciation and extinction rates and T equal to the age of the root node of the observed phylogeny.

3. Compute the observed statistics W 2 and Z C . 4. Generate 9999 trees using the algorithm from Stadler (2011b) with λ and μ as speciation and extinction rates, and the number of species (N) and the age of the root node (T ) both from the observed phylogeny.

5. Compute the statistic (W 2 and Z C ) for the 9999 trees. These give the null distributions of these tests. The P-values are the numbers of these 9999 values which are greater than the observed values computed at step 3 divided by 10,000 (since the observed values are considered as following the null hypothesis).

Note that in this procedure the tested model is that λ and µ are constant through time. If a time-dependent model were tested, it would be necessary to use the appropriate estimators of λ(t) and µ(t) in step 2 and use them in the simulation in step 4. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 5 10 15 20 0 5 10 15 Observed Predicted Branch length (Ma) Density 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 b) q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q 0 10 20 30 Observed Predicted

Fig. 1 .

 1 Fig. 1. (a) A tree with N = 7, (b) its branching times (or node times, or node depths, or node heights) ordered in decreasing size, (c) its branching intervals (or coalescent intervals) ordered through time, and (d) its branch lengths ordered with the seven terminal branches first (note the cherry ones which are duplicated) and then the five internal ones.

Fig. 2 .

 2 Fig. 2. A tree showing cherry and outer terminal branches.

Fig. 3 .

 3 Fig. 3. The probability densities of the number of singletons (α) and of clusters (ω) in a coalescent tree of size N = 10. All graphs are on the same scale. The column on the right-hand side gives the different coalescent events (i).

Fig. 4 .

 4 Fig. 4. Expected numbers of surviving (S t ) and extinct (E t ) lineages through time for a phylogeny with N = 200 species, speciation rate λ = 0.1 and extinction rate µ = 0.05.

Fig. 5 .

 5 Fig. 5. Distributions of cherry (left) and outer (right) terminal branch lengths. The histograms show the observed distributions for two trees simulated with N = 200, T = 50, and (a) λ = 0.1, µ = 0, or (b) λ = 0.2, µ = 0.1. The curves show the predicted distributions for each type of branch and each combination of parameters.

Fig. 6 .

 6 Fig. 6. Same as Fig. 5 except that N = 2000.

Fig. 7 .

 7 Fig. 7. (a) Distribution of cherry branches (histogram), predicted density (curve) under a model of constant speciation and extinction rates, and QQ-plot of the observed and expected quantiles for the Carnivora phylogeny (Nyakatura and Bininda-Emonds, 2012). (b) Same thing for outer branches.

Fig. 8 .

 8 Fig.8. Observed (empirical) and predicted CDFs for both types of terminal branches for the Carnivora phylogeny[START_REF] Nyakatura | Updating the evolutionary history of Carnivora (Mammalia): a new species-level supertree complete with divergence time estimates[END_REF].

Figure 1 :Figure 3 :Figure 4 :

 134 Figure1: (a) A tree with N = 7, (b) its branching times (or node times, or node depths, or node heights) ordered in decreasing size, (c) its branching intervals (or coalescent intervals) ordered through time, and (d) its branch lengths ordered with the seven terminal branches first (note the cherry ones which are duplicated) and then the five internal ones.

Figure 5 :

 5 Figure 5: Distributions of cherry (left) and outer (right) terminal branch lengths. The histograms show the observed distributions for two trees simulated with N = 200, T = 50, and (a) λ = 0.1, µ = 0, or (b) λ = 0.2, µ = 0.1. The curves

Figure 7 :

 7 Figure 7: (a) Distribution of cherry branches (histogram), predicted density (curve) under a model of constant speciation and extinction rates, and QQ-plot of the observed and expected quantiles for the Carnivora phylogeny (Nyakatura and Bininda-Emonds, 2012). (b) Same thing for outer branches.

Table 1 :

 1 Mean branch length of trees simulated with N = 200 and T = 50. The values within parentheses are the predicted means using the equations presented in this paper. The simulations were replicated 10,000 times.

	λ	µ	Terminal	Internal
			cherry	outer
	0.1 0	3.321 (3.333) 8.294 (8.327) 4.902 (4.934)
		0.05 3.609 (3.859) 9.630 (10.422) 5.990 (6.469)
		0.09 3.519 (4.486) 9.777 (12.248) 6.349 (7.302)
	0.2 0.1	1.936 (1.931) 5.364 (5.341) 3.890 (3.501)
		0.15 2.065 (2.137) 6.092 (6.407) 4.721 (4.004)
		0.19 2.053 (2.375) 6.288 (7.477) 5.052 (4.244)

Table 2 :

 2 Same as Table 1 except that N was random.

	λ	µ	Terminal	Internal
			cherry	outer
	0.1 0	3.334 (3.333)	8.339 (8.327) 4.933 (4.934)
		0.05 3.836 (3.859) 10.460 (10.422) 6.829 (6.469)
		0.09 4.165 (4.486) 12.228 (12.248) 8.524 (7.302)
	0.2 0.1	1.937 (1.931)	5.354 (5.341) 3.824 (3.501)
		0.15 2.129 (2.137)	6.424 (6.407) 5.123 (4.004)
		0.19 2.300 (2.375)	7.506 (7.477) 6.297 (4.244)

Table 3 :

 3 Same as Table 1 except that T was random. The predicted values were calculated conditioned on the value of N = 200.

	λ	µ	Terminal	Internal
			cherry	outer
	0.1 0	3.334 (3.332)	8.337 (8.657)	5.008 (5.191)
		0.05 3.862 (3.885) 10.697 (11.002)	7.683 (7.942)
		0.09 4.640 (4.656) 15.511 (15.686) 17.046 (17.111)
	0.2 0.1	1.930 (1.958)	5.347 (5.513)	3.854 (3.968)
		0.15 2.141 (2.166)	6.482 (6.624)	5.610 (5.725)
		0.19 2.393 (2.425)	8.506 (8.563) 11.016 (10.958)

Table 4 :

 4 Results of the goodness-of-fit tests of the constant-rate model to the Carnivora data using two types of branches (cherry and outer) from Nyakatura and Bininda-Emonds's (2012) tree. W 2 : Cramér-von Mises test; Z C : test from[START_REF] Zhang | Powerful goodness-of-fit tests based on the likelihood ratio[END_REF].

		W 2	P	Z C	P
	Cherry	18.85 0.0021	131.50 0.0014
	Outer	9.15 0.0553	40.23 0.0204

Acknowledgements

I am grateful to two anonymous reviewers for their constructive comments. This is publication ISEM 201x-xxx.

Appendix A.

A statistical test of the fit of a birth-death model to a phylogenetic tree may be performed by comparing the expected CDF with the empirical (i.e., observed) cumulative distribution function (ECDF). A traditional approach to perform such a comparison is to define discrete intervals over the possible values of the observed variables (branch length in our situation), count the number of observations falling within each interval, compute with the CDF the expected numbers for the same intervals, and compare the two tables of counts with a χ 2 -test. In practice, this approach presents some difficulties when the counts may be small which is the case in phylogenetic trees where typically long branches are less frequent than the short ones (see figures above). A more powerful approach is to use tests that compare directly the curves defined by the ECDF and the CDF. One of these tests is the Kolmogorov-Smirnov test which considers the largest difference between these two curves. Other tests use all information from both curves. Two of these tests have been used to assess the distribution of branching times [START_REF] Paradis | Testing for constant diversification rates using molecular phylogenies: a general approach based on statistical tests for goodness of fit[END_REF]: the Cramér-von Mises test and the Anderson-Darling test. Recently, some alternative new tests based on the same principle have been developed [START_REF] Zhang | Powerful goodness-of-fit tests based on the likelihood ratio[END_REF][START_REF] Esteban | New goodness-of-fit tests based on sample quantiles[END_REF][START_REF] Zhao | Some new goodness-of-fit tests based on stochastic sample quantiles[END_REF]. The When applied to the Carnivora phylogeny, the goodness-of-fit tests rejected very strongly the null hypothesis that the distribution of cherry branch lengths followed the constant-rate model (Table 4). For the outer branches, the Cramér-von Mises test was not significant (P = 0.055) and Zhang's Z C was slightly significant (P = 0.020). The plots of the predicted CDFs against the ECDFs show more clearly the contrast between the two types of branches: the two curves are remarkably different for the cherry branches while they are more similar for the outer ones (Fig. 8).

Appendix B. Supplementary material

Supplementary data associated with this article can be found in the online version

Cherry branches

Outer branches

Figure 2: A tree showing cherry and outer terminal branches. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q qq q q qq q q q qq qq q q qq q q q q qq q q q q q q q q q q q q q q q q q q qq qq q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

Cherry

Observed Predicted q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq qq qq q q qq q qq qq q q q q q qq q qq q q q q q q q q q q q q q q q qq qq q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q qq q q q qqq q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0


Emmanuel Paradis
email: emmanuel.paradis@ird.fr

Thierry Gosselin

Niklaus J Gr Ünwald

Thibaut Jombart

S T Éphanie Manel

Hilmar Lapp

Towards an integrated ecosystem of R packages for the analysis of population genetic data

Keywords: allelic data, NGS data, R, variant call format

Introduction

As a scientific field, population genetics, despite its relatively young age, occupies a central place in biology. It was considered early on, well before the foundation of molecular genetics, as the way forward to solve the forces behind the evolution of species [START_REF] Fisher | The genetical theory of natural selection[END_REF]. Today, many pressing issues in understanding or predicting biological responses are investigated through population genetics approaches, for example the adaptation of populations of animals or plants facing human-mediated global changes such as habitat destruction or climate warming [START_REF] Fordham | Better forecasts of range dynamics using genetic data[END_REF][START_REF] Merilä | Climate change, adaptation, and phenotypic plasticity: the problem and the evidence[END_REF]. Addressing these and other applied questions benefits greatly from an integrative approach that combines analyses of genetic data with geographical and/or ecological data.

The ongoing revolution in sequencing technologies has had a dramatic impact on the questions population geneticists and molecular ecologists can tackle. Genotypic data are now readily available for many loci for many individuals from many microbe, plant, or animal species [START_REF] Luikart | The power and promise of population genomics: from genotyping to genome typing[END_REF][START_REF] Ellegren | Genome sequencing and population genomics in non-model organisms[END_REF]. However, acquiring massive amounts of data without an interoperating ecosystem of tools to manipulate, explore, and analyse them properly for a wide variety of research questions results in -paraphrasing John [START_REF] Naisbitt | Megatrends: ten new directions transforming our lives[END_REF] -"drowning under data and starving from knowledge". Much progress in addressing this challenge for population genetics has arguably been facilitated through the de facto convergence on a single computer language: the open-source statistical analysis and programming platform R [START_REF] Ihaka | R: a language for data analysis and graphics[END_REF]. The main advantages of R for users as well as for developers of analytical methods in population genetics are described elsewhere in this special issue [START_REF] Kamvar | Developing educational resources for population genetics in R: an open and collaborative approach[END_REF][START_REF] Paradis | Linking genomics and population genetics with R[END_REF].

Evolutionary ecologists have been among the earliest adopters of R for analysing molecular data: ape [START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF]) was first released on CRAN in August 2002, andade4 (Dray et al. 2007) followed in December of the same year. These and other early adopters, such as adegenet, first released in April 2007 [START_REF] Jombart | adegenet: a R package for the multivariate analysis of genetic markers[END_REF], created the beginning of an ecosystem that made it increasingly attractive for population geneticists to continue enriching it with packages implementing new methods.

Recognizing the importance of nurturing this ecosystem, and that it faces challenges from its rapid yet mostly organic growth, the National Evolutionary Synthesis Center (NESCent) held the Population Genetics in R Hackathon on March 16-20, 2015, in Durham (USA) at the Center's headquarters. Hackathons are intensive coding events that bring together groups of people who normally do not meet, to work collaboratively and face-to-face on a shared objective [START_REF] Lapp | The 2006 NESCent phyloinformatics hackathon: a field report[END_REF][START_REF] Groen | Science hackathons for developing interdisciplinary research and collaborations[END_REF]. For the event held at NESCent, the objective was to address interoperability, scalability, and workflow building challenges for users and developers of R packages for population genetics. Such work usually receives little reward in the academic incentive system, and thus tends to fall behind without dedicated initiatives [START_REF] Howison | Incentives and integration in scientific software production[END_REF][START_REF] Prins | Toward effective software solutions for big biology[END_REF]. In order to promote the corresponding efforts of not only the participants of the event but also the community as a whole, and to raise awareness of the importance of a high-quality, interoperable, and well-documented ecosystem of analysis tools, this special issue "Population Genomics in R" highlights several products of the hackathon in combination with similar recent works from others in the population genetics community.

Hackathons have become increasingly popular, including in scientific computing, to facilitate a range of objectives such as resource adoption, community building, or tool innovation [START_REF] Trainer | Community code engagements: summer of code & hackathons for community building in scientific software[END_REF] Their results are outside the scope of this special issue and have been reported separately [START_REF] Trainer | How to hackathon: socio-technical tradeoffs in brief, intensive collocation[END_REF]).

Summary of special issue "Population Genomics in R"

Most papers in this special issue present new packages or significant improvements over existing ones. Two papers deviate form this theme. [START_REF] Kamvar | Developing educational resources for population genetics in R: an open and collaborative approach[END_REF] A motivation shared by several packages presented in this issue is to provide tools that better facilitate users' work, particularly handling their data, which may involve complex sampling designs. strataG provides user-friendly tools to handle allelic data from stratified populations, including the possibility to read a variety of data file formats and perform a series of different analyses [START_REF] Archer | stratag: an R package for manipulating, summarizing and analysing population genetic data[END_REF]. apex facilitates the manipulation of sequence data from multiple genes with tools to display them and explore incongruence among them [START_REF] Jombart | apex: phylogenetics with multiple genes[END_REF]. genepopedit is a collection of tools which helps in the manipulation of large multilocus molecular data sets and integrates with a variety of other packages [START_REF] Stanley | genepopedit: a simple and flexible tool for manipulating multilocus molecular data in R[END_REF].

Another difficult task users often face, especially with multilocus data sets, is data visualization, which is one of the chief objectives of two packages presented in this issue.

minotaur implements several measures of outliers calculated from high-dimensional genomic data and their visualization with a user-friendly interface [START_REF] Verity | minotaur: a platform for the analysis and visualization of multivariate results from genome scans with R Shiny[END_REF].

pophelper, which is both an R package and a web server, provides tools for the visualization of population structure, including the outputs of external applications [START_REF] Francis | pophelper: an R package and web app to analyse and visualize population structure[END_REF].

Handling or analysing (very) big data sets in a scalable manner is increasingly a challenge with the exponential growth in available data. Addressing this challenge is a common goal among several papers. [START_REF] Knaus | VcfR: an R package to manipulate and visualize VCF format data[END_REF] The many different file formats used by population genetics software has been an ongoing interoperability difficulty (discussed in, e.g., [START_REF] Lischer | PGDSpider: an automated data conversion tool for connecting population genetics and genomics programs[END_REF]. Many of the packages presented in this issue have the ability to read many different file formats [START_REF] Archer | stratag: an R package for manipulating, summarizing and analysing population genetic data[END_REF][START_REF] Francis | pophelper: an R package and web app to analyse and visualize population structure[END_REF][START_REF] Paradis | Linking genomics and population genetics with R[END_REF]. However, the general trend among R packages, and other software as well, has been to address the data exchange format issue by adopting unified file formats: FASTA for DNA sequences and VCF for genotypic data.

Testing for or quantifying natural selection in populations is one of the main objective of population genetics. Two of the papers in this special issue present packages that test for selection from genomic data, and use innovative algorithmic and implementation approaches to achieve substantial performance improvements over previous releases. pcadapt [START_REF] Luu | pcadapt :an R package to perform genome scans for selection based on principal component analysis[END_REF]) uses new algorithms for multivariate analysis in high-dimensional tables (see also [START_REF] Paradis | Linking genomics and population genetics with R[END_REF], whereas rehh [START_REF] Vitalis | REHH 2.0: a reimplementation of the R package REHH to detect positive selection from haplotype structure[END_REF]) uses multi-threading (light-weight parallelism) for the analysis of SNP data.

Data simulation under more or less complex scenarios has become an increasingly important task in population genetics. Two papers in this special issue present packages designed for making this task easier. SkeleSim simulates genetic data with a user-friendly interface to help users to set parameters or choose sample sizes, as well as tools to summarize outputs [START_REF] Hoban | SkeleSim: an extensible, general framework for population genetic simulation in R[END_REF]; it uses a similar user interface as minotaur. phylodyn provides functions to simulate phylogenies under a wide range of coalescent models including heterogeneous sampling [START_REF] Palacios | phylodyn: an R package for phylodynamic simulation and inference[END_REF]. This package also implements a variety of inference tools under Bayesian nonparametric coalecent.

Finally, [START_REF] Rousset | The summary likelihood method and its implementation in the Infusion package[END_REF] present the summary likelihood method of statistical inference as an alternative to the approximate Bayesian computation (ABC) method, which is also based on summary statistics when the full likelihood cannot be computed.

Impacts of the hackathon and this special issue

The NESCent hackathon that gave rise to many of the tools and products reported in this special issue was held with the kind of non-tangible objectives in nurturing the community and fostering collaboration whose achievement will only truly manifest in the long term. The collection of articles in this special issue will, hopefully, be a milestone for the future progress of R in population genetics.

It is obviously too soon to pinpoint any such long-term successes. In addition to the tangible outcomes reported by [START_REF] Trainer | How to hackathon: socio-technical tradeoffs in brief, intensive collocation[END_REF] and the products described in this special issue, the hackathon nonetheless had a variety of impacts, some smaller and some larger. It helped create new collaborations on open source tools for population genetics, as evidenced by the author teams of the hackathon-related papers in this special issue; it introduced tools (such as the package hierfstat) and their developers and users for the first time to collaborative code development and public version control on GitHub; and allowed participants to share knowledge and know-how. It also resulted in the revival of a previously existing but barely used community mailing list for population genetics in R (https://stat.ethz.ch/mailman/listinfo/r-sig-genetics).

It will remain to be seen whether these outcomes will have a lasting impact on population genetics as a field. However, the field will continue to face challenges -and opportunities -from the deluge of yet more massive amounts of genetic data, and these are likely to be addressed more effectively by a community well equipped to collaborate, whether across projects, institutions, or continents.

Open questions and future prospects

Population genetics is facing a number of challenges and opportunities from the next-generation sequencing revolution, both technical and scientific. On the technical side, how to efficiently deal with the massive amounts of data generated by next-generation sequencing will remain an ongoing problem. This includes the question of how data generated by different technologies are best analysed in combination. On the scientific side, the different drivers of genomic selection in different scenarios of population evolution have traditionally been investigated separately, using different genetic markers. Today's possibilites for genetic data collection offer the opportunity to assess how different portions of a genome are linked and evolve under different selective pressures in different environments. It will be particularly exciting to see how R and its packages will evolve to meet these challenges in the years to come.

  . Of note, the NESCent hackathon was among those chosen by a Carnegie Mellon University-based research group for studying how and with what success such events use different mechanisms to balance their different objectives.

  present how new internet-based development platforms enable the community to collaborate on promoting teaching of and education in population genetics methods. Paradis et al. (2017) synthesize the recent progress in analysing genomic data for population genetics, and present an overview of the available packages and how they integrate into a common programing environment.

  present vcfR, a package for handling variant call format (VCF) files[START_REF] Danecek | The variant call format and VCFtools[END_REF], including tools to read, write and visualize VCF as well as FASTA (DNA sequences) and GFF (annotation) files.[START_REF] Paradis | Linking genomics and population genetics with R[END_REF] present tools from the pegas package[START_REF] Paradis | pegas: an R package for population genetics with an integrated-modular approach[END_REF]) to scan VCF files, select the loci to read, and analyse them using basic R operations. To handle large quantities of data efficiently, both packages use optimized C/C++ code interfaced with R.[START_REF] Wringe | parallelnewhybrid: an R package for the parallelization of hybrid detection using NEWHYBRIDS[END_REF] use parallel code execution to efficiently detect hybrids from multilocus data in their package parallelnewhybrid.

  However, for the summary likelihood method of Rousset et al. the user does not need to formulate priors on the distribution of the parameters. The package Infusion provides a generic implementation of this method, and Rousset et al. illustrate its use with a coalescent model of population change.

Acknowledgements

We are grateful to the Editors of MER for accepting the idea of this special issue, the authors of each article for their contribution, Armando Geraldes for his assistance, and Shawn Narum and Karen Chambers for editorial support. The authors were participants in the hackathon, and are indebted to NESCent (NSF #EF-0905606) for hosting and supporting the event. This is publication ISEM 2016-236.