An AFLP clock for absolute dating of shallow-time evolutionary history – too good to be true?

Amplified fragment length polymorphism (AFLP) is a well-established molecular method in phylogeography, shallow phylogenetics and population genetics. AFLPs are DNA-fingerprint markers that are widely distributed throughout the genome and consist largely of non-coding DNA (Meudt & Clarke 2007). Data analysis methods are continuously being developed to refine possible inferences from the dominant data produced by this method (Bonin et al. 2007; Foll et al. 2008). However, compared with DNA sequences, the classical markers of phylogeography (Avise 2000), a major drawback of AFLP is the lack of understanding of mutation rates and processes. In a recent study in Molecular Ecology, Kropf et al. (2009) proposed an ‘AFLP clock for the absolute dating of shallow-time evolutionary history’. A molecular clock based on AFLPs, which would allow placing in time recent subdivisions among phylogeographical groups, would be a big step forward and would allow answering many open questions about the recent history of species. However, reading the work of Kropf et al. (2009) in detail, it seems that a full solution to this challenge is still to come.

Kropf et al. (2009) used three data sets of alpine plants, Gentiana alpina, Kernera saxatilis and Atocion rupestris (synonym Silene rupestris), to demonstrate and calibrate the proposed clock. In each of these three species, distinct mountain phylogroups were described from each of three mountain regions (Alps and Massif Central, Pyrenees and Sierra Nevada). Assuming that the phylogroups arose postglacially, Kropf et al. (2009) used paleoecological evidence to estimate times since divergence among these groups and then relate divergence times to genetic distances. They used the slope of a linear regression to infer a rate of AFLP divergence, and proposed maximum and minimum rates to estimate a time interval for separation events. Finally, they applied the proposed rate to additional data sets from Minuartia biflora and Nigella degenii. In this comment, we evaluate the approach proposed by Kropf et al. (2009), both theoretically and empirically. First, we address the issue of linearity, an important prerequisite for the use of a genetic distance as molecular clock, and re-analyse some of the data presented by Kropf et al. (2009). We also comment on the estimation of confidence intervals. Second, we investigate the generality of the approach proposed by Kropf et al. (2009) by applying it to two additional data sets where information about divergence dates is available from macrofossils. The two data sets are from recent phylogeographical studies of Salix herbacea (Alsos et al. 2009) and Cassiope tetragona (Eidesen et al. 2007). Third, we apply the proposed rate to date the colonization of the arctic archipelago Svalbard after the last ice age for eight species (Alsos et al. 2007), as Kropf et al. (2009) did for M. biflora. In this case, a maximum divergence time is given by glacial history. These data are further used to address the AFLP clock’s sensitivity to different levels of genetic diversity in the phylogeographical groups compared. Genetic diversity within the populations in Svalbard differs indeed considerably among these eight species. We also date the separation of the Scandinavian populations from those of the Alps, and the separation between the Alps and the Pyrenees for the species that occur in these regions.

Linearity and confidence intervals

A first step in determining whether genetic distances measured from AFLP data can be used as a molecular clock is to test whether AFLP divergence is linearly correlated with time. Kropf et al. (2009) used Nei’s standard genetic distance (D_N72; Nei 1972) to quantify AFLP divergence. D_N72 estimates between all pairs of populations belonging to different mountain phylogroups were plotted against estimates of the assumed minimum time since divergence between mountain phylogroups. A linear regression of genetic divergence against time was estimated for each species and for all species together. As values resulting from pairwise comparisons are not independent, the significance of the regression slopes was tested using a Mantel test. A Mantel test does not, however, test for the linearity of an association; it only tests for a positive correlation between two matrices and ‘can give valid probability levels for any observed association’ (Mantel 1967, p. 209). It is important to distinguish between testing the significance of a regression slope, i.e. the slope is positive and significantly different from zero, and testing for linearity, i.e. the relation is best described as linear and not for example as exponential or following any other shape. From Figs 1 and 2 of Kropf et al. (2009), it can be observed that the relation between divergence and time appears in fact to be not linear in G. alpina and A. rupestris, as well as in the pooled data set. In all the three cases, the majority of D_N72 estimates between the Alps/Massif Central and the Pyrenees lay somewhat below the linear regression line, and most comparisons involving the Sierra Nevada are located above the regression line, indicating that a nonlinear relation may provide a better fit to the data.

