Geography of speciation affects rate of trait divergence in haemulid fishes

1. Background

The relation between present-day patterns of biodiversity and evolutionary processes (i.e genetic divergence, speciation, extinction, range-size variation) has generated major discussion points among evolutionary biologists and ecologists [1–5]. Many studies [4–10] have documented the connection of geographical range with speciation and extinction in marine taxa, and yet marine groups are under-represented, comprising only a small fraction of all currently published studies about the patterns and processes of speciation [11].

It has been recognized that speciation works on ecological and evolutionary timescales not easily captured by any single approach [12–15]. Advancements in genetic techniques and evolutionary models have allowed biological studies to understand mechanisms and processes involved in the isolation of populations and whether these groups can become sufficiently isolated to become species or regain genetic connectivity in the marine realm [16–20]. Additionally, genetic sequence data are now used to generate robust phylogenies of marine lineages leading to an increasing interest in using these data to investigate the role of geography in marine speciation [21–25].

The relative strength of biological and geographical isolating mechanisms can vary over the course of speciation and between different modes of speciation [15]. Owing to the expected importance of character displacement and reinforcement, models of speciation predict a strong effect of geography on the evolution of traits associated with reproductive isolation. Sister species and very closely related species that occur in sympatry may experience diversifying selection on mate choice traits, patterns of resource use, and traits that underlie the ability to make use of resources [26,27]. Indeed, models of sympatric speciation rely on a linkage between assortative mating and disruptive selection on niche traits [28–30].

More generally, there are a number of reasons why speciation may be a period in lineage history with elevated rates of trait evolution [26,31,32]. Whether adapting to a new region following a dispersal event, or in the absence of gene flow following vicariance, isolation promotes adaptive adjustment to the current geographical distribution. But while palaeontologists have explored patterns of change in rates of morphological evolution through much deeper time than seen in contemporary ecology, it has proven difficult to compare the impact of speciation versus other factors in contributing to the rates of phenotypic change in recent evolutionary radiations (e.g. [33]).

In this study, we ask how the range overlap affects rates of trait evolution in sibling species of New World haemulid fishes, a prominent group of marine shorefishes that has undergone a substantial radiation on coral reefs of the Caribbean and Eastern Pacific regions. Our recent development of a complete phylogeny of New World haemulids reveals 21 sister species pairs, with estimates of the divergence time for each pair. For each sister pair, we estimate trophic overlap and the rate of divergence of trophic morphology, locomotor morphology, overall body shape, and colouration. Range overlap is strongly bimodal with eight pairs of fully sympatric sister species and 13 allopatric pairs. Given our expectation that rates of trait evolution are high during speciation, we predict and test for a negative relationship between the rate of trait divergence and time since divergence between sister species. We also ask whether rates of divergence are affected by range overlap, predicting that sympatric pairs will experience higher rates of divergence.

2. Methods

The phylogeny used in this study was pruned from a larger phylogeny of 107 haemulid taxa and 44 out-group species [34]. Thus, the present study focused on the phylogeny of the monophyletic radiation of all 61 species of haemulids that occur in the Western Atlantic and East Pacific (figure 1). This phylogeny was constructed from a dataset of three mitochondrial and four nuclear genes, and extends results from previous studies [35,36]. For methodological details on phylogenetic reconstruction, see [34]. Specimens were collected throughout their range including regions of sympatry. Preserved specimens were stored and measured at the collections of Centro Interdisciplinario de Ciencias Marinas, Centro de Investigaciones Biológicas del Noroeste, and California Academy of Sciences.

• Download figure
• Open in new tab
• Download PowerPoint

(a) Species pairs

Our approach was to measure trait divergence between sister species, and using the estimated divergence time, we calculated the rate of trait evolution. The phylogenetic analysis revealed 21 sister species pairs. For each of these 42 species, we measured geographical range, habitat preference, diet, body shape, colouration, and a set of functional morphological traits associated with feeding and locomotion.

(b) Geographical overlap

