Clownfishes evolution below and above the species level

The difference between rapid morphological evolutionary changes observed in populations and the long periods of stasis detected in the fossil record has raised a decade-long debate about the exact role played by intraspecific mechanisms at the interspecific level. Although they represent different scales of the same evolutionary process, micro- and macroevolution are rarely studied together and few empirical studies have compared the rates of evolution and the selective pressures between both scales. Here, we analyse morphological, genetic and ecological traits in clownfishes at different evolutionary scales and demonstrate that the tempo of molecular and morphological evolution at the species level can be, to some extent, predicted from parameters estimated below the species level, such as the effective population size or the rate of evolution within populations. We also show that similar codons in the gene of the rhodopsin RH1, a light-sensitive receptor protein, are under positive selection at the intra and interspecific scales, suggesting that similar selective pressures are acting at both levels.


Introduction
Understanding the evolutionary process necessitates the integration of multiple biological scales that are rarely studied together. Population biologists have focused their efforts on the study of the variations in allelic frequencies and the mechanisms of evolution of quantitative traits at the population level, what we call microevolution [1][2][3]. By contrast, palaeontologists and phylogeneticists have been interested in the dynamics of diversification and the rate of phenotypic evolution above or at the level of the species, what is commonly referred to as macroevolution [4,5]. These different timescales challenge our understanding of evolutionary biology and question whether the mechanisms of microevolution could explain the rate of evolution at the macroevolutionary scale. A famous illustration of the difference in rates of evolution observed between timescales is the paradox of stasis, i.e. the discrepancy between the slow morphological evolution in the fossil record and the fast evolution of organisms measured near present [6,7].
Linking micro-and macroevolution remains one of the greatest current challenges in evolutionary biology [8][9][10][11][12][13][14] and relatively few empirical studies have provided mechanisms that can explain the two evolutionary scales; such as natural and sexual selection in stick insects [15], selection related to the beak morphology in Darwin finches [16], or sexual preferences related to the colour of cichlids fishes [17]. When comparing micro-and macroevolution dynamics, the rates of evolution are usually predicted to be slower at the macroevolutionary scale than at the microevolutionary scale [12,13]. A first explanation to this observation is that most of the genetic and morphological variation at the intraspecific level can be lost through evolutionary time due to the extinction of locally adapted populations [2]. For instance, at the molecular level, we expect that genetic polymorphism (neutral or under selection) observed at the present time can be lost before being fixed in the species; hence the rate of molecular evolution should be higher within species than among species [18]. At the morphological level, it is also possible that morphological evolution related to environmental fluctuations on short timescales may never accumulate on large timescales because of stabilizing selection [8,13]. At the same time, it has also been proposed that rates of evolution at macroevolutionary scale could potentially be inferred from rates of evolution at microevolutionary scale, as macroevolution corresponds to several rounds of microevolution [6,12,[19][20][21][22]. Some authors have even proposed that rates of evolution at the species level might be inferred from parameters estimated at the intraspecific level, such as effective population size [12]. Finally, it has been proposed that if there was a continuum of divergence from populations to the species level [23], similar selective pressures (ecological or evolutionary factors) might act at both micro-and macroevolutionary scales. All these predictions still remain to be tested empirically with datasets encompassing morphological and genetic data for a reasonable number of species and individuals within species [13].
Here, we compared the rate of evolution of clownfishes at both intra and interspecific scales. We first constructed two related datasets combining morphological, molecular and ecological features of clownfishes (Pomacentridae, genera Amphiprion and Premnas) at the species-level and within the species Amphiprion clarkii. Then, using newly generated sequences of one gene potentially under selection (RH1, implicated in dim light vision in deep sea [24,25]) and five morphological traits important for fish ecology (electronic supplementary material, S1), we tested the three following predictions: (i) rate of evolution should be accelerated within species compared to between species. (ii) Genetic and morphological evolutionary rates at the species level may be inferred with the parameters estimated below the species level, such as the effective population size or the microevolutionary rate. (iii) The major selective forces (i.e. related to water depth here) should be acting similarly at both the population and the species level.

(a) Comparison between micro-and macroevolutionary rates
Using our two datasets of clownfish species and A. clarkii individuals (figure 1), we showed that, in line with theoretical expectations [6,18], rates of molecular evolution at the micro-and macroevolutionary scales (r micro and r macro ) were substantially different, with higher evolutionary rates within species than among species. This difference was strong for the gene RH1, for which we estimated higher rates of molecular evolution within A. clarkii (median r micro ¼ 7.

(b) Predicting macroevolution from microevolution
Our analyses support the hypothesis that macroevolutionary rates can be predicted from microevolutionary rates for both molecular and morphological data (figures 2 and 3). Most of the intraspecific variation linked with neutral or locally adaptive alleles, is likely lost at the macroevolutionary timescale because such alleles are rarely fixed in a large number of populations of a species [2,19]. The reproductive isolation leading to speciation, and the founder effect likely associated with speciation, favours the fixation of these alleles and enables them to keep their imprint at macroevolutionary scales. To test this hypothesis with molecular sequences, we artificially simulated the fixation of one of the haplotypes present in the A. clarkii species by pruning all but one A. clarkii individual from the intraspecific tree and we estimated the rate of evolution along the single branch leading to A. clarkii. This subsampling procedure mimics the loss of intraspecific (neutral or non-neutral) variation of haplotypes through long evolutionary timescale [2,19]. Using this approach, we obtained similar estimates between the predicted rates of molecular evolution of RH1 estimated with one individual of A. clarkii (r macro_P ¼ 7.66 Â 10 24 For the five morphological traits, we predicted the rate of evolution at the species level using simulations based on parameters estimated at the intraspecific level, i.e. effective population size and the intraspecific variance of each trait (electronic supplementary material S2 and table S4, figure 3). For all traits, we found no significant difference between the empirical and the predicted rates of evolution across a posterior distribution of 1000 trees ( p . 0.05, electronic supplementary material, table S4) and regardless of the model selected (Ornstein-Uhlenbeck: OU or Brownian motion: BM). These results indicate that our estimations are robust to uncertainties in model testing, tree topology and branching times (electronic supplementary material, table S4). The simulations also retrieved the best fitting model previously selected on empirical data for four of the five traits (electronic supplementary material, table S4). We show with additional analyses that the small size of the species level phylogeny is not spuriously influencing our model selection (between BM and OU, pmc analysis, electronic supplementary material, figure S1), or affecting the size of the confidence intervals around the rate of evolution estimates ( fitcontinuous analysis [26], electronic supplementary material, table S5). Additional specificity tests also showed no support for a link between micro-and macroevolutionary rates when the effective population sizes are larger or smaller than those estimated from the empirical data (rates are significantly different from the empirical ones, electronic supplementary material, S2 and table S2). This result indicates that the match found between micro-and macroevolutionary rate is not due to a lack of statistical power.

(c) Comparison between selective pressures between micro-and macroevolutionary scales
At the molecular level, the McDonald-Kreitman test confirmed that the RH1 sequence was under positive selection ( p , 0.05, electronic supplementary material, table S6). At the microevolutionary scale, two amino acid sites (L154I and A299S) of the RH1 gene were found to be under strong positive selection across the posterior distribution of 1000 trees (dN/dS . 1 and Bayes empirical Bayes probability (BEB) . 0.9, table 1). At the macroevolutionary scale, the same two sites and two additional sites (S164S and L88F) showed a dN/dS ratio greater than one in the RH1 gene, with high BEB probability (greater than 0.9 for 1000 trees). Two sites (L88F and L154I) were associated with water depth at macroevolutionary scale (Wilcoxon signed-rank tests, p , 0.01). At the microevolutionary scale, the relationship between water depth and the two sites under selection (L154I and A299S) was not significant.

Discussion
For decades, researchers have faced the challenge to interpret macroevolutionary dynamics in the light of microevolutionary mechanisms. Only few studies have compared the dynamic of evolution below and above the species level [13,14]. Here, we propose an empirical attempt to compare and link evolutionary rates at both micro-and macroevolutionary scales, and we show that macroevolution can be, at least to some extent, predicted by microevolution. First, the molecular rates estimated at the microevolutionary scale were strikingly larger than those at the macroevolutionary scale, which is congruent with the trends observed in other studies [13,14] and corroborates well-known evolutionary patterns already described by the paradox of stasis [6].
Second, our results from both molecular and morphological data, suggest that parameters estimated at the intraspecific level, such as the intraspecific variance and the effective population size, can predict the rate of evolution at interspecific level. For molecular data, we used randomizations to simulate the loss of intraspecific genetic variance during the process of evolution, as most of the local adaptations of diverging populations will be deleted through introgression events or bottlenecks associated with the speciation process [2,19]. By pruning lineages and simulating the fixation of a single intraspecific variant at the microevolutionary scale, we were able to predict the rate of molecular evolution estimated for all clownfish species from the rate evolution of A. clarkii individuals. These results are consistent with the hypothesis of 'ephemeral divergence', previously proposed by Futuyma [2], which stipulates that most of the intraspecific molecular (and trait) variation can be considered as 'ephemeral' and should not be observable at higher taxonomic levels. Another link between micro-and macroevolutionary rates was proposed more recently by Hansen & Martins [12], who hypothesized that the tempo of molecular and morphological evolution of species should be directly related to parameters estimated at the intraspecific level, such as the effective population size [22]. Using simulations, we predicted rates of macroevolution based on the trait variance, the generation time and the effective population size estimated at the intraspecific scale, which was consistent with Hansen & Martins' hypothesis. Our results were robust to different models of evolution and phylogenetic uncertainties (topology and dating). Overall our analyses demonstrate that the rates of evolution of both molecular and morphological traits between species of clownfishes can be explained by parameters estimated at the intraspecific level, which suggest that likely similar mechanisms are at play across these two scales. Third, we detected positive selection in the rhodopsin RH1 gene (light-sensitive receptor protein) for the same two amino acid positions (L154I and A299S) at both micro-and macroevolutionary scales. RH1 is ubiquitous and well-studied among vertebrates [24], such as cichlids and sharks [24,27,28]. Several key amino acid positions are known to change the spectral specificity of the protein and thus allow organisms to adapt to different light intensity regimes [28]. We identified four amino acid positions of RH1 gene under positive selection among the clownfish species and two of those four were also present in the A. clarkii populations. One   of the sites positively selected at both scales (L154I) was also associated with water depth at macroevolutionary scale. We hypothesize that this association with water depth is also present at microevolutionary scale (even if we do not detect it) because in A. clarkii L154I covaries with another position (A299S) well known for being linked with light absorption [29][30][31]. We conclude that both L154I and A299S are two excellent candidates for adaptation to water depth. Our results suggest that molecular changes related to adaptation at the individual level likely cascades to the evolutionary process recorded at the macroevolutionary scale, but further investigation will require much larger sample sizes both in terms of individuals and number of polymorphic sites (e.g. such as in genome wide association studies).
Our approach, applied here on clownfishes, is exploratory and future studies will benefit from using larger molecular data sets for both evolutionary scales and expanding the scope to other traits likely related to adaptation to water depth, such as fin morphology. Our study also simplifies the process of evolution and predicts rates of macroevolution assuming that each incipient species is formed by one haplotype during speciation event. Further studies should investigate if speciation events involving a large amount of genetic diversity will affect the macroevolutionary rates on the long term. As more empirical transdisciplinary studies are expected in the coming years, the next decade will likely see a substantial improvement in the understanding of the links between the different scales of the evolutionary process. Recent literature proposes several directions for future research. First, one of the simplest ways to deal with the problems related to the micro-macroevolution boundary is to take into account intraspecific variance into the reconstruction of phenotypic evolution above the species level [32]. Second, as our results suggested that the evolutionary process from the populations to the species is a continuum (that is tractable numerically), we believe that the remaining gap in our understanding of the overall evolutionary process may not be strictly related to the dichotomy between micro-and macroevolution but rather to other scales, such as the link between individuals and populations levels inside species. To study evolution at multiple scales simultaneously, important developments have been made recently at the crossroad between macroevolution and ecology, such as the individual-based macroevolutionary models related to the unified neutral theory of biodiversity [33,34]. Third, a better understanding of the mechanisms of evolution might also be gained from the field of quantitative genetics focusing on the changes in allelic frequency or phenotype at each generation [35], or on the branching process of populations in function of ecological conditions and population size [36]. Finally, future challenges also concern the understanding of how the phenotypic variance evolves through time [23] in the adaptive landscape [37], and its potential impact on diversification.

Conclusion
Despite considerable discussion about the potential bridges between the micro-and the macroevolution in the last decades [9,12,38], relatively few studies have proposed an empirical test of theoretical predictions [13,14]. We suggest here that macroevolutionary rates might be inferred using microevolutionary rates in clownfishes. Furthermore, we also propose that the selective processes may act on the same genes (and even codons) at both micro-and macroevolutionary scales. While further studies are necessary to establish the generality of our findings across different organisms, our study on clownfishes provide insights into understanding of the genetic and phenotypic differentiation within and across species.

Material and methods (a) Field sampling, morphological and environmental data
For the microevolutionary scale, we chose A. clarkii because we had a good knowledge of its distribution from previous fieldwork [31]. For this species, we were able to collect molecular and morphological information for a total of 53 individuals in the centre of the range, in the Indonesian coral reefs (on the coasts of Bali and Manado, Sulawesi; table 2, electronic supplementary material, figure S5). The Indo-Australian Archipelago region has been previously shown as the main centre of genetic and morphological diversity in clownfishes [31], consistently with the peak of species diversity found in this region for other fishes [39].
Once the fishes were spotted during scuba dives, we caught the specimens using clove oil and nets. We recorded the water depth and clipped the anal or caudal fin of the fish for genetic analyses. We also took a picture of the left side of the body before releasing it back into the host anemone (electronic supplementary material, S3, and figure S2). We used these pictures to measure five main axes of body shape variation associated with key aspects of fish ecology: body and head ratios, peduncle factor, eye height and snout angle. These traits were selected because they correspond to standard morphological measurements related to functional aspects of the reef fish's ecology, such as locomotion, sensory abilities or feeding behaviour [40][41][42][43] (electronic supplementary material, S1).
For the macroevolutionary scale, we obtained samples from previous fieldwork (Bali, Indonesia; and Madagascar), and from loans from aquariums or research institutions [31] for one individual per species for 26 species of clownfishes including A. clarkii (87% of the 30 species of clownfishes). Mean depth for each clownfish species was calculated using the minimum and maximum depth collected in fishbase (http://www.fishbase.org/ [44]; electronic supplementary material, table S3), and our field data for A. clarkii (table 2). For morphological information, pictures of the left side of specimens were obtained from museum specimens. In order to obtain intraspecific variance on morphological traits per species, we collected pictures in museum and specialist collections for several individuals per species (total of 307 individuals, with an average of nine individuals per species; min ¼ 2 and max ¼ 31, electronic supplementary material, table S1). Our analysis should not be biased by plasticity between species given that we estimate mean phenotypes using several individuals per species and that clownfish species are morphologically distinct.

(b) Molecular data
We extracted the DNA of A. clarkii individuals using the DNeasy Blood and Tissue kit (Qiagen GmbH, Hilden, Germany), and we also used already extracted DNA of all other clownfish species obtained from a previous study [31]. Using standard protocols, we amplified fragments of the cytochrome B (cytB) and mtDNA control region (CR) for A. clarkii individuals and the rhodopsin gene (RH1) for A. clarkii individuals and all the other clownfish species (electronic supplementary material, S4). All newly generated sequences have been deposited in the GenBank database (table 2 and electronic supplementary material, table S3).  For the macroevolutionary scale, we used 1000 phylogenetic trees of all clownfish species directly from the posterior distribution of trees built with seven nuclear markers of a recently published study [31]. For the microevolutionary scale, we reconstructed the relationships between A. clarkii individuals using BEAST 2.1.3 [45] by concatenating the alignments from the cytB and CR markers. Relationships between individuals within A. clarkii were reconstructed using a coalescent prior with constant population size and a strict molecular clock with a uniform prior (electronic supplementary material, S5 and S6). The strict molecular clock was selected with the clock model selection implemented in MrBayes 3.2 [46]. We also chose the best model of substitution using phyml.test function in ape R package. We obtained a distribution of 1000 ultrametric trees of A. clarkii individuals with a root node at a relative date of 1. This distribution revealed that only the three deeper nodes were highly supported in the tree of A. clarkii individuals (Posterior probabilities ¼ 1). Recent nodes were less supported, which is expected given that they represent relationships between closely related individuals. To time-calibrate each A. clarkii tree, we used several individuals of A. clarkii also present in the dated phylogenetic tree at the species level [31]. We used the crown node of all A. clarkii individuals that was in common between the two trees to rescale all the branches of the A. clarkii tree. We combined each of the 1000 A. clarkii trees with each of the 1000 species level trees. Each combination of tree was chosen randomly ensuring that the resulting combined trees were not biased toward certain topologies. We thus obtained a timecalibrated distribution of 1000 topologies encompassing both A. clarkii individuals and all the other clownfish species. In the following analyses, micro-and macroevolutionary datasets were considered separately except in the analysis of the molecular rate of evolution where there were combined.

(d) Estimation of rates of evolution at micro-and macroevolutionary scales
We first estimated whether the rates of molecular evolution were different within and between species. We compared the rate of evolution of the RH1 gene between A. clarkii and all the other clownfish species as a proxy for molecular evolutionary rate (electronic supplementary material, S7). Because RH1 was potentially under selection, we estimated the rate of substitution by estimating the branch length of a RH1 tree, while constraining the topology with the dated phylogeny built using neutral genes (cytB and CR) using the optim.pml function of the phangorn R package [47]. We then estimated the mean rate of molecular evolution by calculating the ratio between the sum of the branch lengths of RH1 tree and the sum of the branch lengths of the tree built with neutral genes. This analysis was replicated separately over the two distributions of 1000 trees at micro-and macroevolutionary levels. At microevolutionary scale, the estimates of molecular rates of evolution were obtained directly from the 1000 trees of A. clarkii individuals. At macroevolutionary scale, the molecular rate of evolution was obtained from the 1000 species-level trees, but given that each of the 53 A. clarkii sequence may represent the A. clarkii species in the tree, we ran the analysis for the 53 possible species level trees resulting in 53 000 analyses (¼ 53 Â 1000). We then estimated if rates of evolution were significantly different between the micro-and the macroevolutionary scales by comparing the 95% confidence intervals of the rates, i.e. when the confidence intervals were not overlapping evolutionary rates were significantly different. We compared micro-and macroevolutionary rates only for molecular data, while morphological rates of evolution were not estimated at the intraspecific level given that individual variation in morphological traits may be due to other factors than evolution (e.g. plasticity, but see the supplementary analyses in electronic supplementary material, S9). At the macroevolutionary scale, the rate of evolution for each of the five morphological traits was estimated on each the 1000 species-level trees for BM and OU models using the mvMORPH R package [48]. We assigned to the A. clarkii species, for each trait, the mean trait value of the 53 A. clarkii individuals. The best fitting model was determined using the Akaike information criterion corrected for sample size (AICc). We additionally analysed the traits using the phylogenetic Monte Carlo method (from the pmc R package [49]) to assess if the phylogenies at both scales had sufficient statistical power to distinguish between BM and OU processes. We also ran the two models with fitContinuous (from the geiger R package [26]) to quantify the amount of parameter uncertainty due to the rather small size of the phylogeny. The pmc and fitContinuous analyses were both ran on a consensus tree built using Treeannotator [50] from the distribution of 1000 species-level trees.

(e) Predicting macroevolution from microevolution
We used two different approaches, one for molecular sequences and one for morphological traits, to predict the rates of evolution between clownfishes species from rates of evolution inside A. clarkii.
For RH1 molecular sequences obtained at both microevolutionary and macroevolutionary scales, we expect that a large amount of the polymorphism observed at present will not be maintained through long evolutionary timescales. Considering that each A. clarkii individual could be an haplotype that might get fixed in this species by drift, we can artificially simulate macroevolution by pruning all but one A. clarkii individual in the combined tree containing all clownfish species and A. clarkii individuals. This approach allowed us to predict the species-level rate of evolution of A. clarkii by estimating the rate on the terminal  2). We repeated the analysis 53 times, keeping each time only one of the 53 A. clarkii individuals, to account for the variability in the choice of individual selected. This molecular evolutionary rate estimated from the branch of the A. clarkii was then compared with the rate of evolution estimated from the tree containing all the species of clownfishes. For morphological traits, we used simulations to determine whether parameters estimated at the microevolutionary scale can predict the rate of evolution of the mean species trait at the macroevolutionary scale. Our simulations were based on three parameters: (1) the trait variance s 2 , (2) the generation time t, and (3) the population size N e . These parameters were either estimated from our empirical data or obtained from the literature (electronic supplementary material, S2). The simulations can be seen as long-term approximation of the approaches of Jones et al. [51,52] and Revell [19], which allow us to use normally distributed phenotypic traits at the population level to evolve macroevolutionary rates of evolution (see electronic supplementary material, S2). Our simulations consist of drawing at each generation N e individuals from a normal distribution (defined by the mean and the variance of the species trait). While the intraspecific variance (s 2 ) is assumed to be constant, the mean is obtained as the empirical mean of the trait values of the previous generation. The mean of the trait is thus evolving stochastically from one generation to another [22]. By repeating this process across many generations along the species tree, we generated changes in the species trait values at macroevolutionary timescales (see electronic supplementary material, S2 for more details). We then fitted OU and BM models on the simulated species trait using the mvMORPH package [48] and we ranked the models by AICc. We performed 1000 simulations for each trait using the distribution of 1000 species-level trees. We first compared if the best fitting model for simulations was the same than the best fitting model for the empirical data. Then, for each trait, we compared the rate of evolution (s 2 ) inferred from the 1000 simulations with the range of rate values obtained from the empirical trait (electronic supplementary material table S4, figure 2). To do this, we computed the probability that our microevolution simulations correctly predict the macroevolutionary rate, e.g. how many simulations provided a rate of evolution falling in the 95% confidence interval of empirical rate values. We also ran specificity tests, to show that our microevolutionary simulations do not predict rates of macroevolution when the population size is artificially smaller or larger than the population size estimated from A. clarkii data (electronic supplementary material, table S2).
Our approach relies on the fact that clownfishes are a very homogeneous group of fishes with species sharing a large number of morphological and behaviour traits in common, such as their symbiosis with anemone that structures their life history [53] or their very high self-recruitment rates [54]. A. clarkii as a representative of clownfishes species was chosen because its populations reflect some diversity in host usage and distribution range and it is not a specialist species with a too narrow range (such as A. chrysogaster, A. omanensis, A. mccullochi).
(f ) Comparison between selective pressures at microand macroevolutionary scales