To test for linearity, we used an anova to compare a linear model with time as a continuous variable (identical to the regression used by Kropf et al. 2009) to a model with time as a factor with three levels (t0, t1 and t2). As pairwise comparisons are not independently sampled items, the significance levels obtained in an anova are inflated. Therefore, we assessed the significance of the F-values with a permutation test: populations were randomly assigned to regions 10 000 times; the two linear models and the anova were calculated for each of the randomized datasets and the proportion of F-values equal to or larger than the observed F was taken as the significance level (p_perm). The same permutation approach was used to assess whether the parameters of the linear models were different from 0. All analyses were carried out in R 2.9.0 (R Development Core Team 2009), if not mentioned otherwise. The model with time as factor was significantly better than that with time as a linear variable for A. rupestris and for the pooled data set (Table 1), showing that the relation between divergence and time was not linear. The anova was not significant for the two other data sets, where the number of comparisons was considerably lower. The data of Kropf et al. (2009) provide only three time points, one of them being 0, to assess linearity, and thus limit the generality of possible inference. Nevertheless, our analysis revealed that the results of Kropf et al. (2009), all together, do not support the linearity of the relationship between the degree of AFLP divergence and time of isolation.

Looking at the literature, it seems that the linearity of relations of genetic divergence to time is rarely tested for directly (see Kumar 2005 for a review; Espinasa & Borowsky 1998 for RAPDs). However, Beerli et al. (1996) tested parametrically for deviations from linearity when discussing a molecular clock for allozyme frequencies. For sequence data, the applicability of a molecular clock is usually determined by model choice procedures within the context of a phylogenetic analysis (e.g. Drummond et al. 2006). But also for other types of data, such as AFLP, it is important to show that a genetic distance has a linear relation to time, at least in a certain time frame, before it can be calibrated as a molecular clock and recommended to estimate divergence times.

When estimating a linear regression, the precision of the regression slope is given as a standard error (SE) calculated from the sum of squared residuals (Sokal & Rohlf 1995). Kropf et al. (2009) transformed this SE to a standard deviation (SD) by dividing it by the square root of the sample size (number of pairwise comparisons). Such a transformation is valid between the SE of sample mean values and the SD in a sample, but it is neither correct nor meaningful for the SE of a regression slope (Sokal & Rohlf 1995). Further, the divergence rate estimated from the regression slope (r = 0.037 D_N72 per 10 000 years) was converted to a time per D_N72 unit rate (as 1/r = 0.27 million years Myr/D_N72). To estimate a maximum and minimum for this rate, the ‘SD’ of the regression slope (inferred as mentioned above) was inverted in the same way as the rate. This inverted ‘SD’ was then added to the estimated rate in time per D_N72 unit to provide a maximum: ‘slow rate’ = 1/r + 1/SD = 0.517 Myr/D_N72, respectively subtracted from it to provide a minimum: ‘fast rate’ = 1/r – 1/SD = 0.024 Myr/D_N72. Inverting the SD and the rate separately, before adding or subtracting them to create a confidence interval (CI), is not correct, because the resulting CI will be inversely proportional to the magnitude of the original SD. All together, the CI obtained by Kropf et al. (2009) is, in our opinion, much too wide. The 95% CI estimated from the regression slope in a conventional way (±1.96*SE) represents an underestimation of the real CI, because the points used to estimate the regression are not independent:

95% CI: 0.0291–0.0447, resulting in 0.22–0.34 Myr/D_N72.

An ad hoc solution may be to consider larger CI, such as

99% CI: 0.0266–0.0473, resulting in 0.21–0.38 Myr/D_N72.