Species ranges and distribution records were obtained from published sources [13–24]. Ranges were drawn onto equal area maps using the free geographic information system (GIS) software GRASS (Geographic Resources Analysis Support System). These maps were not intended to be highly exhaustive reconstructions of species distributions; rather they are approximations, which we used to calculate species geographical overlay. Note that distributions in this context may be seen as large-scale overlaps; however, we confirmed that all sympatric species pairs identified in this study co-occur at the microhabitat scale and have been seen during our dives within a few metres of each other. Range overlap was defined as the percentage of the more restricted species range overlapped by its more widespread sister [37]. Species were coded allopatric if ranges did not overlap or if overlap was less than 20%, whereas if species ranges overlapped more than 90% then they were categorized as sympatric.

(c) Geography of speciation

We explored the ‘age–range correlation’ for the New World haemulids, using our time-calibrated phylogeny and the values of range overlap. We plotted range overlap against time since divergence for the 21 sister pairs. Range overlap for the deeper nodes combined areas of species distributions and were calculated following the range merging procedure as in Barraclough & Vogler [38]. If a significant age–range relationship is present, then younger pairs are expected to be informative on the geography of speciation. Intercept values of 0.5 or higher and a negative slope are indicative of a predominantly sympatric model. Whereas an intercept of less than 0.5 and a positive slope suggests that the most frequent mode of speciation is allopatric [39,40].

(d) Niche overlap

To estimate niche overlap between sister species, we compiled a dataset of ecological traits using published information [41–50]. We focused on two main axes of niche differentiation: habitat and diet. Habitat included nine categories of habitat substrate type. Dietary overlap was quantified from the use of 13 discrete categories for prey items. Both diet and habitat overlap were calculated as the number of items shared by the two species divided by the sum of total unique items included.

(e) Morphological and functional traits

We quantified body shape using geometric morphometric methods. We selected 23 anatomical landmarks and 13 semilandmarks, the latter describing head profile and lateral line contour. The x,y coordinates of these landmarks were collected from digital lateral-view photographs of three to five adult individuals per species. A total of 183 pictures of 61 haemulid species were used. In geometric morphometrics, landmarks are analysed assuming biological homology [51,52]. Semilandmarks were used to represent homologous curves by sets of points, establishing a geometric homology. Landmark configurations for each specimen were aligned using a generalized Procrustes superimposition and average species shape was calculated using the Geomorph [53] package. TPSrelw v. 1.5.1 was used to perform a relative warp analysis. This corresponds to a principal component analysis (PCA) of variation in shape. Finally, a matrix of scores of all species in each of the relative warps was generated.

Additionally, we analysed a set of functional traits related to feeding (11 traits) and locomotion (13 traits) [54]. All linear measurements were log transformed, while masses were cube root transformed prior to analysis. We performed a phylogenetic size correction and phylogenetic PCA [55] to reduce dimensionality of each dataset. The broken stick method was used to determine how many principal components were retained for further analysis. All dataset manipulations and statistics were performed in R (R Development Core Team 2012) using the ape [56], geiger [57], and phytools [53] libraries.

(f) Body colouration

We selected one photograph per species from which to measure colouration and digitized 73 anatomical landmarks and semilandmarks using TPSDig v. 1.2.1. Haemulids are not known to show sexual dichromatism. The landmarks captured the overall body shape of the species. We superimposed all of these images, fitted all images to the consensus configuration, and then averaged the unwarped images. This procedure resulted in identical (shape and size) images for each species. These images were modified in Photoshop by averaging pixel colours using a predetermined palette, based on dominant haemulid colours. Each pixel on the body had a numerical value corresponding to a combination of three independently stored colour channels: red, green, and blue. These three matrices were extracted and concatenated. Colour extraction and manipulation of edited images was done in R with the pixmap library.

(g) Standardized contrasts

We used phylogenetic standardized contrasts of trait values between sister species to estimate the rate of trait evolution [58,59]. Contrasts were calculated for body colouration and the first phylogenetic principal component (PC) of body shape, trophic traits, and locomotor traits using 1000 trees sampled randomly from the posterior distribution of time-calibrated phylogenies. Contrasts were calculated as the Euclidean distance between species in the trait space and standardized by dividing by the square root of the sum of the branch lengths (time) separating the species [60]. The variance of standardized contrasts is an estimate of the Brownian rate parameter for that trait on the phylogeny when the branch lengths are in units of time [59] and thus individual standardized contrasts estimate the rate of trait evolution across the section to the tree they are calculated on. We used a linear regression to test for a significant relationship between standardized contrasts and the divergence time between sister species among sympatric and allopatric pairs separately. We used ANOVA to test for overall differences between allopatric and sympatric pairs in the rate of trait evolution. To account for phylogenetic uncertainty and uncertainty in divergence times, all analyses were performed on a set of 1000 trees sampled from the posterior distribution of trees (electronic supplementary material). However, results and statistics reported here are based on the maximum clade credibility tree.

3. Results

Overall, the combined molecular data yielded a robust phylogenetic hypothesis with strong support (Bayesian posterior = 1.0) for more than two-thirds of all nodes. Divergence times between sister species exhibit a range of ages from a maximum of 17.5–10.3 Ma between Genyatremus cavifrons and G. pacifici to 1.6–0.3 Ma in Haemulon bonariense and H. parra.

(a) Geography of speciation

Twenty-one pairs of sister species were identified with only three of them being supported by less than 0.95 posterior probability. Two distinct clusters of pairs were found in degree of range overlap (figure 2). Among the 13 allopatric pairs, seven were geminate pairs across the Panamanian isthmus and two previously unidentified pairs each include one member from the Galapagos Islands and one from the Caribbean.

• Download figure
• Open in new tab
• Download PowerPoint

A regression of arcsine-transformed range overlap against time since divergence performed on the entire phylogeny was significant; however, the adjusted R² was extremely low (adjusted R² = 0.085; F-statistic: 6.385 on 1 and 57 d.f., p-value = 0.01). The intercept of 0.59 ± 0.025 s.e. was not significantly different from 0.5, with most of the recent pairs showing an even distribution between 0 and 1 in range overlap. A significant positive slope of 0.026 was found. Under an allopatric model, sympatry tends to increase with node age, suggesting that range changes have occurred with increasing times since divergence in this group. Since the intercept is not different from 0.5 (indicative of sympatry) but the slope is slightly positive (indicative of allopatry), we consider this result inconclusive on the issue of geography during species formation. However, a plot of range overlap against time of divergence resembled closely the simulated data under a null model of random ranges and 50% sympatry [9].

Among the 13 allopatric pairs, 11 showed zero range overlap, with the Panamanian isthmus being the most common geographical barrier, separating six geminate pairs. These species pairs ranged in age from 2.95 to 14.12 Ma. An interesting pattern can be seen in two Galapagos endemics each of which has a sister species in the Western Atlantic and not in the adjacent eastern Pacific as might be expected. Further discussion can be found in Tavera et al. [34].

Seven of the eight sympatric pairs are found within Haemulon. From these seven, five are broadly distributed throughout the Western Atlantic from eastern USA south to Brazil; including the Gulf of Mexico and the Caribbean and the remaining two pairs are distributed through the eastern Pacific, from Mexico to Peru. Time since divergence for these seven pairs range from 0.946 to 5.68 Ma (95% highest posterior density (HPD): 0.377–7.9682 Ma) (table 1).

The average age of the allopatric pairs (5.53 Ma; s.d. 3.94) was not significantly different from that of the sympatric pairs (4.04 Ma; s.d. 2.44) as indicated by the Welch two-sample t-test (t = 1.066, d.f. = 18.98, p-value = 0.3) in spite of a greater overall range of allopatric ages: allopatric pairs ranging from 1 to 14.2 Ma (95% HPD: 0.48–17.56 Ma) versus 1.8 to 7.8 Ma (95% HPD: 0.37–10.37 Ma) for the sympatric pairs. Results were consistent after running the test independently on 1000 trees sampled from the posterior distribution and after removing two pairs that did not have full posterior probability support (see electronic supplementary material).

(b) Geography and ecology