A more correct approach would be to estimate a CI by bootstrapping either over populations or over AFLP loci. Note that inverting the CI of a rate always results in an asymmetric interval around the inverted rate.

Genetic divergence against time in Salix herbacea and Cassiope tetragona

We further investigated the potential of the approach proposed by Kropf et al. (2009) by applying it to two data sets, where estimates of divergence times were available from macrofossil findings. As with Kropf et al. (2009) in the case of Minuartia biflora, we assumed that populations became rapidly isolated after colonization of remote areas such as islands.

In the amphi-Atlantic S. herbacea, five distinct main phylogeographical groups were identified: E Canada/W Greenland, the Pyrenees, the Alps/Carpathians, one E Atlantic and one W Atlantic group. Four of the main groups (except the Pyrenees) were further subdivided into regional subgroups (Alsos et al. 2009). Macrofossil data were available to date approximately the divergence between the following groups: The postglacial isolation of S Scandinavia from the Alps was estimated to be 12 500 calibrated years before present (cal yr), as the most recent macrofossil records from the area between the current populations are dated to 12 160–12 795 cal yr and 12 400–12 880 cal yr (Denmark) and 13 400–12 500 cal yr (Germany/Baltic Sea) (Alsos et al. 2009 and references therein). It is likely that both Svalbard and Iceland were colonized from N Fennoscandia/Russia. The divergence between N Fennoscandia/Russia and Svalbard was estimated to 7900 cal yr, the date of the oldest recorded fossil from Svalbard (Birks 1991). The oldest fossils found in Iceland are dated 10 200–9300 cal yr (divergence dated to 9800 cal yr; Alsos et al. 2009 and references therein). The divergence between W Greenland and NE Canada (the likely source for colonization of W Greenland) was dated to 9000 cal yr (Alsos et al. 2009 and references therein). We also dated the divergence between the Alps and the Pyrenees to 15 800 years, the estimate used by Kropf et al. (2009), although S. herbacea might have been isolated in the Pyrenees for a longer time (Alsos et al. 2009).

The genetic structure in the circumpolar C. tetragona ssp. tetragona showed a strong east–west trend, with a partition in five main phylogeographical groups (the Siberian group, the Beringian group, the Canadian group, the E Canadian/W Greenlandic group and the E Greenlandic/Scandinavian group; Eidesen et al. 2007). The genetic pattern, together with glacial history and taxonomy of this species, supports Beringia as the main source region for (re-) colonization. However, the strong genetic differentiation between the Siberian group and the remainder suggests that the last westward expansion from Beringia must pre-date the last glacial maximum. A separate Siberian refugium, at least during the last glaciation, is supported by both glacial history and fossil evidence (several fossil finds of C. tetragona are reported from E Siberia, the oldest dated to 58 400 ¹⁴C yr and from Taimyr ca 27 000 ¹⁴C yr; Kienast et al. 2001; Eidesen et al. 2007 and references therein). We therefore assume that the Beringian group and the Siberian group must have been separated for at least 60 000 years. The last expansion from Beringia eastwards, into Canada, Greenland and Scandinavia, was probably postglacial. The genetic pattern indicates that migration occurred through Canada along a northern route, and then southwards along both coasts of Greenland. The oldest postglacial fossil finds in this area are from Ellesmere Island (8500 cal yr), Northwest Greenland (8000 cal yr) and Northeast Greenland (8200 cal yr; Eidesen et al. 2007 and references therein), suggesting that the split between the E Canadian/W Greenlandic group and the E Greenlandic/Scandinavian group occurred at least 8000 years ago.

D _N72 was calculated between all pairs of populations belonging to different phylogeographical groups for which time since divergence could be estimated (see Supporting Information, Appendix S1 for a list of populations). We used the software AFLPsurv (Vekemans 2002) instead of TFPGA version 1.3 (Miller 1997) used by Kropf et al. (2009), because our data sets were too large for TFPGA. D_N72 was estimated from allele frequencies calculated with the square root method (Lynch & Milligan 1994) assuming Hardy–Weinberg equilibrium. In addition, as with Kropf et al. (2009), we calculated D_N72 between all population pairs within phylogeographical groups to estimate D_N72 at time 0. D_N72 values were plotted against time and a linear regression was calculated for each species.

In S. herbacea, there was no monotonous increase of D_N72 with divergence time (Fig. 1). D_N72 between the Alps and S Scandinavia was lower than the estimates for splits attributed to more recent dates and D_N72 between the Alps and the Pyrenees was much larger than that expected from a linear relationship. Some of this variation may be a result of the fact that in this example, minimum divergence times (Alps-Pyrenees and Alp-S Scandinavia) were mixed with maximum divergence times (colonization of formerly glaciated areas dated from the oldest fossil appearances). The estimated regression slope was 0.0243 per 10 000 years (SE = 0.0026, p_perm < 0.001) and 0.0169 per 10 000 years (SE = 0.0026, p_perm < 0.001), when excluding the comparison with the Pyrenees, thus rather close to the rates estimated by Kropf et al. (2009). For C. tetragona, the regression slope was 0.0086 per 10 000 years (SE = 0.0004, p_perm < 0.001), which was considerably lower than Kropf et al. (2009)’s rate. Although for C. tetragona the relation looked rather linear (Fig. 1), a model with time as a factor was significantly better (F = 47.54, p_perm < 0.001) than a model with time as a continuous variable, indicating nonlinearity in this case also.

These results suggest that the increase of D_N72 with time is in general not linear. They also show that the proposed rate, which fitted rather well in the examples used by Kropf et al. (2009), is far from universal – a possibility Kropf et al. (2009) were well aware of. In addition to variation because of inexact time estimates and violations of the assumption of constant population sizes, the lower rates estimated in this study may be a result of the longer generation times of the dwarf shrubs C. tetragona and S. herbacea compared with that of the herbaceous plants used by Kropf et al. (2009), resulting in slower genetic drift and less mutations per time unit.

Divergence rate and genetic diversity

To further examine the performance of the proposed AFLP clock and to address its sensitivity to differences in levels of genetic diversity, we applied it to eight of the nine data sets analysed in a comparative study addressing immigration to Svalbard (Alsos et al. 2007). Five of these species are among the most thermophilous species in Svalbard; thus in situ glacial survival could be excluded. It is likely that the majority of the thermophilous species colonized Svalbard during the warm period of the Holocene, which lasted approximately from 9500 to 4000 years ago (Alsos et al. 2007 and references therein). In some of the species, genetic diversity in Svalbard was much lower than that in other regions, whereas in others the diversity levels were similar. We did not include the ninth species, Saxifraga rivularis, in this analysis, because it was difficult to delimit potential source regions for the colonization of Svalbard in this species, and because we could not exclude glacial survival in Svalbard (KB Westergaard, Tromsø University Museum, Tromsø, unpublished). Divergence time was estimated between the populations in Svalbard and those in the region from where postglacial colonization most probably took place. In addition, we estimated divergence time between the Alps and Scandinavia, and between the Alps and the Pyrenees for species occurring in these regions. As with Kropf et al. (2009), we used TFPGA to estimate D_N72 between the regions.

D _N72 values between 0.0006 and 0.107 were obtained between Svalbard and the different source regions (Table 2). The very low value for Arabis alpina is explained by the fact that in this species, all individuals sampled in formerly glaciated areas were nearly identical and there was virtually no genetic diversity in the whole region (Ehrich et al. 2007). In Vaccinium uliginosum, three different populations in Svalbard are likely to have originated from two different sources (Alsos et al. 2007). Therefore, two divergence values were estimated. In the population originating from Greenland, only four different clones were sampled. This small sample size may explain the large D_N72 estimated. The calculation of D_N72 is indeed based on estimating allele frequencies from the AFLP pattern with the square root method. To obtain reliable estimates from this approach, rather large sample sizes are required (Bonin et al. 2007). Excluding A. alpina and the Svalbard–Greenland comparison for V. uliginosum, the remaining values resulted in time estimates between 4050 years in C. tetragona and 15 930 years in Rubus chamaemorus. Time estimates varied thus almost by a factor four. For some species, times were close to what is indicated by the fossil record (S. herbacea and D. octopetala: the oldest fossil in Svalbard is 7900 cal yr old), or by reconstructions of climate conditions (sparse arctic vegetation from approximately 9500 cal yr and more thermophilous species from 9000 cal yr). For C. tetragona, the estimated time was astonishingly short, whereas the large value obtained for R. chamaemorus was clearly wrong, considering that Svalbard was still deeply glaciated 15 000 years ago and that R. chamaemorus is one of the most thermophilous species growing in the archipelago today.