Diet overlap was significantly lower in sympatric pairs (mean = 0.43) than allopatric pairs (mean = 0.72) (Welch two-sample t-test; t = 3.3, d.f. = 17, p-value = 0.0038) (figure 3). However, there was no relationship between diet overlap and age of divergence within either sympatric or allopatric pairs or when the groups were combined (maximum clade credibility tree: allopatric pairs: p-value: 0.151; sympatric pairs: p-value: 0.2533). The same pattern was found for habitat overlap; however, the latter did not show significant differences between the two groups (Welch two-sample t-test; t = −0.24, d.f. = 17, p-value = 0.80) or relationship with age of divergence (maximum clade credibility tree: allopatric p-value = 0.44; sympatric p-value = 0.29) (figure 3). These results were confirmed when analyses were repeated on the 1000 trees from the posterior distribution (see electronic supplementary material, figure S3).

• Download figure
• Open in new tab
• Download PowerPoint

(c) Geography and phenotypic evolution

Contrasts for body colouration PC1, body shape PC1, locomotion PC1, and feeding PC1 all decreased significantly with divergence time in sympatric pairs (p-values: 0.00027; 0.0028; 0.039; 0.019, respectively) (figure 4). Only body colouration contrasts showed a significant relationship with divergence time in the allopatric pairs (p = 0.0032), although the trend was towards decreasing contrasts with age across allopatric pairs in all traits. Declining rates of trait evolution with time since divergence in the sympatric pairs are consistent with speciation involving a period of enhanced rates of trait evolution. The ANOVA indicated no overall difference in rates of trait evolution between sympatric and allopatric pairs for any of the trait categories.

• Download figure
• Open in new tab
• Download PowerPoint

4. Discussion

A number of major macroevolutionary issues revolve around the temporal pattern of evolutionary change through the history of lineages. For many reasons, it is expected that the rate of ecological and phenotypic evolution will not be constant through time as both the biotic and abiotic factors that cause niche and trait evolution likely change in response to a wide range of time-dependent factors. Phylogenies and comparative methods allow one to evaluate variation in the rate of trait change across speciation events in the phylogeny, but it has been difficult to evaluate changes in the pace of trait change during the time between speciation events [61]. By focusing on species pairs of varying time since divergence, we find evidence that the time period near speciation exhibits elevated rates of trait evolution in sympatric sister species of New World haemulids. The inferred rate of evolution declined with increasing time since speciation, consistent with an interpretation in which sympatric species diverge most rapidly during or immediately following speciation, and thus have lower net rates of divergence the longer they have been separated.

Although the present-day geography of sister species is not conclusive evidence of their geography during speciation, we found no overall age difference between sympatric and allopatric sister species, suggesting that our results are not an artefact of an increasing probability of range overlap among older sister species pairs [27]. The observed level of sympatry (38%) among New World grunt sister taxa was within the range reported from other animals (33–45% mammalian; 35% flycatchers; 50% Drosophila; 50% marine invertebrates) but low if compared to some other coral reef fishes (64%) and with plants (80%) [62–64]. This result is consistent with the view that allopatric speciation is frequent in the sea, even in organisms with wide dispersal capabilities [15,65]. The importance of biogeographic barriers in separating sister species in coral reef fishes has been extensively studied [63,66–73] demonstrating that allopatric speciation is an important engine in marine fishes diversification, although it appears not to be the only mode of speciation in the sea [63,74–76] as sometimes suggested [77]. Sympatric speciation, on the other hand, has always been contentious. Nevertheless, it has been demonstrated as one of the most probable modes of speciation for many marine fishes [74,76,78] and marine invertebrates [62]. We found two sympatric species pairs to be younger than or close to 1 Ma, which disagrees with previous suggestions that 2 Ma is the minimum time required for sister species to develop sympatric distributions following allopatric divergence in hermit crabs [79] and 4 Ma in other coral reef fishes [80]. Indeed, Quenouille et al. [80] found that after speciation, some species expanded their range and became sympatric, while others appear to have experienced limited or no range expansion, with the second of these dynamics being the most frequent. However, hybridization and introgression cannot be easily ruled out in these young pairs which is expected to confound attempts to reconstruct the history of divergence [81].