Divergence times estimated between the other regions varied even more. The oldest divergence times were inferred in A. alpina (Table 2), whereas the shortest divergence time was estimated between the Alps and Scandinavia in Betula nana with 3510 years. Assuming that the alpine populations of B. nana have been isolated from the populations in Scandinavia since the end of the last glaciation in central Europe (approximately 10 000 years), even the time estimated with the ‘slow rate’ of Kropf et al. (2009), 6720 years, was considerably lower than that expected. However, it is possible that there has been some gene flow among these areas in the early Holocene through surviving patches of B. nana in Europe, north of the Alps.

From the comparisons between these eight species, it is evident that the level of intra-population diversity influences the estimates of D_N72. Genetic diversity within populations was in general low in A. alpina compared with that in other species (Ehrich et al. 2007), and in this species, the longest divergence times overall were estimated. In S. herbacea, the divergence between the Pyrenees and the Alps was particularly large (Fig. 1), and the population from the Pyrenees was one of the populations with the lowest diversity in that data set (Alsos et al. 2009). The species with the highest time estimates between Svalbard and the source region were those with the most severe reduction in genetic diversity between the source region and Svalbard (Table 2). The relation between time and the proportion of average intra-population genetic diversity in Svalbard to that in the source region was significant (linear model, parameter estimate = −16572, SE = 4127, t = −4.016, P = 0.01; excluding A. alpina). Repeated immigration to Svalbard during the Holocene (Alsos et al. 2007) could explain the lower time estimates obtained for species with a small reduction of genetic diversity in Svalbard compared with that in the source region, but cannot explain the very large time estimates obtained in other species. Several of the populations in Svalbard, as well as S. herbacea from the Pyrenees, are small and isolated outposts of the species’ range. In small and isolated populations, as well as in populations descending from a small number of founders, genetic drift is enhanced, explaining low genetic diversity and large estimates of divergence times. In such cases, a molecular rate based on frequency differences is unlikely to be applicable.

Conclusions

Altogether our results indicate that the claimed absolute dating of shallow-time evolutionary history based on AFLP data is too good to be true. Of course, AFLP band frequencies are subject to genetic drift and, in general, AFLP divergence, quantified with any genetic distance measure, increases with time. This general relationship seems, however, not sufficient to allow the calibration of an AFLP clock. We have shown that the relation between AFLP divergence, as measured by D_N72, and time is not linear, both for some of the data sets presented by Kropf et al. (2009) and for two additional data sets. Although the large interval between the ‘fast rate’ and the ‘slow rate’, together with the uncertainty of phylogeographical histories, makes it difficult to find examples where estimated dates are not plausible, we obtained several surprising results when we applied the proposed rate to eight data sets from a comparative phylogeographical study. Our analyses have also shown that the suggested genetic measure, D_N72, is clearly influenced by intra-population genetic diversity, biasing inferred divergence times.

Considering the processes leading to genetic differentiation between two isolated phylogeographical groups, including both the complexity of real population histories and the molecular mechanisms underlying the variation revealed by the AFLP method, the proposed approach converting a summary statistic to a divergence time is clearly an oversimplification. In the future, it would be interesting to further investigate AFLP mutation processes and rates, and to attempt combining such knowledge with phylogeographical modelling.