Shallow marine environments, such as the coral reefs and near shore habitats inhabited by haemulids, move with a frequency that is much greater than the rate of new species formation. As others have pointed out [11], global warming and cooling cycles have changed sea levels by as much as 150 m several times in the Pleistocene, meaning that the species ranges will have had to move with the moving habitats and, thus, present distribution may be a relatively poor predictor of the geography of populations at the time of species formation.

The rapid phenotypic divergence of young sympatric sister species is consistent with a role for disruptive selection and character displacement acting during time periods close to divergence. Further, the very young age of some of the sympatric sister species pairs in our study is suggestive of the possibility that these species were formed in sympatry, in which case, disruptive selection on ecological traits is clearly expected [82]. Finding this pattern in functional morphological traits involved in feeding and locomotion, as well as overall body shape and body colouration further suggests that ecological differentiation in sympatric species is most intense during and shortly after speciation. Our results differ from observations in some other groups, for example, there was a dominance of similarity in ecological and reproductive traits found in closely related but sympatric plant species within the California Floristic Province [64,83].

Interestingly, although sympatric sister species appear to undergo elevated rates of phenotypic divergence at the time of speciation, the overall rate of divergence between pairs was not generally greater in sympatric pairs than allopatric pairs. Perhaps, allopatric ranges encompass a greater diversity of environments, and so some diversity of selective pressures, which might increase adaptive responses and rates of divergence [84]. While sympatric species pairs showed less diet overlap than allopatric pairs, our study does not find strong evidence that sympatry of sister species results in elevated rates of phenotypic divergence. Allopatric pairs showed a trend towards declining rates of divergence since speciation, similar to that seen in the sympatric pairs, but this relationship was only significant in colour pattern. The combination of decreasing rates of divergence with time in sympatry and no difference in rate of divergence between allopatric and sympatric pairs suggests that ecological interactions between sister species may not be a dominant cause of rates of trait evolution in young sister pairs. Nevertheless, this conclusion is tentative because of the relatively low number of young allopatric and sympatric pairs. Indeed, the highest rates of divergence in colour pattern, the overall body shape, and feeding traits were seen in the youngest sympatric pairs, suggesting that a larger sample of these young pairs may confirm an enhanced rate of divergence in sympatry.

Our results partially agree with those of Pitteloud et al. [85] in Pyrgus butterflies in which they found higher rates of niche evolution along climatic dimensions for sister lineages that diverge in sympatry. Haemulids are an ecologically diverse group of reef and shore fishes in the tropical Western Atlantic and the Eastern Pacific. The feeding ecology of species ranges from benthic sediment sifting to mollusc crushing, piscivory, and mid-water zooplankton feeding. Trophic ecology is reflected in functional morphological features of the feeding apparatus and the locomotor system [34,54]. We find the lowest diet overlap in the youngest sympatric sister species, indicating that these lineages diverged rapidly in patterns of prey use (figure 4), in association with the highest inferred rates of divergence in feeding traits and among the highest rates in locomotion traits (figure 4). These results are consistent with an important role for ecological interactions in shaping early patterns of divergence in resource use and morphology of recently diverged sympatric sister species of New World haemulid fishes.

5. Conclusion

Although our study does not support differences in the rate of divergence between allopatric and sympatric pairs, we found that recently diverged sympatric pairs exhibit higher levels of diet divergence and higher rates of trait evolution, diverging most rapidly during or immediately following speciation. This pattern is consistent with a role for disruptive selection and character displacement during these speciation events. Whether these cases represent in situ sympatric speciation or secondary contact remains unsettled, but the results suggest that in these species pairs that are presently sympatric, particularly rapid trait evolution occurred early in their divergence.

Data accessibility

The dataset supporting this article is available in the Dryad Digital Repository: https://doi.org/10.5061/dryad.5qd0gd0 [86] and the R code has been included as part of the electronic supplementary material.

Authors' contributions

J.J.T. and P.C.W. conceived and designed the study, made fish measurements, and drafted the manuscript. J.J.T. collected specimens, carried out the molecular laboratory work, and analyses.

Competing interests

The authors have no competing interests.

Funding

This research was supported by NSF grant no. DEB-1061981 and IOS-0924489 to P.C.W. and Universidad del Valle project CI 71028 to J.J.T.

Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4397888.