How speciation mechanisms shape biodiversity on Earth is a long-standing question in biology [1–6]. Because speciation depends on the response of species to abiotic conditions and biotic interactions, understanding diversification processes relates directly to species ecology [7]. While divergent selection pressures are considered as the major drivers of lineage differentiation, the importance of shifts in ecological niche may vary according to the geographical mode of divergence. Show In sympatric speciation, reproductive isolation is typically caused by differential adaptation of populations to contrasting environmental conditions [8–11], resulting in ecological speciation [12–15]. Sympatric speciation thus relies on adaptation to distinct ecological niches, which creates barriers to gene flow among populations and so allows lineage divergence [5,16]. In contrast, it is often assumed that the geographical isolation of populations in allopatry enables genetic drift to produce differentiation of sister lineages [17,18]. However, even though allopatric speciation has often been shown to be neutral, it can potentially be associated with ecological divergence [19]—driven by natural selection—when the geographical areas occupied by the diverging lineages show distinct environments [20–22]. Since they are based on different geographical contexts, sympatric speciation is also considered to differ from allopatry because it can directly arise as a consequence of competition for resources [23]. This hypothesis is supported by several models [24], as well as by cases of recently diverging lineages [25,26] and ongoing ecological speciation taking place under natural conditions [27–29]. Over small spatial scales, habitat heterogeneity in the form of climatic variation or differences in community composition may initiate niche shifting in response to competition. Variation in the performance of populations under contrasting environmental conditions may ultimately promote divergence, as shown along a depth gradient in fishes [30]. In contrast to sympatric processes, many events of speciation in allopatry are not associated with any form of ecological divergence [31], resulting in a strong signal of niche conservatism [32,33]. Even when speciation is partly driven by ecological divergence, ecological niche shifts are nonetheless expected to be less frequent than in sympatry. Because physical separation does not allow the differential action of competition, this driver of divergence is absent in allopatry. The alternative view would be that in allopatry we would expect a greater possibility of distinct opportunities available, as we would probably expect both abiotic and biotic conditions to be more distinct than in sympatry. While ecological niche shifts have been demonstrated in many studies using single-taxon approaches, the examination of the signature of ecological speciation at larger evolutionary and geographical scales remains sparse in the literature [34–37]. This study aims at improving our understanding of the geographical context of speciation associated with ecological niche shifts, at a genus-wide evolutionary level and over large spatial scales. Using the Palaearctic butterfly genus Pyrgus, we test the hypothesis that the evolutionary rate of the abiotic component of the ecological niche is higher in cases of sympatric versus allopatric speciation. We investigate variation in the evolutionary rates of eight climatic variables for sympatric and allopatric speciation events using increasingly larger delineations of geographical areas (e.g. mountain massifs) and comparing strict and broad definitions of speciation modes to test the influence of geographical areas and speciation modes. We propose an integrative approach combining next-generation sequencing of the mitochondrial genome and ribosomal DNA, with phylogenetic, biogeographic and ecological analyses. We focus on an almost exhaustive sampling of the genus of skipper butterflies Pyrgus (Hübner, 1819; family: Hesperiidae, subfamily: Pyrginae, tribe: Pyrgini; [38]), represented by ca 37 species in the Palaearctic. Climatic factors associated with temperature and precipitation are critical determinants of ecological niche establishment in insects [39,40]. Thus, it is expected that such ecological factors should be associated with differences in butterfly life-history traits such as growth, development and survival at adult or pre-imaginal stages [41,42]. Whereas biotic components related with host-plant use are typically considered as key ecological factors to explain lineage divergence in phytophagous insects [43,44], this is not the case in Pyrgus butterflies as most of the species are generalists on Potentilla spp. or on various host plants [45]. Here, we targeted all divergent pairs of sister lineages (more than 90% of the Palaearctic taxa being included in our study), that evolved within contrasting geographical contexts, to test if evolutionary rates of the climatic niche are greater when divergence produces lineages with spatially overlapping ranges, than when divergence occurs between regions. Sampling includes 36 of the 37 accepted Pyrgus species in the Palaearctic. Specimens were obtained through either field collection, museums or private donations. Preservation conditions of samples varied according to their origin (see electronic supplementary material, appendix S1, table S1). Specimens from museums or private collections were all dry-pinned and collected between 1929 and 2012. Field samples were collected from 1993 to 2013 and stored in 90% EtOH at −18°C. The sampling was broadened with published and unpublished COI sequences from additional specimens to constitute a final molecular dataset of 132 specimens. For all the studied species, we selected several specimens from different locations across their distribution. Spialia sertorius (Hoffmannsegg, 1804, tribe: Carcharodini), Erynnis montanus (Bremer, 1861, tribe: Erynnini, GenBank accession number KC659955), Eagris sabadius (Gray, 1832, tribe: Tagiadini, GenBank accession number FJ817723), and Celaenorrhinus humbloti (Mabille, 1884, tribe: Celaenorrhinini, GenBank accession number FJ817832.1) were used as outgroup taxa in the phylogenetic reconstructions [46]. Species identification was performed following de Jong [47], based on morphology and upon examination of genitalia when necessary, as Pyrgus species are relatively undifferentiated in overall appearance (figure 1a).  Figure 1. (a) Illustration of Pyrgus species. (i) Pyrgus carthami, (ii) Pyrgus malvae, (iii) Pyrgus alveus, (iv) Pyrgusfoulquieri. Most Pyrgus species are only characterized by slight variations in wing patterns. Genitalia analysis is often required to confirm morphological determination, which nonetheless remains particularly difficult within the alveus species complex characterized by several taxonomic ambiguities. (b) Chronogram with tips cartooned to the species level, displaying phylogenetic relationships among Pyrgus taxa. Colours on branches indicate inferred sympatric (pink) versus allopatric (green) speciation events used for the climatic niche analyses in the sensu lato speciation mode assignment considering 15 biogeographic areas. Node supports higher than 0.5 are displayed above branches.
Gathering of species occurrences required for biogeographic and ecological niche analyses was carried out based on exhaustive species’ range distributions. We used all available precise coordinates from the samples included in this study, in combination with occurrence data from the Butterfly Tissue and DNA Collection at the Butterfly Diversity & Evolution laboratory (Institute of Evolutionary Biology, Barcelona, Spain) and from the Global Biodiversity Information Facility [48]. For taxa distributed in Asia, we used all known accessible databases and monographies (see electronic supplementary material, appendix S1, figure S1). We considered only occurrences with a precision of at least 5 km (see below), yielding a total of 16 177 records with a median value of 146 occurrence points per species. DNA was extracted from butterfly legs. Museum specimens were extracted using the QIAamp DNA Micro kit (Qiagen, Hombrechtikon, Switzerland) whereas extraction of fresh samples was automated using a high-throughput DNA-extraction robot (Biosprint 96 workstation—Qiagen, Hombrechtikon, Switzerland). Mitochondrial sequences of nine collection samples were obtained by direct multiplex sequencing (DMPS) following a modified protocol of Stiller et al. [49] and Meyer & Kircher [50]. A detailed description of the procedure is available in electronic supplementary material, appendix S2. Sequencing was performed with an Illumina MiSeq 2 × 150 bp protocol (Lausanne Genomic Technologies Facility, Switzerland). Mitochondrial genomes and nuclear ribosomal DNA regions were analysed in 56 specimens: shotgun libraries for 24 fresh Pyrgus specimens were prepared using the Nextera DNA Sample-Prep-kit (Illumina, San Diego, USA) and further sequenced by Illumina HiSeq 2 × 100 bp (Fasteris, Geneva, Switzerland); 32 Asian and rare European specimens, mostly from private collections, were shotgun-sequenced using a modified version of Tin et al. [51]—the main modification includes a phosphorylation of the denatured extracted DNA allowing the ligation of the same P1 adaptors used for the DMPS technique (electronic supplementary material, appendix S2)—and sequenced by Illumina MiSeq 2 × 300 bp (Lausanne Genomic Technologies Facility, Switzerland). As historical shotgun libraries are more prone to sequencing errors (historical samples are usually characterized by fragmented DNA, i.e. less than 300 bp [52,53]), a large overlap between paired sequenced reads allowed for more accurate sequence assessment. Both shotgun and DMPS raw reads were demultiplexed and cleaned with Trimmomatic [54]. DMPS libraries were assembled with MIRA 4 [55], using the complete mitochondrial genome of the closely related species Erynnis montanus (GenBank accession number KC659955) as a reference. For each specimen, contigs were aligned independently using CLUSTAL W [56]. Alignments were merged together with COI sequences by pairwise profile alignment and manually checked with BioEdit 7.0.4.1 [57]. For HiSeq and MiSeq shotgun libraries, de novo assemblies were performed using SPAdes [58]. The nuclear reference required for alignment was built using the longest SPAdes contigs that matched along ribosomal DNA using BLAST [59]; the latter included 18S ITS1, 5.8S, ITS2 and 28S regions. The Erynnis montanus complete mitochondrial genome served again as reference for mitochondrial alignment. Contigs were individually mapped against references using the Geneious mapper algorithm [60] with the medium sensitivity option set with two iterations. Multiple alignment was performed using MAFFT 7.017 [61] with the G-INS-I algorithm, 200PAM/k = 2 scoring matrix, gap open penalty of 1.53 and offset value of 0.123. The mitochondrial genomes as well as the nuclear ribosomal DNA regions were aligned and subsequently annotated using respectively the complete and fully annotated sequences of the mitochondrial genome (mtDNA) of Erynnis montanus (GenBank number KC659955) and the nuclear ribosomal DNA region of Papilio xuthus (GenBank number AB674749] as references. The alignment and annotation of the mitochondrial and nuclear DNA regions were done in Geneious v. 8.1.7 (Biomatters, Auckland, New Zealand) using MUSCLE [62] to conduct the alignments and the implemented Annotate and Predict functions (a similarity threshold of 80% was applied to predict and transfer the annotations). The annotations were subsequently manually checked and the reading frames adjusted if necessary. The full matrices comprising shotgun-based rDNA and mtDNA—also encompassing DMPS and COI sequences—were produced by pairwise profile alignment using ClustalX [63] and manually corrected with BioEdit 7.0.4.1 [57]. Characters of phylogenetic relevance were selected using BMGE 1.1 [64], yielding a matrix of 20 508 nucleotides—6174 bp and 14 334 bp for the ribosomal DNA and mitochondrial genome, respectively. The two aligned and fully partitioned supermatrices (mitochondrial and nuclear alignments) were concatenated and used to conduct the phylogenetic inferences. The accession of Erynnis montanus was used as the most external outgroup taxon. Partitioned maximum-likelihood (ML) analyses were conducted using RAxML [65,66]. We used the RAxML-HPC v.8 tool implemented on the CIPRES portal (http://www.phylo.org/). Analyses were fully partitioned by providing files containing the set of partitions and were run using the GTRCAT model applied to each partition, as well as 25 rate categories as suggested by Stamatakis [65]. Node support was performed using aBayes calculation in PhyML [67,68], for which the GTR + G analysis on one single partition produced a topology identical to the one produced by RAxML. Finally, because the biogeographic and ecological niche analyses require an ultrametric phylogenetic tree [69], we first performed a dated phylogeny using BEAST 1.8.3 (http://beast.bio.ed.ac.uk/) based on the COI alignment, by applying a classic mitochondrial DNA molecular clock of 2.3% divergence per million year [70]; second the inferred root date was used on the combined partitioned best RAxML tree using a penalized likelihood approach, as implemented in the R ape package [71]. The smoothing value (λ = 100) was established using the cross-validation routine implemented in the ape package and the relaxed model of substitution rate variation among branches was applied (see ape function chronos). In this study, the biogeographic scenario of Palaearctic species of Pyrgus will serve as a basis to estimate and compare the evolutionary rates between splits in allopatry versus sympatry along the climatic niche. To infer the biogeographic history of Palaearctic Pyrgus species, we defined biogeographic areas based on natural geographical boundaries and the distribution of the genus as follows. In order to examine putative effects of the spatial scale at which sympatry is defined, we considered three biogeographic divisions of the Palaearctic, ranging from narrow to broad area definitions (electronic supplementary material, appendix S1, figure S1). We first used 11 areas that were defined as follows: (A) North Africa and the Iberian peninsula (including the Pyrenees), (B) Western Europe including the Alps and surrounding territories limited by the North Sea, the Balkans and including the Italian pre-Alps, (C) Italy, Sardinia, Corsica and Sicily, (D) Eastern Europe including the Balkans and the Carpathians, (E) Anatolia, the Caucasus, Crimea, Southwestern Russia and Mesopotamia, (F) Ural Mountains, (G) Fennoscandinavia and Baltic regions, (H) Central Asia as well as Altai and Cashmere mountain massifs, (I) Central Himalayan range, Eastern Tibet and Sichuan (China), (K) Western Manchuria (China), Mongolia and Siberia (Russia), and (J) South-East Russia and Japan. A finer partitioning resulted in 15 divisions of the Palaearctic, while the broader delineation recasts areas into seven regions. For a detailed description of area partitioning in 15 and seven areas, please refer to the electronic supplementary material, appendix S1, figure S1. Biogeographic scenarios were inferred using the dispersal–extinction–cladogenesis (DEC) likelihood model [72,73], extended from dispersal–vicariance analysis [74]. The inference of ancestral areas along a phylogenetic tree considers two processes: (i) anagenetic (internode) evolution and (ii) cladogenetic range evolution (at nodes) [73]. Anagenetic evolution is governed by a transition matrix (Q-matrix). Q-matrices, designed independently for each biogeographic division of the Palaearctic, were adjusted to reflect area connectivity by assigning high dispersal probabilities to 100 (=1) when areas were adjacent and lower dispersal probabilities to 10−3 (=0.001) when they were not, with a few cases where they were set to 10−1.5 (= 0.032) for intermediate cases where areas were not strictly adjacent (electronic supplementary material, appendix S1, table S2). We did not consider the intraspecific level because coalescent processes may prevent lineage sorting and the formation of well-segregated groups. Cladogenetic range processes at each node were assessed by performing DEC analysis. Analyses were performed using the RASP software [75] providing the ancestral area and its probability of assignment at each node. Based on the DEC ancestral areas reconstructions and the corresponding events route, we attributed speciation modes (i.e. sympatry versus allopatry) for each node leading to pairs of sister species (i.e. at the level of the most recent common ancestors of these sister lineages). We firstly applied a broad definition of sympatry, including dispersal events during the subsequent anagenetic process and cases of more complex geographical context (e.g. divergence events implying more than one common area between sister species). Following Buerki et al. [76], instances of peripheral isolation were registered under allopatry. A more stringent definition of sympatry was also applied by selecting only the nodes unambiguously showing a split either in sympatry or in allopatry. Climatic niche analyses were conducted for broad and strict definitions of sympatry and for the three partitioning types of biogeographic areas. Based on species occurrences gathered from all known records of Pyrgus, we extracted climatic data using environmental layers from Hijmans et al. [77] on http://www.worldclim.org, with a grid size of 5 arc-minutes resolution. We selected eight BioClim variables among the 19 available. Related to temperature and precipitation, these variables were chosen according to their ability to describe basic general features of the environment important for butterfly ecology [78,79]. Temperature was represented by the annual mean temperature (Bio1), the maximal temperature of the warmest month (Bio5), the minimal temperature of the coldest month (Bio6), the mean temperature of the warmest quarter (Bio10) and the mean temperature of the coldest quarter (Bio11). Precipitation was described by the annual precipitation (Bio12), the precipitation of the warmest quarter (Bio18) and the precipitation of the coldest quarter (Bio19). Note that we used here a species-based rather than a tip-based approach (i.e. with splits of interest selected upstream of the intraspecific level) to comply with requirements of algorithms that compute evolutionary rates along phylogenies [80] and to increase the sampling for a more accurate description of the ecological niche. Following climatic data extraction of specimen localities, we established ecological-niche evolutionary rates using phylogenies with one single-tip per species, each of them summarized by mean values for each of the eight BioClim variables. To account for branch-length uncertainty, we randomly sampled from the dated phylogenetic tree a total of 100 single-tipped-species trees with identical topology but different relative branch lengths. We compared evolutionary rates for climatic factors (each of the eight BioClim variables) along branches of those 100 trees using a maximum-likelihood method implemented in the brownieREML function of the Phytools R package [81]. This function uses restricted maximum likelihood to fit the Brownian rate variation model [69] in order to compare different rates of continuous trait evolution and test whether a single-rate model is favoured over a multiple rates model that allows multistate trait evolution. In the present analysis, the alternative model considers two evolutionary rates expected to differ between pairs of sister species that evolved in sympatry versus allopatry, as previously registered according to the biogeographic analyses. The multiple continuous traits of evolution correspond in this case to differential evolutionary rates of ecological niche and are reflected by the different climatic factors. Log-likelihood values of the estimated parameters were used to calculate the Akaike information criterion (AIC), which allowed selection of the best model (i.e. single versus multiple rate model). When AIC indicated a significantly better fit of the multiple versus simple model to the data (p value < 0.05), we extracted the values of the estimated evolutionary rates (sigma2) and identified rate shifts for both speciation modes (i.e. sympatry versus allopatry). Significance levels were assessed by randomly assigning allopatric versus sympatric speciation modes to the pairs of sister species and repeating the Brownie analysis, in order to generate evolutionary rates expected under a null-hypothesis. Based on the results obtained in those randomized datasets, the ratio of evolutionary rates in sympatry over allopatry was computed and compared to the ratio of the empirical estimates using Wilcoxon sign-rank tests and applying Bonferroni correction to account for multiple testing. To test whether there was significant climatic variation between biogeographic areas we performed a principal component analysis (PCA) on Bioclim extracted climatic data with the ade4 R package [82] and then compared the differences among the means on each of the two first PCA axes using one-way ANOVA and post-hoc Tukey test, with the biogeographic region as a grouping factor. All of the protein-coding genes (13) and rRNA regions (two: 12S and 16S) of the mtDNA genome were included in the supermatrix, but we only managed to successfully annotate 12 of the 22 tRNA regions (54.5%). The aligned nuclear ribosomal DNA matrix included the following five partitions: 28S rRNA, ITS1, 5.8S rRNA, ITS2 and 18S rRNA. The combined supermatrix thus contained 32 partitions and 20 508 characters (14 334 from the mtDNA genome and 6174 from the nuclear ribosomal region; see below). We did not find any incongruence between the mtDNA and ribosomal DNA topologies (when considering aBayes supports higher than 0.7); however, the combined partitioned inference mainly reflects the information provided by the mtDNA genome, since the inference from the nuclear ribosomal DNA region shows little support (data not shown). We obtained mtDNA and rDNA sequences for 56 of the 58 shotgun-sequenced specimens—assemblies were not satisfactory for the representatives of P. alpinus and P. centaureae. Our dataset therefore comprised samples from 34 of the 37 Palaearctic species described and consisted of 132 specimens (including five outgroups), with variable sequence lengths ranging from 121 bp to 20 508 bp. DMPS of historical and fresh samples yielded sequences from 121 bp to 7124 bp, MiSeq shotgun libraries of historical and fresh samples from 542 bp to 15 597 bp and HiSeq shotgun libraries of fresh samples from 19 397 bp to 20 508 bp. The median number of nucleotides analysed per specimen (excluding those for which only COI was available) was 11 744. The phylogeny largely retrieves relationships described by de Jong [47]. The monophyly of all species was supported with aBayes support higher than 0.7 and 30 among the 52 nodes describing among-species relationships within Pyrgus (i.e. 58%) exhibited support higher than 0.95 (figure 1b). Exceptions are: (i) nine species represented by a single specimen, (ii) seven taxa belonging to the P. alveus complex, a group with debated species circumscriptions, and (iii) P. melotis, a species sharing an identical COI sequence with specimens from P. malvae. The position of the outgroups was confirmed with a support of 0.87 and the dating analysis estimated the origin of the Pyrgus genus at 4.04 million years ago (see electronic supplementary material, appendix S1, figure S2). The Lagrange biogeographic scenario that takes into account phylogenetic and branch length uncertainty is given in the electronic supplementary material, appendix S1, table S3. Both allopatric and sympatric nodes used in this study were relatively well distributed across the tree depth. Attribution of speciation modes based on the broad definition of sympatry results in 16, 18 and 16 speciation events in sympatry and 10, 8 and 10 in allopatry for the partitioning of areas into seven, 11 and 15 regions, respectively (see an illustration of speciation mode assignment in figure 1b). Sympatric speciation events under the strict definition were inferred at 12, seven and six nodes and allopatric speciation was invoked at four, four and four nodes for biogeographic divisions into seven, 11 and 15 regions, respectively (electronic supplementary material, appendix S1, table S3). Based on the AIC criterion, analyses of the evolutionary rate of the climatic niche made with the brownieREML R function significantly favoured the multiple model that considers different rates of evolution to the single Brownian-motion model. Results of the comparison of the evolutionary rates are shown in table 1, for both definitions of speciation mode assignment and the three definitions of biogeographic areas. When evolutionary rates are inferred from area partitioning that considers 11 or 15 regions, the outcome is consistent across the majority of the climatic variables individually tested and exhibits a significant p-value in most analyses. When using the strict assignment of sympatry, the multiple model applied on Bio6 and Bio11 was not significantly different from the single model. Bio5 and Bio10 appeared to have higher evolutionary rates in allopatry than in sympatry. When using seven areas in biogeographic inferences, evolutionary rates of the ecological niche were always higher for the case of allopatric speciation.
Climatic differences between the biogeographic areas are illustrated in electronic supplementary material, appendix S1, figure S3. The two first axes of the PCA explained 60.6% and 30.0% of the variance, respectively. The first axis was most strongly correlated with the following variables: Bio1 (r = 0.896), Bio5 (r = 0.934) and Bio10 (r = 0.946). The second axis was mostly explained by Bio6 (r = 0.848) and Bio11 (r = 0.785). Results are given for the three definitions of biogeographic areas for the first PCA axis. The percentages of significant p-values from Tukey tests for pairwise comparisons of means were 42.9%, 60.0% and 72.4% for the division into seven, 11 and 15 areas, respectively. The signature of shifts along abiotic ecological dimensions during lineage divergence has largely been left unexplored at the macro-evolutionary scale (but see [34–37]). Ecological niche shift has been shown to be at work in a large number of cases of ongoing speciation. Alternatively, ecological niche assortment—an extension of Losos’ size assortment concept [83] to the multidimensional ecological niche—could also allow lineages to successfully colonize contrasted habitats in sympatry. Despite the fact that our experimental design does not investigate potential consequences of ecological niche assortment at the community level, it demonstrates that the footprint of ecological niche differentiation is visible at a wide phylogenetic scale. Here, we provide evidence, based on the analysis of Palaearctic Pyrgus butterflies, that the evolutionary rate of abiotic dimensions of the ecological niche is usually larger when lineages diverge in regional sympatry than in allopatry, at least when defining relatively narrow geographical areas at which sympatric versus allopatric speciation processes act (see table 1); an outcome that fits expectations of the competitive displacement theory [8] or the fitness advantage of specialization to distinct conditions [84]. However, our results point to the fact that the scale employed to define geographical areas for producing biogeographic scenarios may have a strong impact on the results. When considering broad area definitions (as in the case of our seven-areas dataset), allopatric events show a higher rate of evolution of the climatic niche than sympatric events. This is actually not surprising, as broader areas would less finely translate real geographical barriers. This may be confirmed by the low percentage of significant p-values in the Tukey tests for multiple means comparisons of climatic components among areas (42.9%), which is the lowest among all three datasets, indicating that with a seven-area definition, areas are less dissimilar climatically than when applying a more narrow definition (see electronic supplementary material, appendix S1, figure S3a). As a consequence, speciation events that would be classified as allopatric under a narrower definition of areas, would in that case be recorded as sympatric. Because branches of those splits usually have lower rates of evolution, an incorrect assignment will decrease the average values of evolutionary rates for sympatric speciation events, and blur the distinction between the two modes. Not only does the number of areas considered have an impact on the biogeographic scenario at work, but so does the way we define sympatry. Indeed, considering a strict definition of sympatry leads to a strong decrease in the number of nodes of interest, which may then produce erroneous results because of the stochasticity associated with small numbers of observations. Although sympatric speciation was still associated, generally, with a higher evolutionary rate of the climatic niche for both the 11-area and 15-area datasets, evolutionary rate estimates were higher in allopatry for one out of eight climatic variables for the 11-area dataset. We thus suggest that the definition of biogeographic areas and the number of nodes taken into account should be considered with caution to guarantee the robustness of the analyses and the subsequent interpretation of results. We therefore recommend that studies aiming at applying the same methods to other biota should (i) favour a fine geographical division of biogeographic areas consistent with dispersal characteristics of the study species and the geographical barriers of the study area and (ii) use phylogenies with a high number of nodes in order to reduce stochasticity and, eventually, to restrict node selection to those that display high phylogenetic support and high probability of ancestral areas assignment. In addition, biogeographic inference methods using more than a single input phylogenetic tree (e.g. Bayes–Lagrange, S-DIVA as implemented in RASP [75]) may account for uncertainty in topologies, branch lengths, and therefore in ancestral areas assignment. Notwithstanding with those potential limitations, the general pattern of our result is in line with a large number of specific case-studies [25–28,85] but also evidences the substantial role played by niche shifting in sympatry at the macro-evolutionary level. At a wider taxonomical scale, it also complements the findings from one of the few works on diversification of lepidopteran lineages, in which allopatric divergence without niche shift is shown to be commonly underestimated [86]. Although one might consider that allopatric cladogenesis may require niche shifts to accommodate the contrasting climatic conditions found in the occupied areas (see electronic supplementary material, appendix S1, figure S3b and S3c), these results reveal the importance of abiotic niche differentiation in sympatric speciation. This outcome contrasts with recent work on the geographical context of plant speciation, which has acknowledged that ecological divergence did not vary as a function of range overlap between sister lineages [23,87]. On the contrary, our hypothesis is particularly consistent with conclusions drawn by Bush & Smith [88] in a study where they evidence the role of competition in sympatry as a leading force for niche shifting. We suggest that similar mechanisms played an important role in shaping the evolutionary history of Pyrgus butterflies by promoting high rates of niche shifting in sympatric sister species. Our results represent compelling evidence that the footprints of micro-evolutionary processes driven by abiotic ecological components during sympatric versus allopatric speciation are reflected at the macro-evolutionary scale. Raw files are available from the Dryad Digital Repository at: http://dx.doi.org/10.5061/dryad.sc613 [89]. C.P., L.P., S.B. and N.Al. elaborated the study, analysed the data and wrote the manuscript. N.Ar. led sequence assembly and alignment. T.S., A.M.-Y. and L.F. designed wet laboratory procedures. R.V., V.D., J.H.-R., E.B., Y.C. and I.K. provided samples and data. All authors contributed in the form of discussions and suggestions, and approved the final manuscript. The authors declare no competing interests. This work was supported by an SNSF Professorship PP00P3_144870 awarded to N.A., by grants CGL2013-48277-P, PR2015-00305 (MINECO) and CGL2016-76322-P (AEI/FEDER, UE) to R.V. and a Marie Curie International Outgoing Fellowship within the 7th European Community Framework Programme to V.D. (project no. 625997). We are grateful to Catherine Berney for advice on molecular issues, to Yves Gonseth and Brent Emerson for insightful discussions, to the Lausanne Genomic Technologies Facility and Fasteris SA for the sequencing job, to Jindrich Majer, Jana Slancarova and Zdenek Fric for providing data and literature sources, to Russell Naisbit from Albion Science Editing for correcting and editing the manuscript, and to three anonymous reviewers for constructive comments on a previous version of this manuscript. FootnotesElectronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3726859. References
Page 2Bacterial cells secrete numerous extracellular factors to favourably modify their environment. These include hydrolytic enzymes, protective polymeric matrices for biofilm formation and biosurfactants that aid motility. The benefits of such exoproducts can accrue both to the producing cell and to neighbouring cells and are therefore termed ‘public goods’ [1]. Public goods are costly for individual cells to produce, and cooperating populations are consequently at risk of social exploitation by non-producing ‘cheats’ [1,2]. In theory, cheats can outcompete cooperators, because they do not incur the cost of public goods production, but derive benefits from the cooperation of others. Whether cooperation persists over evolutionary time in the face of the advantages of cheating is largely dependent on aspects of population structure that act to align individual interests [3]. Many cooperative behaviours seen in bacteria are regulated at the population level by cell-to-cell communication or quorum sensing (QS) systems [4–6]. Cells produce and release QS molecules to regulate the production of a range of public goods which aid in scavenging for nutrients, providing scaffolding for biofilms and facilitating motility. Because these cooperative secretions can be key determinants of successful growth or persistence, there has been considerable interest in the impact of QS on ecological competition between different genotypes or strains of bacteria [7–9]. For example, mutant genotypes which do not respond to QS molecules, and consequently produce fewer or no public goods (even though the loci that directly encode these public goods are intact), have been shown to act as social cheats both in vitro, in vivo and in biofilms [8,10–13]. In addition to regulating public goods production, QS molecules have been shown to have non-signalling effects, such as immune modulation, cytotoxicity, redox potential and iron binding [14,15]. The impact of these indirect effects by QS molecules on social competition has not previously been explored, and so here we empirically demonstrate how production of a QS molecule can alter the social landscape of a seemingly unrelated trait, siderophore production. Pseudomonas aeruginosa is a Gram-negative opportunistic pathogen which employs a multi-layered QS system to regulate a number of public goods, many of which are important for virulence [4,5]. One well-defined P. aeruginosa QS signal is the pseudomonas quinolone signal (PQS) [16]. PQS is a member of the 2-alkyl-4(1H)-quinolone family of molecules and acts as a QS molecule in the classical sense, in that it interacts with a specific receptor protein, and sets in motion a regulatory cascade leading to increased production of toxins and biofilms [16–18]. PQS also has other biological properties that are distinct from signalling: these include balancing redox reactions, aiding in competition with other species and interacting with cell membranes [17,19,20]. In addition, PQS has iron-chelating activity, though it does not act as a true siderophore, because it does not directly ferry iron into the cell [21,22]. It has therefore been suggested that PQS may act as an iron trap, aiding in the sequestration, but not in the membrane transport of iron [22]. Moving iron from either a host or the environment into the cell is often achieved by the production of dedicated iron scavenging molecules known as siderophores [23]. Pseudomonas aeruginosa produces two major siderophores, pyoverdine and pyochelin. Pyoverdine has been experimentally demonstrated to be a public good [24], which is exploitable by cheats both in vitro and in vivo [25,26]. Here, we test whether the iron-chelating properties of PQS can change the social landscape of siderophore production. We show that PQS (i) increases the production of pyoverdine and pyochelin, and consequentially decreases the fitness of siderophore producers and (ii) increases the relative fitness of siderophore cheats in co-culture with a producing strain. Our findings highlight how direct modification of the environment by one bacterial exoproduct, in this case a QS signal molecule, can indirectly affect the evolutionary dynamics of another social trait. For a rich, iron-replete growth environment, we used lysogeny broth (LB) (10 g l−1 tryptone, 5 g l−1 yeast extract, 10 g l−1 NaCL), and for an iron-limited growth environment we used casamino acids (CAA) medium (5 g l−1 CAA, 1.18 g l−1 K2HPO4.3H2O, 0.25 g l−1 MgSO4.7H2O). We prepared both media in dH2O and supplemented CAA medium with sodium bicarbonate solution to a total of 20 mM. For all experiments, we inoculated single colonies of the relevant bacterial strain into 5 ml LB and incubated at 37°C at 200 r.p.m. for 18 h. We then washed pre-cultures in the appropriate medium, corrected to an optical density of OD600 = 1.0, and inoculated experimental cultures to an initial density of OD600 = 0.01. To study the effects on siderophores of varying concentrations of iron, PQS and its precursor 2-heptyl-4-hydroxyquinoline (HHQ), we used the strain PAO1ΔpqsAH, which is defective in 2-alkyl-4(1H)-quinolone production [22]. To test whether investment in siderophores increased with added PQS or HHQ, we inoculated a washed pre-culture of PAO1ΔpqsAH into 750 µl LB medium containing varying concentrations of PQS and HHQ in microtitre plates and incubated at 37°C for 14 h. Following incubation, we measured the OD600 of resulting cultures, and then filtered the cell supernatants. We measured pyoverdine and pyochelin using excitation/emission assays (ex/em wavelengths 400 nm/460 nm for pyoverdine and 350 nm/430 nm for pyochelin [27] using a Tecan Multimode plate reader). We corrected fluorescence values by subtracting the fluorescence of a sterile medium blank and assuming a 5% leakage from the pyoverdine into the pyochelin channel, as previously described [27]. We estimated per-cell siderophore production as fluorescence (relative fluorescence units) divided by culture density (OD600). To measure the growth (fitness) of monocultures, we inoculated a washed pre-culture of either PAO1 or PAO1ΔpvdD/pchEF into 750 µl CAA medium supplemented with 20 mM NaHCO3 containing no addition, or supplementation with either 50 µM PQS, or 100 µg ml−1 (1.25 mM) transferrin in 48-well microtitre plates and incubated at 37°C for 14 h. To study the effect of PQS on the competition between siderophore producers and non-producers, we used wild-type PAO1 and a mutant that was defective in pyoverdine and pyochelin production and labelled with a constitutive luminescence marker (PAO1ΔpvdD/pchEF CTXlux). For competition assays, we pre-cultured, washed and density corrected both strains and mixed them to a ratio of approximately 99 : 1 (producer : non-producer). We incubated these in 5 ml iron-limited medium (CAA) in the presence and absence of 50 µM PQS for 24 h at 37°C with agitation at 200 r.p.m. To measure relative abundance of the strains, we plated the co-cultures before and after incubation, and counted total colonies and luminescent colonies. Relative fitness was calculated using the formula w = p1(1 − p0)/p0(1 − p1), where p0 and p1 are the proportion of non-producing mutants in the population before and after incubation, respectively [28]. The effect of PQS and HHQ supplementation on growth, and siderophore production were all analysed using the ordered heterogeneity approach [29]. This allows for the evaluation of an ordered alternative hypothesis but does not require the fitting of curves. We chose this approach because our question is about the effect of increasing concentrations of PQS but without any concern for the exact shape of these relationships. For each test, we calculated the test statistic rsPc = rS × (1 − p). p is the p-value from an ANOVA of raw data with concentration of PQS or HHQ fitted as a categorical variable, and rS is the absolute value of Spearman's rank correlation coefficient between the means of the relevant independent variable for each level of PQS or HHQ, and the concentration of PQS or HHQ. The relative fitness of a siderophore non-producer in iron-limiting conditions, and the effect of PQS supplementation on the relative fitness, were examined using t-tests. All statistical analyses were performed using R 3.0.2 [30]. To test whether PQS increased production of siderophores, we measured the amount of pyoverdine and pyochelin in cultures of a PQS-deficient mutant (PAO1ΔpqsAH) grown in LB, and supplemented with exogenous synthetic PQS at varying concentrations. We found that PQS reduced growth and increased the per-cell concentrations of pyoverdine and pyochelin in a concentration-dependent manner (figure 1a–c, ordered heterogeneity tests: growth F6,28 = 179, p < 0.001, rs = −0.964, rsPc = 0.964, p < 0.05; pyoverdine F6,28 = 1311, p < 0.001, rs = 0.643, rsPc = 0.643, p < 0.05; pyochelin F6,28 = 2222, p < 0.001, rs = 0.75, rsPc = 0.750, p < 0.05). As PQS plays a role in cell–cell communication, it is possible that reduced growth and induced iron scavenging are the result of QS-dependent regulation of gene expression. To exclude this possibility, we repeated the experiments in the presence of HHQ, the immediate precursor of PQS [31]. HHQ does not bind iron, but maintains a signalling role in cell–cell communication [22,31]. We found that increasing concentrations of HHQ did not significantly affect growth, pyochelin or pyoverdine production (figure 1d–f; ordered heterogeneity tests: growth F6,28 = 9.1, p < 0.001, rs = 0.214, rsPc = 0.214, p > 0.05; pyoverdine F6,28 = 6.6, p < 0.001, rs = 0.25, rsPc = 0.250, p > 0.05; pyochelin F6,28 = 24.4, p < 0.001, rs = −0.071, rsPc = 0.071, p > 0.05). Overall, we conclude that it is the iron-chelating activity of PQS, and not its signal function, that triggers an iron starvation response: cells increase their production of iron scavenging siderophores and either a metabolic burden or slower uptake due to the presence of a chelator leads to poorer growth. Figure 1. PQS causes iron starvation in P. aeruginosa cultures. (a) Increasing concentrations of exogenously added PQS decrease the growth of a PQS mutant (PAO1ΔpqsAH), in iron-rich conditions and increase the production of the iron scavenging molecules (b) pyoverdine (PVD) and (c) pyochelin (PCH). (d,e,f) Iron starvation effects are not seen with the addition of HHQ, the biosynthetic precursor to PQS that does not bind iron. Error bars represent the standard deviation of five independent measurements.
The production of siderophores is a social trait that can be exploited by non-producing cheats [23]. We therefore predicted that intra-specific social competition over iron would intensify in the presence of exogenous PQS, due to the greater pool of siderophores available and the concomitant cost to producer growth. First we looked at the effect of PQS on monocultures of a PAO1 wild-type and a strain defective in the production of both pyoverdine and pyochelin (PAO1ΔpvdD/pchEF) in iron-limited media. The PAO1ΔpvdD/pchEF strain reached slightly higher optical densities than PAO1 in CAA media (p < 0.001), but we found that PQS reduced the fitness of both strains which is consistent with the iron chelation effects of PQS (figure 2a). We compared the PQS effect against the effect of transferrin, a chelator previously used in iron-limited media siderophore experiments [25,32]. We found similar reductions in growth which shows that the PQS iron-chelating effect is comparable with that of transferrin (figure 2a). Consistent with existing work on the social dynamics of siderophore production, we found that the PAO1ΔpvdD/pchEF mutant functioned as a ‘cheat’ in iron-limiting conditions, having a relative fitness greater than 1 when grown in co-culture with the wild-type (figure 2b; t1,10 = 3.32, p < 0.01), although the small increase in the fitness of the mutant (figure 2a) could partially explain this finding. In line with our hypothesis, the addition of PQS significantly increased the relative fitness of the mutant (figure 2b; F1,10 = 95.4, p < 0.001). Figure 2. PQS increases the relative fitness of a siderophore cheat. (a) Monocultures of either PAO1 wild-type (WT) or a double pvdD/pchEF mutant (mutant) grown in CAA with either no supplementation or supplementation with either 50 µM PQS or 100 µg ml−1 transferrin. Error bars represent the standard deviation of five independent measurements. (b) A siderophore non-producing mutant gains a relative fitness advantage in co-culture with a siderophore producer in iron-limiting conditions. When 50 µM PQS is added to the culture this relative advantage increases due to increased siderophore output of the producer and subsequent increase in exploitation by the non-producer. The dashed line indicates the value of relative fitness (w = 1) at which both producer and non-producer have equal fitness. The box-plots indicate the median (line), the interquartile range (box) and the extreme values (whiskers).
Here we show, for the first time, that environmental modification via a QS molecule affects the selection for public goods that are not, as far as we are aware, directly regulated by QS. Specifically, we show that the iron-chelating properties of PQS lead to increased production of costly siderophores and consequently, increased relative fitness of a siderophore cheat. We found that the addition of synthetic PQS to cultures of P. aeruginosa results in a concentration-dependent decrease in bacterial fitness (growth) and an increase in the production of the siderophores pyoverdine and pyochelin (figure 1). The biosynthetic precursor of PQS, HHQ (which does not bind iron), had only a small effect on the production of pyoverdine but no effect on the production of pyochelin or on growth (figure 1). We hypothesized that this effect of PQS would enhance the relative fitness payoff of siderophore non-producing cheats in competition with the wild-type. Consistent with this hypothesis, we show that when siderophore production is increased by PQS in an iron-limited environment, this leads to an increase in the relative fitness of a cheating mutant (figure 2). The increased relative fitness of a cheat in the presence of PQS is likely due to a combination of the increased availability of siderophores to exploit, and the increased costs paid by siderophore-producing cells. Our findings complement and build upon previous work, which showed that when less iron is available to cells, this results in greater production of siderophores, and an increase in the relative fitness of cheats [33]. In previous work, the authors artificially modified iron levels in the growth medium [33]. Our work differs in that we show that direct modification of iron levels in the environment by a QS molecule can alter selection for siderophore production. Overall, our work builds upon a growing body of experimental studies exploring the complexities of cooperation in P. aeruginosa, an organism that is an excellent laboratory model for applying and testing and extending social evolution theory [1,2,6,8,11,13,34,35]. Microbes produce a diverse array of public goods, and little is known about how social traits interact with each other either directly or indirectly [36], although recent work has shown an interconnection between pyoverdine and pyochelin production [24]. Put another way, to what extent does the production of one social trait affect the social dynamics of another trait(s)? Existing examples include (i) the direct regulatory effect of communication on the production of public goods [7–9] and (ii) the genetic linkage of traits via pleiotropy [37]. Future work in this area should continue to highlight and demonstrate which traits are social in microbes [38,39], but also begin to focus efforts on how apparently discrete traits interact, and how this affects population ecology and evolution within environments. This will require experiments that reveal the fitness effects of trait linkage, and also experiments to unravel the mechanisms by which traits are linked. Given that we have shown PQS production enhances P. aeruginosa vulnerability to siderophore cheating, this suggests there are ecological and biological role(s) of PQS beyond its well-documented role as a QS signal [16–18]. 2-alkyl-4(1H)-quinolones (including HHQ) have previously been shown to be produced by several bacterial species, but to date, PQS has only been shown to be produced by P. aeruginosa [18,40], suggesting that PQS may have evolved functions distinct from signalling. One possibility previously suggested, is that PQS-bound iron associates with the bacterial envelope making it easier for dedicated siderophores to shuttle iron into the cell. This could ensure that metabolically expensive siderophores are not easily lost to other cells [22]. Such a mechanism could help to reduce siderophore cheating, but our data show that siderophore cheats flourish when PQS is present. Another role could be in ‘privatizing’ iron for P. aeruginosa when it is in competition with other bacterial species. In this case, PQS-bound iron could reduce the availability of iron for heterospecific competitors while making it easier for pyoverdine and pyochelin to transport iron into cells. Future work focusing on understanding the ecology between bacterial species could help to unravel these interactions. Data and R code for analysing the data are available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.81081 [41]. R.P. and S.P.D. designed the experiments. R.P., F.H., A.C.S. and S.E. collected and analysed the data in the laboratory. R.P., F.H., L.M. and A.C.S. conducted the data analysis. R.P., L.M., P.W. and S.P.D. wrote the paper. All authors revised and approved the manuscript. We have no competing interests. This work was supported by the Natural Environment Research Council (grant number NE/J007064/1); the Human Frontier Science Program (grant number RGY0081/2012); the Wellcome Trust (grant number 095831); and the Biotechnology and Biological Sciences Research Council (grant number BB/J014508/1). We thank Stuart West for comments on the manuscript. Footnotes© 2017 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited. References
Page 3During the past 10 000 years, humans have domesticated over 260 plant, 470 animal and 100 mushroom-forming fungal species [1–3]. Humans have modified these domesticates through diverse conscious or unconscious programmes of artificial selection that required, at least temporarily, reduced gene flow between populations of domesticates and those of their free-living progenitors [4]. Barriers to gene flow have included the isolation of domesticates in discrete garden plots and livestock pens accompanied by programmes of selective breeding, in some cases via asexual propagation, self-fertilization, and the propagation of reproductively isolated polyploid and translocation races. In a subset of cases, barriers to gene flow have also included the separation of domesticates from their free-living conspecifics by allopatry (i.e. by the transport, whether deliberate or incidental, of the domesticates to localities at the peripheries of or completely outside of their ancestral ranges) [5–9]. With the advent of genomics, historical patterns of domestication by humans are the focus of reinvigorated research [10–13]. Here, we examine phylogenetic patterns of non-human domestication of fungi by ants that may provide insights into the dynamics underlying similar processes of domestication in other animal groups such as termites, bark beetles and even bees [14–18]. Fungus-farming ‘attine’ ants (Formicidae: Myrmicinae: Attini: Attina) are a monophyletic subtribe of approximately 250 described New World species that cultivate fungi for food [19]. Although the ants are obligate symbionts, their fungal cultivars (Agaricaceae: Leucocoprineae and Pterulaceae: Pterula) vary in symbiotic commitment in a pattern that is highly correlated with ant and fungal phylogenies [20–23]. The multiple species of fungi associated with the most primitive fungus-farming ants are, so far as is known, facultative symbionts (i.e. capable of a free-living existence outside of the symbiosis). Because these ‘lower’ attine fungi have been shown to be freely interbreeding members of larger, conspecific, free-living populations, they are usually regarded as non-domesticated [19,24,25]. In contrast, the clade of fungal cultivar species associated with so-called ‘higher’ attine ants, including the well-known leaf-cutting ant genera Atta and Acromyrmex, have become polyploid, obligate symbionts, and are no longer capable of living apart from their ant hosts. Higher attine fungi represent the best-confirmed case of domestication in attine ant agriculture and they have become the subjects of intensive study including proteomic and genomic investigation [26–29]. Ant agriculture is hypothesized to have arisen in the wet forests of tropical South America approximately 55–65 million years ago [16,23,27,30–33]. Most subsequent attine agricultural evolution, including the domestication event that produced the ancestor of the higher attine cultivars, is likewise hypothesized to have occurred in South American rainforests because ant-cultivated fungi are thought to be native to such habitats and incapable of surviving elsewhere as free-living organisms [30,34–37]. This ‘out-of-the-rainforest’ hypothesis is consistent with the observation that of all South American habitats (including deserts and seasonally dry habitats), wet Neotropical forests are home to the highest diversity and abundance of species in the fungal tribe Leucocoprineae, from which the ancestral attine cultivars arose [20,24,38–41]. The hypothesis is also consistent with the fact that all known free-living close relatives of attine fungal cultivars, including conspecifics of ant-cultivated fungi, have been collected in the wet forests of Panama and Brazil [24]. The out-of-the-rainforest hypothesis presents a problem, however, for explaining the domestication event that resulted in the origin of higher attine fungi and for ant–fungus coevolution, more generally. If domestication requires reproductive isolation of the domesticate from its free-living progenitor, and if, in wet forests, an ant-cultivated fungal species is genetically connected to a larger population of a free-living fungal species, how could domestication have occurred? Here, we describe the results of a phylogenomic analysis of the fungus-farming ants, in which we used hundreds of ultra-conserved element (UCE) loci to resolve outstanding ambiguities in attine ant phylogeny. This robust phylogeny permits the identification of major evolutionary transitions on specific lineages, including the transition from lower to higher agriculture, with greatly reduced phylogenetic uncertainty. This phylogeny also allows the reconstruction of the ancestral habitats and areas in which these transitions is likely to have occurred. We find that the out-of-the-rainforest hypothesis is supported with regard to the origin of attine ant agriculture; however, contrary to expectation, we find that the transition from lower to higher agriculture is very likely to have occurred in a seasonally dry habitat (cerrado or savannah), inhospitable, at least during the dry season, to the growth of free-living populations of attine fungal cultivars. We suggest that inhospitable habitats favoured the isolation of attine cultivars over the evolutionary time spans necessary for domestication to occur. We also identify the sister group of attine ants with high confidence, providing a new target for investigations into the biological traits that promoted fungus farming. We chose a total of 119 taxa for inclusion in our study, representing a broad range of fungus-farming and non-fungus-farming ant species (electronic supplementary material, table S1). Within the fungus-farming ants, we included 78 species, covering all 16 genera and including the morphologically enigmatic taxa Mycetosoritis asper, M. explicatus and Paramycetophylax bruchi. For Cyphomyrmex and Trachymyrmex, two genera that are known to be non-monophyletic, we sampled broadly across all clades identified in a previous study [23]. Outside of the fungus-farming ants, we included 41 outgroup species, 29 from within the tribe Attini and 12 from outside of this group. We chose genus-level outgroups from within the Attini to match the extensive sampling of Ward et al. [42], missing only the genera Diaphoromyrma, Mesostruma and Talaridris. We employed the UCE approach to phylogenomics [43,44], combining target enrichment of ultra-conserved elements (UCEs) with multiplexed, next-generation sequencing. For UCE enrichment, we used an RNA bait library for Hymenoptera that targets 1510 UCE loci [44]. The laboratory protocol we used closely follows the methods reported in [44], and we provide a detailed description of the protocol in electronic supplementary material, appendix A2. The sequencing centres demultiplexed and converted raw Illumina data from BCL to FASTQ format. Starting with the FASTQ files, we performed all initial bioinformatics steps using the Phyluce v. 1.4 software package [45] and associated programs (see electronic supplementary material, table S2 for all sequencing and assembly statistics). We cleaned and trimmed raw reads using Illumiprocessor [46] and assembled contigs de novo using Trinity v. r2013-02-25 [47]. After assembly, we used several Phyluce scripts to identify and extract UCE contigs, remove potential paralogs, and add in data from two genome-enabled taxa (Atta cephalotes and Acromyrmex echinatior; see electronic supplementary material, table S3). We aligned the UCE loci using MAFFT v. 7.130b [48] and trimmed the alignments with Gblocks v. 0.91b [49,50] using reduced stringency settings. We filtered the master set of alignments for varying levels of taxon occupancy (percentage taxa required to be present in a given locus) and selected the 75% filtered alignment set as the main set for phylogenetic analysis (‘Attine-118T-F75’; see electronic supplementary material, table S4 for all matrix statistics). The Attine-118T-F75 alignment set includes 950 loci and 652 774 bp of sequence data, of which 305 858 sites are informative. See electronic supplementary material, appendix A2 for additional detail on matrix preparation and taxon occupancy filtering. Using the Attine-118T-F75 alignment set, we investigated the effects of inference method and data partitioning on results (see electronic supplementary material, appendix A2 for additional detail). Note that we excluded one taxon from this focal alignment set (Paramycetophylax bruchi) because it received few captured UCE loci (electronic supplementary material, table S2). We did, however, include this taxon in a separate, 119-taxon analysis, described below. For phylogeny estimation, we compared maximum likelihood (ML), Bayesian inference (BI) and species tree (ST) approaches. For concatenated ML analyses, we compared four partitioning schemes using RAxML v. 8 [51]: unpartitioned, partitioned by locus, partitioned by PartitionFinder v. 2.0 using the rcluster algorithm (PF; data pre-partitioned by locus) [52,53], and partitioned by PF v. 2.0 using the kmeans algorithm [54]. For each analysis, we executed a rapid bootstrap plus best tree search (‘-f a’ option), and we used the GTR + Γ model of sequence evolution (for both best tree and bootstrap searches). To address the concern that bootstrap scores can be misleading with phylogenomic data [55], we also performed a jackknifing analysis, in which we randomly sampled subsets of loci multiple times (100 replicates of 100 UCE loci). For BI analyses we used ExaBayes v. 1.4.1 [56], and performed unpartitioned and partitioned searches on the concatenated matrix. For the partitioned searches we used the same kmeans-partitioning scheme that we used with RAxML. We selected kmeans because the tree resulting from the kmeans-partitioned ML analysis had the highest likelihood score and reasonable branch length estimates (electronic supplementary material, table S5). For the partitioned BI searches, we performed two analyses, one with parsimony starting trees and one with random starting trees. We assessed run performance by examining log files with Tracer v. 1.6.0 [57]. We performed ST analyses using the program ASTRAL v. 4.8.0 [58,59]. First, we generated gene trees using RAxML and then we used only the 500 gene trees with the highest mean bootstrap scores (calculated in R v. 3.2.2 [60] using a script modified from [61]) as input into ASTRAL. We included this step to reduce noise introduced by including uninformative loci (see [62]), and we conducted the ASTRAL analysis with 200 multi-locus bootstrap replicates [63]. To explore our data for other potential biases, we generated two additional matrices. First, we used the program BaCoCa [64] to identify any loci exhibiting significant deviations from base composition heterogeneity (χ2 test, p < 0.05). After removing biased loci (36 total), we concatenated the remaining 914 loci for analysis (‘Attine-118T-F75-975’). Second, to control for either base composition heterogeneity or saturation, we converted the concatenated Attine-118T-F75 matrix to RY-coding (‘Attine-118T-F75-RY’). We analysed both of the above matrices unpartitioned with RAxML. To place Paramycetophylax bruchi (excluded due to poor UCE capture) in the attine tree we performed one additional analysis with this taxon included. We generated a new set of alignments filtered at 75% taxon occupancy (‘Attine-119T-F75’), concatenated the loci, and then performed an unpartitioned analysis with RAxML. We generated a time tree for the evolution of fungus-farming ants using BEAST v. 1.8.2 [65]. To calibrate the analysis, we used nine fossil calibrations, and one secondary calibration (electronic supplementary material, table S6). To decrease computation time we (i) used a constraint tree and turned off tree search operators, and (ii) used subsets of UCE loci rather than the entire matrix. For five separate matrices (each with 20 randomly selected UCE loci), we performed two independent BEAST runs, each progressing for 50 million generations. We assessed burn-in, convergence among runs, and run performance by examining log files with Tracer v. 1.6. We generated chronograms for each of the five matrices separately and for all runs combined. For additional detail, see electronic supplementary material, appendix A2. We inferred the biogeographic history of the fungus-farming ants using the R package BioGeoBEARS (BGB) [66,67]. For the tree, we used the BEAST time tree pruned to include only the fungus-farming ants and their sister group. We coded taxa for the following areas: (A) Nearctic; (B) Middle America (including the Caribbean); (C) South America; (D) Afrotropics; and (E) Australasia. Using BGB, we compared six different biogeographic models: DEC [68], DEC + J, DIVALIKE [69], DIVALIKE + J, BAYAREA [70] and BAYAREALIKE + J. For each model we performed a time-stratified analysis using the time periods 0–5 Ma, 5–35 Ma and 35–65 Ma, which correspond to post-closure of the Isthmus of Panama [71], pre-closure of the Isthmus of Panama and pre-glaciation of Antarctica [72], respectively. For additional detail, see electronic supplementary material, appendix A2. We performed trait reconstruction analyses to examine the evolution of attine ant (i) agriculture and (ii) habitat preference. In both cases, we used the same pruned BEAST tree that we used for biogeography. For agriculture, we coded taxa as practising (0) no agriculture, (1) lower agriculture, (2) coral-fungus agriculture, (3) yeast agriculture, (4) higher agriculture or (5) leaf-cutter agriculture (agricultural systems reviewed in [23]). We then used the ‘ace’ function from the R package APE [73] to test three different reconstruction models: equal rates (ER), symmetrical rates (SYM) and all rates different (ARD). For habitat preference, we coded taxa as occurring in: (A) continuously wet habitat (rainforest), which should be hospitable all year long to free-living populations of attine ant fungal cultivars; (B) seasonally dry habitat (e.g. cerrado, savannah, desert, dry scrub), which would be inhospitable, at least for part of the year, to free-living cultivars; or (AB) both wet and dry habitats. We inferred ancestral habitat preference using BGB, because it more realistically treats habitat AB as a combination of habitats rather than as a distinct, third character state. We performed six analyses comparing the same models used for biogeographic inference. We investigated diversification dynamics in the fungus-farming ants using two approaches. We tested for significant shifts in diversification rates across the entire tree using the R package TreePar [74], which is an ML-based program that allows for non-constant diversification rates and incomplete taxon sampling. We also tested for rate shifts among lineages using the Bayesian program BAMM v. 2.5 [75–78] (see also [79,80]) and the accompanying R package BAMMtools [81]. A useful feature of BAMM is that it allows for non-random, incomplete taxon sampling via the input of a sampling probability file, which we incorporated into our analysis (electronic supplementary material, table S7). For additional detail, see electronic supplementary material, appendix A2. All phylogenetic analyses (ML, BI and ST) recovered the fungus-farming ants as a highly supported clade, with the dacetines (referred to here as the subtribe Dacetina and including the genera Acanthognathus, Colobostruma, Daceton, Epopostruma, Lenomyrmex, Mesostruma, Microdaceton) as the sister group (figure 1; electronic supplementary material, figures S1–S11 and table S8), confirming a previous result [42]. Support for the monophyly of Dacetina was high (more than 90%) in all analyses except for the ST analysis (6%). Within the fungus-farming ant clade, there was no topological variability among analyses and nearly all nodes received maximum support. Most relationships are congruent with previous molecular studies [23,32,82], but we recovered a number of novel results within the neoattines. Most importantly, we confidently resolved the positions of major Cyphomyrmex lineages and several enigmatic species of Mycetosoritis for the first time (both genera are non-monophyletic). Relationships among remaining genera within the more inclusive tribe Attini were less well resolved. The tribe and several groups of genera within the tribe were recovered with high support, but most relationships among genus groups were poorly supported. We did not recover the sister group to the Dacetina + Attina clade with high confidence, although most ML analyses suggested that the Blepharidatta group (Blepharidatta, Wasmannia, Allomerus) was sister to Dacetina + Attina. Figure 1. Evolution and biogeography of the fungus-farming ants and their agricultural systems. The tree topology (left and right sides) matches the best tree topology recovered from all partitioned RAxML analyses of the Attine-118T-F75 matrix (950 UCE loci, 652 774 bp). We estimated divergence dates using BEAST and 10 node calibrations (electronic supplementary material, table S6). Left side: numbered nodes received less than 95% support in at least one of six analyses (raxml-rcluster/raxml-kmeans/raxml-jacknife/exabayes-kmeans/raxml-ry-coding/astral) and the colour of the node corresponds to the frequency with which that node was recovered across 10 analyses (electronic supplementary material, table S8; black = 10/10, purple = 8–9/10, green = 5–7/10, sky blue = < 5/10). The asterisk (*) signifies 100% clade support and the dash (-) signifies that the clade was not recovered in the best tree. We mapped the five distinct attine ant agricultural systems (lower agriculture, coral-fungus agriculture, yeast agriculture, higher agriculture and leaf-cutter agriculture) onto the tree using ML-based trait reconstructions. Right side: coloured squares indicate current or ancestral geographical ranges, with the ancestral ranges inferred using the program BioGeoBEARS (DEC + J model). We used the following ranges: (A, blue) Nearctic, (B, green) Middle America, (C, yellow) South America, (D, purple) Afrotropics and (E, orange) Australasia. Coloured branches indicate current and ancestral habitat preference (blue, wet habitat; red, dry habitat; turquoise, wet and dry habitats), with ancestral preference inferred using BioGeoBEARS (BAYAREALIKE + J model). Both sides: The bars at the bottom of each chronogram provide stem- and crown-group ages for each agricultural system. Dotted lines correspond to the 95% HPD of the BEAST divergence date estimates. The wavy, light grey line depicts average global temperature (adapted from [72]). For reference, major global events are highlighted on the geological timeline (EECO, Early Eocene Climatic Optimum; TEE, Terminal Eocene Event; MMCO, Mid-Miocene Climatic Optimum). The three vertical black bars on nodes mark rate shifts identified by BAMM (all rate increases).
We recovered roughly identical divergence dates from the five different data subsets that we analysed and from the combined averaging of all results (electronic supplementary material, table S9). Additionally, our ancestral reconstructions of ant agriculture returned essentially identical results across the three different models that we tested, with the ER model being statistically favoured (electronic supplementary material, table S10 and figures S12–S14). Our date and trait inferences are largely congruent with previous studies [23,32], except that we recovered slightly older dates for the origin of higher fungus farmers. The Attini sensu lato evolved 66 Ma (56–76 Ma HPD) at the K–Pg boundary followed by what was a rapid radiation of all of the major genus groups and subtribes. The fungus-farming ants originated as lower fungus farmers at the end of the Palaeocene sometime between 61 Ma (52–70 Ma) and 57 Ma (48–66 Ma), and the two major clades within fungus-farming ants, the ‘palaeoattines’ and ‘neoattines’, evolved as lower fungus-farmers more or less simultaneously at 49 Ma (39–58 Ma) and 50 Ma (42–59 Ma), respectively. The transition from lower agriculture to pterulaceous (coral-mushroom) agriculture (Apterostigma pilosum group) occurred sometime between 21 Ma (17–26 Ma) and 16 Ma (11–20 Ma). Within the neoattines, the transition from lower agriculture to yeast agriculture (Cyphomyrmex rimosus group) occurred between 33 Ma (27–39 Ma) and 14 Ma (9–19 Ma). The evolution of higher agriculture occurred between 31 Ma (26–37 Ma) and 27 Ma (22–33 Ma) and the leaf-cutting ants (although not necessarily leaf-cutting agriculture [22,27,83]) evolved most recently during the Miocene between 19 Ma (15–24 Ma) and 18 Ma (14–22 Ma). BioGeoBEARS model comparison favoured the DEC + J model as the best-fitting model for the biogeographic data (figure 1; electronic supplementary material, figures S17, S18 and table S10), and most models inferred similar scenarios of ancestral range evolution (electronic supplementary material, figures S15–S26). The fungus-farming ants originated in South America and have maintained a strong presence in the region, as evidenced by the fact that most clades also have South American origins (59 out of 80 nodes). The ancestors of both yeast-farming attines and coral-fungus-farming attines evolved in South America. Several species and clades, however, dispersed from South America into Middle America and less often into the Nearctic. One major clade inferred to have had an origin in Middle America is the Cyphomyrmex wheeleri group, which includes one species that occurs far into North America (C. wheeleri). Notably, the other major clade to originate in Middle America is nested inside of the higher fungus farmers and includes the leaf-cutter ants. This clade is composed of the Trachymyrmex intermedius group, the T. septentrionalis group, which is an exclusively Nearctic lineage, and the leaf-cutting ant genera Atta and Acromyrmex. The leaf-cutter ants also originated within Middle America with later dispersal into both South America (multiple times) and the Nearctic (at least twice). For habitat preference, BGB model comparison favoured the BAYAREALIKE + J model, which produced results nearly identical to the next-most-favoured model, DEC + J (figure 1; electronic supplementary material, figures S27–S38 and table S10). The ancestral fungus-farming ants evolved in continuously wet rainforest habitat; however, the two major clades of fungus farmers, the palaeoattines and neoattines, show very different trends, with the palaeoattines diversifying predominately in wet habitat, and the neoattines shifting into and diversifying predominately in dry habitat. Within the neoattines, there were two radiations that occurred predominately in wet habitat, one in the Cyphomyrmex rimosus plus C. wheeleri groups and one in the Trachymyrmex intermedius group. In addition to the general trend of neoattines diversifying in dry habitat, we found that both higher agriculture and leaf-cutter agriculture originated in dry habitat, contradicting the prevailing view that most evolution and domestication in the fungus-farming ant–fungus mutualism occurred in South American rainforest. The tree-wide diversification rate analysis using TreePar found one significant rate shift across the fungus-farming tree at 6.5 Ma (p < 0.01, electronic supplementary material, table S11), with the net diversification rate showing a slight decrease. The among-lineages diversification rate analysis using BAMM found evidence for three rate increases (figure 1; electronic supplementary material, figures S39–S43), with the three-shift model receiving the highest posterior probability when compared with competing models (electronic supplementary material, table S12). The rate shifts occurred along the branches leading to Myrmicocrypta, Apterostigma (minus A. megacephala), and the neoattines minus several depauperate lineages (Cyatta, Kalathomyrmex, Mycetarotes). Counter to our expectations, none of these shifts directly corresponded with major shifts in ant agriculture; however, the neoattine shift is loosely correlated with climate change (see discussion below). Our model-based results (figure 1) support the long-held (but previously untested) hypothesis that fungus farming originated in the rainforests of South America [30,34–37,84,85]. However, in stark contrast to the consensus view that most subsequent evolution also occurred in South American rainforests (although see [86]), we found that the ancestors of most major attine lineages, including the ancestors of higher fungus farmers and leaf-cutter ants, probably evolved in dry or seasonally dry habitat (cerrado, savannah, desert etc.). This result is significant because it suggests a possible mechanism for the origin of the symbiotically obligate, highly derived, higher attine fungi. As in human-mediated domestication, the domestication of fungi by ants almost certainly required, at least temporarily, the genetic isolation of fungal cultivars from their free-living, wild-type conspecifics. Based on the results of our habitat-preference analyses, we hypothesize that the dispersal of fungus-farming ants into dry or seasonally dry habitats may have been the key factor driving fungal isolation and subsequent coevolution. Reduced fungal outcrossing by ant dispersal into dry habitat is supported by the observations that leucocoprineaceous fungi, the group that includes most attine ant cultivar species, are most abundant and diverse in rainforests [40,41], and that all free-living conspecifics of attine cultivars have been collected only in the rainforests of Brazil and Panama [20,24]. Unable to escape the humid confines of their host nests [87,88], selection would eventually have driven higher attine fungi to become obligate ant nest specialists. In addition to shifts into dry habitat, our results suggest fungus-farming ants expanded their range into Middle America and the Nearctic (figure 1) multiple times, with the earliest dispersal out of South America occurring around 27–22 Ma—well before the traditional estimate for the closure of the Isthmus of Panama (approx. 3 Ma) [71] (but see also [89–91]). Intriguingly, the sister group to the leaf-cutter ants is a lineage of completely Nearctic species in the Trachymyrmex septentrionalis group and the ancestor of the leaf-cutting ants probably originated in Middle America. If a shift into dry, inhospitable habitats provided the conditions for strict ant-fungus coevolution, it might also be the case that dispersal out of South America and into peripheral areas played a part as well. Movement into Middle America, and especially into the Nearctic, would have further isolated fungal cultivars from their parent populations and free-living progenitors. In the case of higher fungus farmers, geographical shifts would have the effect of isolating the ants from access to alternative species of higher attine fungi under cultivation by other higher attine species in the ancestral range. Our results present a more detailed picture of the overall biogeographic history of fungus-farming ants. The tribe Attini, to which the fungus farmers belong, originated at about the same time as the K–Pg mass extinction event. Our phylogeny suggests that there was probably a burst of diversification at this time as evidenced by the short internodes subtending major clades and the fact that we could not confidently resolve relationships among most genus groups, even with genome-scale data. The fungus-farming ants originated shortly after the K–Pg extinction event between 61 Ma and 57 Ma, possibly during the post-extinction-event recovery period [92] and shortly before the early Eocene climatic optimum. It is during this time that the ants began their symbiosis with leucocoprineaceous fungi and lost their ability to synthesize arginine [27,33], committing them to agrarian life. Transitions from lower to other agricultural systems (five systems total; see figure 1) occurred after the terminal Eocene event (TEE) at approximately 35 Ma. This event involved a major drop in global temperature, which brought glaciers to Antarctica and began the expansion of drier habitats throughout the New World [93,94]. Although dry habitats such as grasslands have been present since the early Eocene [94], several studies have noted significant expansions of C4 grasses starting around 30 Ma [95–97], remarkably close in time to the origin of higher agriculture. It is also notable that we inferred a diversification rate increase in the lineage that includes all neoattines, minus a grade of several depauperate groups (Cyatta, Kalathomyrmex, Mycetarotes and Mycetosoritis hartmanni), at about the same time as the TEE. Given that most genus-level lineages and all derived agricultural systems originated after the TEE, it is plausible that this global cooling event spurred both attine ant lineage diversification and ant–fungus coevolution. A critical, long-standing problem for understanding the evolution of ant agriculture has been identifying the non-fungus-growing sister lineage of the fungus-farming ants. The identity of the sister group is important because it could provide critical information about the behavioural, physiological or ecological precursors to fungus farming. Twelve different lineages have been variously proposed and, of these, morphological studies have supported the cryptic leaf-litter ant Blepharidatta brasiliensis as the sister group [98–100]. More recent molecular studies, however, have provided low to moderate support for a sister-group relationship between the fungus-farming ants and the subtribe Dacetina [23,32,42]. Using phylogenomic data and fundamentally different analytical paradigms, we found nearly unequivocal support for a sister-group relationship between the fungus-farming ants and the subtribe Dacetina. The dacetines form a group of largely specialized predatory ants that occur in the Neotropical (Acanthognathus, Daceton, Lenomyrmex), Afrotropical (Microdaceton) and Australasian (Colobostruma, Epopostrma and Mesostruma) regions. Most of the genera have elongate, trap-jaw mandibles. Biological information for the dacetines is surprisingly limited, with the best-studied genera being Acanthognathus and Daceton [101–108]. Importantly, these ants have no known associations with fungi. The Dacetina and Attina do, however, share a few morphological and behavioural traits [104,109]. Most notably, both sister clades are highly specialized, the Dacetina specialized predators and the Attina specialized agriculturalists. Although these are behavioural extremes, it is possible they share a common historical origin. Given that both of these lineages evolved from generalized hunter–gatherers during the ‘nuclear winter’ that followed the K–Pg extinction event (figure 1) [92,110], it is possible that a reduction of generalized resources drove them to specialize, one on live prey and the other on fungi. As noted by Janzen [110], the best survivors of the K–Pg nuclear winter were probably those whose food did not directly depend on immediate photosynthesis. Specialization on fungi by the Attina may have been driven both by a reduction in generalized prey as well as by a proliferation of fungi [111]. The importance of animal and plant domestication to the rise of modern human civilization is well understood. Through a variety of mechanisms, including selective breeding and genetic isolation, humans increased food yield and nutritional value, fuelling the growth of civilizations. Similarly, in fungus-farming ants, the evolution of agriculture and subsequently of high-yield domesticated crops resulted in the rise of the ecologically dominant leaf-cutting ants, which form ‘superorganism’ colonies composed of millions of individuals. Unlike human agriculture, however, an explanation for how ants unconsciously domesticated their fungal cultivars has been highly uncertain. Our results provide the first evidence that fungal domestication occurred in dry habitats, suggesting an explanation for how fungal cultivars may have evolved into domesticated mutualists obligately dependent on their ant hosts. For all newly sequenced samples, raw sequence reads and Trinity contig assemblies representing UCE loci are available from the NCBI Sequence Read Archive and GenBank, respectively, under BioProject accession PRJNA379607. Additional data, including alignments, alignment supermatrices, tree files, partitioning files, BEAST xml files, trait and biogeography data files, and tables are available from Dryad (http://dx.doi.org/10.5061/dryad.d7pr6) [112]. M.G.B., T.R.S. and S.G.B. conceived of the study. M.G.B. carried out the molecular lab work and performed the analyses. M.G.B. and T.R.S. drafted the manuscript. A.J., J.S.-C. and M.W.L. helped carry out the molecular lab work and revise the manuscript. B.C.F. and S.G.B. helped with data analysis and manuscript revision. All authors gave final approval for publication. We declare that we have no competing interests. M.G.B. was funded by a Peter Buck postdoctoral fellowship at the Smithsonian Institution and NSF grant no. DEB-1354996 (Project ADMAC). A.J. and J.S.-C. were funded by Peter Buck predoctoral fellowships at the Smithsonian Institution. T.R.S., J.S.-C. and A.J. were supported by NSF grant nos. DEB-0949689 and DEB-1456964. S.G.B. was supported by NSF grant no. DEB-1555905. T.R.S. was supported by the National Museum of Natural History (NMNH) Small Grants Program and the Smithsonian Institution Scholarly Studies Program. We are grateful to the following individuals for provision of ant specimens: Alan Anderson, Carlos Roberto Brandão, Jacques Delabie, David Donoso, Rodrigo Feitosa, Brian Fisher, John Longino, Caué Lopes, Jonas Maravalhas, Claudia Medina, Ulrich Mueller, Eugenia Okonski, Christian Rabeling, Heraldo Vasconcelos and Philip Ward. We thank Christian Rabeling for providing us with unpublished estimates of Mycocepurus species diversity. We thank Eugenia Okonski for laboratory assistance. For sequencing assistance we thank Joe DeYoung at UCLA and Peter Schweitzer at Cornell. Lab work was conducted at the Smithsonian NMNH Laboratory of Analytical Biology (LAB), and phylogenetic analyses were performed using computing resources at the Smithsonian Institution and the CIPRES Science Gateway [113]. We thank two anonymous reviewers for helpful suggestions about the manuscript. FootnotesElectronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3726862. © 2017 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited. References
Page 4Optical variation among aquatic ecosystems has strong effects on fundamental ecological processes [1]. In shallow waters, a major driver of optical variation is the concentration of dissolved organic carbon (DOC) [2,3]. Essentially, by affecting the light environment, DOC drives a variety of ecological processes, from consumer foraging and habitat use [4], to the thermal stratification of lakes [5]. But perhaps the greatest effect of DOC is its effect on primary production. By strongly absorbing photosynthetically active radiation, DOC inhibits primary production at high concentrations [6–8] indirectly reducing zoobenthic production [9], the abundance of predatory fishes [6,10] and their somatic growth [11]. Yet, the effects of optical variation due to DOC on evolutionary processes are less well understood, despite the fact that many aquatic organisms employ visual signals during courtship and territorial displays. Two types of mechanisms, evolutionary and ecological, appear likely to drive DOC-mediated effects on visual signals and sexual communication. First, light limitation can presumably reduce the effectiveness of visual signals, ultimately shaping their evolutionary trajectory through natural selection on signal performance (e.g. transmission). These selection-based explanations for signal–environment correlations typically draw upon evolutionary models such as sensory drive [12,13] to provide a mechanistic basis for divergence. Indeed, empirical evidence lends support to these evolutionary mechanisms by showing that spatially heterogeneous optical conditions can drive divergence in sensory systems and signal attributes [14–16]. Second, light conditions can also limit signal production and maintenance, representing a more proximate, ecological, basis for signal divergence. We believe this inhibitory effect can arise via two nutritional pathways, pigment limitation and energetic limitation. The development of many colourful signals in animals depends on the acquisition of plant-derived pigments (i.e. carotenoids) [17]. As discussed previously, DOC is likely to limit the availability of these algal-derived pigments by stifling primary production at high concentrations. Whether bottom-up limitation of sexual signals could be driven by DOC is unknown. However, existing evidence for canopy-driven limitations on carotenoid supply and signal development in guppies suggests it may play an important role [18]. Alternatively, because signals are energetically costly to produce and maintain, light limitation may simply have negative effects on signal production by driving reduced male condition. Since most elaborate signals appear to have some degree of condition dependence, such a link between low mean condition and low mean signal elaboration seems likely. Our primary objective in this study was to better understand the effects of optical variation on interpopulation divergence in sexually selected visual traits. We explored this effect by sampling a series of 21 Bahamian mosquitofish populations on Abaco Island (Gambusia bahamasminimus) across an optical gradient driven primarily by DOC. We focused on three male traits used during courtship displays and subject to female preference: dorsal fin redness, anal fin redness, and gonopodia length [19,20]. While DOC was the focal driver in our study, other ecological costs such as predation are potent evolutionary agents (reviewed in [21,22]) that also drive the adaptive evolution of sexual signals in this system [20,23,24], and many others [25–28]. Therefore, we included predation surveys in our study to assess this effect in concert with the spectral environment, and in doing so we offer an assessment of their individual and combined effects on sexual signals. This study was conducted on Abaco Island, The Bahamas. Mosquitofish populations on Abaco are of an undescribed species, tentatively named Gambusia bahamasminimus (previously G. hubbsi), within a Bahamian clade consisting of at least three phenotypically and ecologically similar species [29]. Bahamian mosquitofish are small, live-bearing, fish found in a variety of aquatic habitats [24,30,31]. Males exhibit enlarged, orange-red coloured, condition-dependent dorsal fins subject to sexual selection via female preference [20,23,24]. Males also possess anal fins with orange-coloured membranes and an elongate sperm-transfer organ formed from modified fin rays 3–5, the gonopodium [23]. Fin coloration in Bahamian mosquitofish has heritable components, but also seems affected by diet [23]. For this study we sampled Gambusia populations between June and August in 2011. To encompass a wide range of spectral conditions we sampled a variety of aquatic ecosystems throughout Abaco including blue holes, freshwater marshes and mangrove-lined tidal creeks. From each population we collected 4–71 adult males, measured their standard length (SL), gonopodium length (GL) and photographed them live for colour analysis following Giery & Layman [24]. Fin coloration was measured in RGB colour space using Adobe Photoshop CS5. RGB defines colour in three-dimensions where R indicates red, G indicates green and B indicates blue. Following Endler [32] we use RG, an index of signal coloration (R − G)/(R + G) that estimates colour (hue) along an axis from red (1) to yellow (0) and green (−1) which encompasses the range of colour (yellow-red) expressed in Bahamian Gambusia (figure 1; electronic supplementary material). Figure 1. Male Bahamian Gambusia vary widely in coloration among populations. Individuals in panel (a) represent the diversity of coloration exhibited among populations on Abaco Island, The Bahamas. Panel (b) depicts the diversity of optical environments characterizing mosquitofish habitats. (Online version in colour.)
To describe the optical conditions in each site we passed a beam of light between two fibre optic cables fitted with co-limiting lenses positioned 10 cm apart and submerged in a basin of water collected from each site. Water colour was calculated as the amount of light at particular wavelengths reaching the spectrometer (Jaz, Ocean Optics, FL) after being transmitted through a sample, proportional to a distilled water standard (e.g. T = water sample/control). We restricted transmission measures to 370 and 570 nm because these wavelengths correspond with peak sensitivities for short and long-wavelength cones in poeciliid fishes (e.g. [33,34]), and piscivorous fishes ([35], electronic supplementary material). Gross patterns in water colour, i.e. the spectral shift that typifies our study system (reddening), were estimated with an index, TRB. This index is the difference in transmission between long (red) and short (UV/violet) wavelengths [(T570 − T370)/(T570+T370)] with larger positive TRB values indicating long wavelength dominated environments and large negative values indicating environments dominated by short wavelengths (hereafter called ‘red’ and ‘blue’ environments, respectively). We sampled DOC concentrations between June and November 2012. Water samples were passed through precombusted (400°C for 4 h) Whatman GF/F filters. Filtered samples were analysed in situ with a Turner Aquaflor fluorometer and recorded in relative fluorescence units (RFU). RFUs were converted to mg C l−1 using field-collected samples measured for TOC concentrations following method 5310B [36]. Two to four samples were collected from each study site across the sample period. Repeatability was high among samples within sites (r = 0.8) indicating consistent DOC concentrations across the sampling period [37] so DOC measures were averaged for analysis (electronic supplementary material, table S1). Annual variation in DOC was not assessed in this study. However, previous study on Abaco shows that water colour, a product of DOC concentration in part [2], is largely consistent from year to year (see table S1 in [30]). To assess the contribution of DOC to site-level variation in optical environments we regressed TRB against DOC (log-transformed). We also generated qualitative estimates of predator risk (predatory fish present/absent) for each site using visual surveys following Giery & Layman [24]. We began our analysis by checking for heterogeneous trait allometries between predation regimes [19,29]. We used linear mixed models (LMM) in which focal signalling traits were dependent variables. SL(log-transformed), predation, and their interaction were fixed effects. Population and population*SL were random effects in order to allow population-level heterogeneity in slopes and intercepts. There was no significant effect of predation on the slopes of the SL – GL, SL – RGd, or SL – RGg relationships (p = 0.4; p = 0.56; p = 0.72, respectively). Predator-dependent allometry was not considered further. To test for effects of water colour and predation on signal morphology we employed linear mixed models. We ran separate models for each focal trait: GL, RGd and RGg. Initial visual inspection of our raw data (e.g. RGd on TRB) suggested a quadratic effect of TRB for fin coloration. While not part of our original prediction, we chose to include quadratic TRB terms in all models. Independent variables included log-transformed standard length (SL), water colour (TRB), predator regime, the interaction between TRB and predator regime, a quadratic term ( We also assessed the effect of geographical distances on the similarities of site conditions (TRB, predator) and size-corrected morphological traits (RGd, RGg and GL) using Mantel tests. p-values were generated from 9999 simulations. DOC concentrations were highly variable among sites (electronic supplementary material, table S1), but typical for coastal ecosystems in the region [39]. Optical environments (TRB) were also extremely variable among sites (electronic supplementary material, table S1) and strongly associated with DOC (F1,19 = 129.9, p < 0.001, R2 = 0.87; figure 2). Piscivores such as snappers (Lutjanus spp.), needlefish (Stongylura spp.) and great barracuda (Sphyraena barracuda) were detected coexisting with 14 mosquitofish populations (electronic supplementary material, table S1). The predator status of each population has been maintained for at least several years, see Giery & Layman [24]. Importantly, a difference between predator and no-predator habitats in optical conditions TRB was not detected (t = −1.4, d.f. = 19, p = 0.20). Figure 2. Dissolved organic carbon (DOC) strongly absorbs short and medium-wavelength light generating an asymptotic relationship between DOC and water redness (note the log-transformed x-axis). Non-predator populations are displayed in black predator populations in red (or grey). (Online version in colour.)
We captured, photographed and subsequently measured 879 male mosquitofish (electronic supplementary material, table S2). Visual signals expressed by male mosquitofish varied considerably among populations (electronic supplementary material, table S2). Our reduced LMM explaining interpopulation variation in fin colour (RGd and RGg) included linear and quadratic water colour terms (table 1). Yet, the only significant effect in each model, besides the body size covariate (SL) was the quadratic term (p < 0.03) indicating a unimodal relationship between water colour and the coloration of dorsal and anal fins (figure 3). For gonopodium length, all focal terms besides SL were removed during model reduction indicating little-to-no effect of observed ecological variation on gonopodium length (figure 3 and table 1). Figure 3. Water redness (TRB) appears to have a strong nonlinear effect on mean (+ standard error) sexual signal coloration for 21 populations of Bahamian Gambusia: mean dorsal fin redness (RGd), top panel; and mean anal fin (RGg) redness, middle panel. No effect of water colour was detected for gonopodium length (GL), bottom panel. Predation was dropped from all models during model selection, however populations subject to predation are displayed here in red (or grey), non-predator populations in black. All traits are partialed for the effect of body size (Log SL). (Online version in colour.)
Geographical distance was not correlated with morphological traits (r = 0.06, p = 0.26) or ecological characters (r = 0.01 p = 0.46) indicating that phenotypic and ecological variation among sites are not driven by geographical proximity alone. Simple mantel tests performed for each morphological and ecological variable also failed to reveal a significant geographical effect (electronic supplementary material, table S3). High concentrations of DOC removed a substantial fraction of short-wavelength radiation from aquatic habitats and generated a strong gradient in water colour (figures 1 and 2). In turn, spatial heterogeneity in the optical environment appeared to drive phenotypic divergence in sexual signals used by Bahamian mosquitofish. We found no evidence for a predation effect in this study suggesting a relatively weak influence of predation pressure on sexual signals in our system—a result consistent with some of our previous findings [30], but inconsistent with others [23,24]. Despite this contribution, the most notable result from our study was not that optical environments are more important for sexual signal expression, but rather that this effect was non-monotonic (i.e. ‘hump-shaped’). Most interpretations of light-mediated signal divergence in fishes favour evolutionary processes such as sensory drive to explain interpopulation variation in signal coloration (e.g. [40–42]). Although the direction of divergence appears to vary among species (negative or positive slopes along optical gradients), linear or near-linear relationships appear typical for continuous traits and are commonly regarded as indicators of an adaptive response to optical conditions [40,43]. However, the form of the response we found differs substantially from these existing empirical data and their underlying theoretical prediction. We believe that this difference is due, at least in part, to bottom-up ecological controls on signal elaboration, but we can only speculate without further examination of underlying mechanisms. Following, we discuss several hypotheses about how ecological variation might drive the observed pattern of phenotypic differentiation in order to stimulate new research in this field. We begin with evolutionary mechanisms and follow with ecological ones. Our first two hypotheses invoke a role for functional thresholds in the evolution of visual communication. If phenotypic (signal) variation is constrained, selection on signal efficacy may produce thresholds at which selection for maximum conspicuity (e.g. via chromatic contrast) favours complex patterns of phenotypic divergence. For example, if signal coloration is constrained between yellow and red, then across a gradient of spectral conditions (e.g. blue to red), optimal signal design via chromatic contrast requires a nonlinear relationship between signal and environment to maximize conspicuity: e.g. yellow signals in blue environments, red signals in green environments, and yellow signals in red environments. The role of constraints in the adaptive diversification of visual signals has not been a focus of study, yet such an effect seems compatible with sensory drive hypotheses [13,44], with empirical studies lending support to this. For example, Fuller [45] shows that colour-polymorphic bluefin killifish (Lucania goodei) maximize conspicuity rather than maintain specific hues across a DOC-driven water colour gradient. The result is that blue-finned fish are common in red environments, and red-finned fish are common in blue environments. Such a pattern suggests that sexual selection may act on signal conspicuousness via chromatic contrast when the optical environment is highly variable. While such a strategy seems well-established, it has never been shown for fishes that are not colour polymorphic, such as Bahamian mosquitofish. Therefore, a role for evolutionary constraints in these cases remains unexplored. It is also plausible that the yellow coloration exhibited by fishes inhabiting the most red optical conditions simply reflects lowered investment in signal production. Disruptive transmission environments can relax selection on sexual signal production [46,47]. If the transmission environment is antagonistic to effective signalling, investment in signal production may constitute a substantial fitness cost [48] which could explain the predominance of yellow signals (i.e. less carotenoid rich) that we observe in disruptive signal environs. While we can only speculate as to whether divergent signal coloration across the DOC gradient is adaptive, the idea of colour signals ‘tuned’ to suit environmental variation generates appealing hypotheses with abundant support from a variety of study systems. Alternatively, the effect of DOC on sexual signal coloration may be entirely ecological. That is, DOC-mediated light limitation may simply suppress sexual signal production via bottom-up controls. Light-limitation in oligotrophic aquatic ecosystems (like those in The Bahamas) strongly regulates productivity [6,9]. These effects can propagate up trophic levels, affecting zooplankton abundance, and reducing individual growth rates and population sizes of vertebrate consumers (i.e. fish; [6,10,11]). Therefore, at high DOC concentrations, a bottom-up reduction in primary productivity and pelagic zooplankton biomass could in turn limit the availability of carotenoids and/or the energy available for the production of red sexual signals [17,18]. Interestingly, several researchers also have shown unimodal responses in fish population density [10] and growth rate [11] across a gradient of DOC. They attributed this response to an ecological threshold at which DOC switches from a stimulator of lake productivity to an inhibitor via shading effects. Similarities in the response to increasing DOC, although superficial, point to a nonlinear, mechanistic link between light limitation, ecosystem productivity, growth, and the production of sexual signals in our system. Indeed, when we explored this link further using male Gambusia growth data collected previously from nine of our study sites [31], we found that growth rates also followed a unimodal pattern along the DOC gradient (electronic supplementary material, figure S1) providing corroborative evidence for bottom-up regulation of fish secondary production at high DOC concentrations. While merely suggestive, these supplemental data provide additional evidence for strong, bottom-up, ecological control of signal production via nutritional pathways. And while evolutionary hypotheses such as sensory drive provide attractive explanations for observed correlations between signal and optical environment, we suggest that one should not overlook the ecological effects that accompany variable optical conditions and the potential for bottom-up regulation of sexual signals. What do these results, and those of previous studies on Bahamian mosquitofish, suggest about the importance of predation for the evolution of sexual signals in this system? Here, we find no significant effect of predators on interpopulation variation in sexual signals, a result shared with a previous archipelago-wide study [30], yet inconsistent with others [23,24]. These inconsistencies highlight heterogeneous prey responses to variable predation pressures in these systems, a tendency seen in guppy sexual signals as well [49,50]. Underlying causes of variable responses to predation remain unclear, however several explanations seem likely. First, antipredator adaptations may compensate for increased cost of signalling under high predation risk [22,51,52]. Second, predators might also regulate bottom-up drivers of trait variation [18] through trophic cascades — a hypothesis in line with the light limitation hypothesis detailed above. Essentially, these inconsistencies show that much about the effects of predators on sexual selected traits remains to be explained. Nevertheless, it seems clear that detectible effects of predators, in this system and others, are variable or at least hard to predict based on simplified scenarios of ecological costs [25,53]. At a more fundamental level, our results provide evidence that ecological gradients can generate complex patterns of phenotypic variation. In turn, nonlinear phenotypic responses to these gradients may have important evolutionary and ecological implications [54]. For example, non-monotonic divergence along environmental gradients could dampen local adaptation if phenotypes at disparate ends of ecological gradients are convergent, as in our system. Therefore, assortative mating could weaken local adaptation by facilitating gene flow among populations [55]. This suggests that adaptive divergence in sexual signals along environmental gradients may not necessarily favour reproductive isolation due to assortative mating by signal phenotype, a common phenomenon in fishes [56], and a hypothetically important mechanism in models of speciation with gene flow (e.g. [57,58]). Despite a wealth of study on a suppressive role for predation regarding the elaboration of sexual traits, we found no such effect in our study. Rather, our data indicate that water colour was a strong driver of variation in the colour of a sexually selected trait. While these findings are important, we note that the more intriguing result is the complex relationship between sexual signal coloration in Bahamian mosquitofish and the optical properties of the environmental. In fact, this study is one of the first to identify non-monotic divergence in sexual signals due to environmental variation. Further study of the underlying mechanism(s) and evolutionary implications of this divergence pattern is needed. Nevertheless, our results clearly illustrate that the coupling of evolutionary and ecological dynamics are important for sexual signalling systems, and given the prevalence of optical heterogeneity in nature, other systems are also likely to exhibit complex divergence patterns. All work was approved by the Animal Care and Use Committee of the University of North Carolina (protocol no. 14-057-A) and The Bahamas Department of Fisheries. Data are available from the Dryad Digital Repository at: http://dx.doi:10.5061/dryad.gh225 [59] S.T.G. conceived of the study, collected and analysed data. S.T.G. and C.A.L. wrote the manuscript. The authors have no competing interests. This research was supported by NSF (OCE #0746164; C.A.L., DEB #1406399; S.T.G.), a FIU SEEDS grant (C.A.L. and S.T.G.), and a FIU DEA fellowship (S.T.G.). We thank Friends of the Environment, D. Hayes, K. Rennirt, and James Richard for field support. B. Langerhans, J. Gilliam, J. Stroud, N. Lemoine, M. Fork, M. Zimova, J. Abbott, and two anonymous reviewers provided constructive comments on this manuscript. FootnotesElectronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3719854. References
Page 5Angiosperm evolution has given rise to an overwhelming diversity of floral morphologies adapted to pollination by a multitude of different vectors. This diversity is mirrored in the high variability of breeding systems and reproductive strategies across angiosperms. Hence, it is hypothesized that floral form and function have important effects on diversification [1–4]. There is an extensive body of literature on floral morphology, pertaining both to extant and extinct taxa [5–11]. However, the distribution of flower morphological diversity across major subclades, let alone across the angiosperms as a whole, has rarely been addressed using an explicitly analytical and synthetic approach [12,13]. Such broad-scale analyses of disparity (morphological diversity) have so far been largely restricted to animal groups [14–17]. Morphospace analyses are used to study macro-evolutionary patterns and trends in disparity within and among clades. While disparity analyses are traditionally conducted on large numbers of traits capturing the overall morphology of particular organisms [18–21], some studies have focused on sets of traits of specific functional, developmental or evolutionary significance (e.g. the morphology of animal jaws in relation to feeding behaviour [22–24]). This latter approach also allows us to account for the fact that different traits might evolve with different modes and rates [25] and at different times in the history of clades [26]. Under the assumption that traits/subsets of traits involved in different functions are subject to different evolutionary constraints and selective regimes, one might expect these traits to show different levels of disparity. Studies on significant subsets of organs are thus necessary to clearly characterize the different drivers and causes underlying the disparity exhibited by a clade; such studies provide an alternative and complementary approach to traditional comparative structural analyses. Most flowers are composed of three main functional modules. From the periphery to the centre, a flower usually comprises one or two sets of sterile organs (perianth), a set of male reproductive organs (androecium) and a set of female reproductive organs (gynoecium; electronic supplementary material, figure S1). In the perianth, sepals commonly protect younger organs during pre-anthetic stages, while petals mainly attract and guide pollinators at anthesis [27]. The function of the androecium is pollen production and presentation, and, more rarely, pollinator attraction. Finally, the main functions of the gynoecium are ovule production, pollen reception and sustenance of pollen tube growth, as well as seed protection and dissemination. The organization and development of these three functional modules are not only the basis of traditional, taxonomic descriptions and comparative analyses of floral structure, but are also the target of modern, molecular developmental (evo-devo) models of floral evolution and development such as the ABCE-model [28,29]. Recent research has also focused on the synorganization of functional units in flowers [11,30,31]. However, the allocation of morphological variation among these three floral components has never been quantified. Here, we tested whether disparity in the flower as a whole is equally reflected in the three functional modules, or whether, by contrast, one part of the flower varies more than the others. At the scale of an entire plant order, where organs can be compared at the organizational and functional level, we expected the perianth to show low variability, owing to the simple structure of petals and sepals. On the other hand, the gynoecium is a very complex structure, achieving numerous functions throughout the flower's life [7,32], and we expected it to show high variability when compared with the rest of the flower. We addressed these issues in the Ericales, a speciose order of angiosperms nested in the asterid clade of core eudicots. Ericales diverged from their sister group in the Early Cretaceous, ca 112 million years ago [33] and encompass 22 families (APG IV [34]; figure 1a), 346 genera, and approximately 11 550 species [41] displaying considerable ecological diversity [42]. In many tropical rainforests, ericalean taxa account for up to 10% of the total tree species diversity [43]. The order includes species of great economic importance, such as tea (Theaceae), kiwi (Actinidiaceae), persimmon and ebony (Ebenaceae), Brazil nuts (Lecythidaceae), sapote (Sapotaceae, Ebenaceae) and a variety of ornamental species such as heather and rhododendrons (Ericaceae), and primroses (Primulaceae). Ericales have a worldwide distribution and a considerable diversity in habit, general morphology, method of nutrient uptake, and in particular, floral morphology [42,44]. This diversity is also reflected by the fact that the identification of non-molecular synapomorphies for the order as a whole has been proven difficult, while detailed comparative studies of floral structure have identified series of potential synapomorphies for various suprafamilial clades [45–47]. Importantly, Ericales also have a comparatively rich fossil record with a series of charcoalified flower fossils from the Late Cretaceous, the geologic period during which the angiosperms began to dominate most terrestrial ecosystems [36,38,40]. As the charcoalification process preserves the three-dimensional shape of floral fossils and only leads to moderate alterations at the morphological level (e.g. shrinkage [48]), most of these fossil flowers are extremely well preserved and can be compared directly with their extant relatives [8,48]. Figure 1. The floral morphospace of Ericales. (a) Phylogenetic relationships among ericalean families; dated tree modified from [33], keeping only nodes (with crown ages given in Ma) that are supported in [35]. Pictures of six of the fossil genera included in the analyses: a Raritaniflora, b Paleoenkianthus, c Glandulocalyx, d Parasaurauia, e Paradinandra and f Pentapetalum. Assigned positions (according to original papers) of the fossils are highlighted in (a) by superscript letters on the family names. g Proposed positions of Actinocalyx (picture not shown). (b) Disparity. In blue: mean pairwise dissimilarity; in orange: maximum pairwise dissimilarity (range) rarefied to 10; in black: number of species according to [34]. Error bars are bootstrapped s.e. (c) Morphospace representation using principal coordinate analysis. Each graph corresponds to the two-dimensional representation of the space. Black dots: species of highlighted major suprafamilial clades or families; grey dots: all remaining ericalean species. (d) Illustration of floral diversity in Ericales: from top to bottom: Satyria sp.* (Ericaceae), Sarracenia flava (Sarraceniaceae), Symplocos pendula (Symplocaceae), Schima superba (Theaceae), Anneslea fragrans** (Pentaphylacaceae), Primula officinalis* (Primulaceae), Cantua quercifolia (Polemoniaceae), Couroupita guianensis* (Lecythidaceae), Impatiens paucidentata (Balsaminaceae), Mitrastemon matudae*** (Mitrastemonaceae). (e) Position of the nine fossil species (black dots) in the morphospace (grey dots). Fossil pictures: a republished with permission of The University of Chicago Press from [36]; b, e and f republished with permission of the Botanical Society of America from [37–39]; c republished with permission of Oxford University Press from [40], permission conveyed through Copyright Clearance Center, Inc.; d by P. Herendeen. Photos by *A. Weissenhofer; **T. Rodd; ***D. Breedlove, included with the authorization of D. L. Nickrent. (Online version in colour.)
For this study, we compiled an extensive dataset of extant and fossil ericalean species and built a floral morphospace based on 37 traits capturing the morphology of their flowers. Our first goal was to quantify patterns of morphospace occupation within and among ericalean families and suprafamilial clades. We tested the hypothesis (i) that various suprafamilial clades do not overlap or only overlap partly in the floral morphospace. This hypothesis derives from the fact that several clades such as, for instance, the ericoid and the primuloid clade, were not considered to be closely related in pre-molecular, largely morphology-based classifications (e.g. [49]), suggesting divergent floral morphologies. At the same time, we hypothesized (ii) that families that are supported as closely related based on comparative floral structure (such as the balsaminoid families or the polemonioid families), will occupy overlapping areas in the morphospace due to their relatively recent common ancestry. We then used our dataset to investigate (iii) whether floral disparity is coupled with clade age and/or species richness. Furthermore, we placed several Cretaceous ericalean fossils in the floral morphospace of extant taxa. Based on their old age and the fact that most of these fossils have been referred to different ericalean lineages, we hypothesized (iv) that they fall in different areas of the total floral morphospace of Ericales. Finally, we compared the disparity of the sterile, male and female parts of the flower, to test the hypothesis (v) that levels of disparity differ according to the biological or ecological function of organ modules, reflecting different evolutionary constraints and selective regimes. We sampled 381 species belonging to 275 genera (accepted in [50]). For each family, we sampled at least one species per genus, to ensure that our sample is representative of the families' floral morphological diversity. When there were less than 10 genera in a family (that was the case for 14 families, see electronic supplementary material, figure S2), we sampled at least 10 species whenever possible (electronic supplementary material, figure S2). For the families Ericaceae (126 genera), Primulaceae (68 genera) and Sapotaceae (60 genera), we sampled at least 50 genera, taking care to cover all major clades identified in phylogenetic/taxonomic studies (e.g. [51–53]). We used original species descriptions or, when available, recent taxonomic revisions, online floras, and other scientific literature, as well as personal observations from living and alcohol collections of the Botanical Garden of the University of Vienna. We scored 37 floral characters describing: general features (e.g. flower size; five characters), the perianth (12 characters), the androecium (13 characters) and the gynoecium (six characters) of the anthetic flower. These characters were chosen for their capacity to characterize the number and position of organs, organ union (organs of the same type), organ fusion (organs of different types) and organ form (for staminodes). Using taxonomic keys as our primary source, we selected all the characters that described the flower and that were applicable throughout the whole order of Ericales. For species producing unisexual flowers, characters of androecium and gynoecium organs were only scored for functional organs and not for sterile organs, such as staminodes and pistillodes. Some characters, such as organ number and size, were frequently described as polymorphic in the literature. As the frequencies of these variations were rarely detailed, we choose to code the most common state, when it was documented. For instance, ‘(4-) 5 (-6) petals’ in a description was coded ‘5’ in our dataset. The remaining polymorphic characters (e.g. ‘4–6 petals’) were coded as such, which represents 274 data entries (2.2% of all data entries; electronic supplementary material, table S2). As most of our analyses do not support polymorphisms, we randomly sampled a matrix (for each analysis), in which each of the polymorphic cells was replaced by a value comprised of the cell range (for numerical discrete data) or by one of the possible states (for binary and categorical data). Given the low amount of polymorphic data entries, this matrix sampling had no effect on the statistics calculated nor on data visualization (data not shown). Data were entered and are stored in the online database PROTEUS [54]. Each data entry, user name, source and putative notes are available in electronic supplementary material, table S2. The detailed description of the characters and character states is given in the electronic supplementary material. All our analyses were performed using the software R v. 3.0.0 [55]. Scripts are available upon request from M. Chartier and S. Gerber. We calculated pairwise dissimilarities between taxa using the mean character difference (here noted D) [56]. Let us have two taxa A and B described by N morphological characters. For a character i, the difference dABi between A and B was calculated in different ways depending on the type of character:
The mean character difference DAB between taxa A and B was finally computed as To illustrate morphological differences between the 22 Ericales families, we visualized the morphospace of Ericales with a principal coordinate analysis (PCoA; [57]) taking as input the original dissimilarity matrix. Calculation and analyses of the morphological diversity (disparity) were carried out from the original dissimilarity matrix, and not from the ordination scores. Disparity within each family was calculated as the mean pairwise dissimilarity, i.e. the mean D per family (here noted Partial disparity (PDiv, the additive contribution of each family to the disparity of the whole order) was calculated following [61]. It is the sum, over each PCoA axis, of the squared Euclidean distances between all taxa from a clade and the centroid of the whole dataset (divided by the total number of species in the dataset). Two groups falling in different parts of the morphospace are morphologically different, whereas they are similar if their distributions in the space overlap. We assessed morphological differences among ericalean families with non-parametric analyses of variance (npMANOVA, sometimes also referred to as PERMANOVA) using the function adonis() from the vegan package in R [62]. We used the original dissimilarity matrix as input, and 10 000 permutations to calculate the distribution of a pseudo F ratio under the null hypothesis. We used the same analysis as post hoc, with a Bonferroni correction. Spearman correlation tests were performed to investigate the links between disparity ( To compare variation among the 37 morphological traits, we averaged dAB, i.e. the differences between taxa, for each character. The resulting values, here noted Dchar, increase with the variation of a character in the dataset. In our dataset, the perianth morphospace is based on 12 characters, the androecium space on 13 characters, and the gynoecium space on seven characters. Because these three morphospaces differ in character composition and in size, their respective disparities cannot satisfactorily be compared directly. We first investigated if the disparity for each of the functional modules increased with total disparity in Ericales. To do so, we performed a Mantel test, to test for a correlation between the disparity matrices (D for each taxa pair) calculated for each of the modules' character sets, respectively, and the disparity matrix calculated for the whole character set. Taxa pairs for which D could not be computed for one of the functional modules, due to missing data, were pruned from the matrices before performing the test. To assess if the perianth, androecium and gynoecium of Ericales are more or less variable than the rest of the flower, we then compared the disparity ( In addition to the 381 extant species, nine ericalean floral mesofossils from the Cretaceous were added: Actinocalyx bohrii (Diapensiaceae; [64]), Glandulocalyx upatoiensis (Actinidiaceae or Clethraceae; [40]), Parasaurauia allonensis (Actinidiaceae; [65]), Paleoenkianthus sayrevillensis (Ericaceae; [37]), Paradinandra suecica (Pentaphylacaceae, Theaceae or Actinidiaceae; [38]), Pentapetalum trifasciculandricus (Theaceae; [39]), Raritaniflora glandulosa, R. sphaerica and R. tomentosa (Ericales; [36]). Morphospace and partial disparity were recomputed for the dataset including these fossils. Our dataset contains 12 512 data entries. In total, 1927 (13.4%) data are missing. The average percentage of missing data is 13.4 ± 10.3 (mean ± s.d.) per taxon and 13.4 ± 12.2 per character. The floral morphospace of Ericales is organized in a continuous cloud (figure 1c; electronic supplementary material, interactive three-dimensional figure S3). The first three principal coordinate axes of the space representation summarized 31.5% of the original variance (15.8%, 8.45% and 7.25% respectively, electronic supplementary material, figure S3), the first two, 24.2% (figure 1c). Both three-dimensional and two-dimensional representations gave a fair approximation of the relative dissimilarity among taxa (Pearson's r = 0.79, p < 0.001 for three axes; Pearson's r = 0.71, p < 0.001 for two axes). In this space, most of the seven suprafamilial clades, plus the families Theaceae, Lecythidaceae and Mitrastemonaceae, occupy distinct neighbouring regions arranged in a mosaic pattern (PERMANOVA: F = 23.24, r2 = 0.36, p < 10−4; table 1 and figure 1c; electronic supplementary material, figure S3). The only exception is the balsaminoid clade, which does not significantly differ from most of the other suprafamilial clades, mainly because two of its three families, Tetrameristaceae and Marcgraviaceae, overlap with most clades in the order (PERMANOVA: see electronic supplementary material, table S1).
Finally, each family cluster overlaps with at least two other families from its own and other suprafamilial clades (PERMANOVA: F = 18.85, r2 = 0.52, p < 10−4, electronic supplementary material, table S1 and figure S3). For instance, Ericaceae significantly differ from only 11 of the 21 other families and distinctly overlap with, e.g. Cyrillaceae and Marcgraviaceae (electronic supplementary material, table S1). The disparity of ericalean families ranged from Four ericalean families together contribute 50% of the Ericales disparity: Lecythidaceae (PDiv = 16.1%), Sapotaceae (PDiv = 14.3%), Primulaceae (PDiv = 14%) and Ericaceae (PDiv = 9.8%, electronic supplementary material, figure S4). We found a positive, but nonlinear, correlation between disparity ( Figure 2. Relation between disparity and (a) species richness (log transformed) and (b) stem age in ericalean families. Blue line: linear regression between disparity and log transformed species number (Spearman's rho = 0.14, p = 0.63). (Online version in colour.)
The nine fossils of Ericales included in our dataset group together near the centre of the morphospace (figure 1e). Their distribution in the space overlaps with the balsaminoids, the ericoids, the Pentaphylacaceae-Sladeniaceae clade and Mitrastemonaceae (electronic supplementary material, table S2) and they contribute 1.8% to the total disparity of Ericales. Note that the percentage of missing data in the fossil dataset was 10.51% (electronic supplementary material, table S2). The most variable characters (i.e. characters with a mean pairwise difference between all taxa Dchar ≥ 0.5) stem from the androecium (anther orientation, filament fusion to corolla and anther attachment) whereas the least variable characters (i.e. characters with a mean pairwise difference between all taxa Dchar ≤ 0.05) describe the perianth (number of petal whorls, petal phyllotaxis, sepal phyllotaxis, and perianth differentiation; electronic supplementary material, figure S5). Disparity for each of the three modules significantly increases with disparity of the whole character set (figure 3a–c): Mantel test for the perianth (71 631 taxa pairs): p < 0.001; for the androecium (72 390 taxa pairs): p < 0.001; for the gynoecium (70 876 taxa pairs): p < 0.001. Figure 3. Differential variation in the three functional modules of Ericales flowers. (a–c) Pairwise dissimilarities (D) for (a) the perianth, (b), the androecium and (c) the gynoecium, plotted against pairwise dissimilarities for the total character set. Black plain lines: y = x. Blue dashed lines: linear regressions between the plotted variables. The heat colour gradient indicates the density of dots in the graph. (d–f) Disparity (in blue) of (d) the perianth, (e) the androecium and (f) the gynoecium, plotted with the distributions of
The perianth varies significantly less than the rest of the flower (mean pairwise distances between all taxa for the perianth organs only: The floral morphospace of Ericales is organized in a continuous cloud, where most of the suprafamilial clades occupy distinct neighbouring regions arranged in a mosaic pattern. In other words, each of these lineages evolved towards a distinct combination of floral morphological traits. Our analysis also reveals two contrasting patterns of trait variation. On the one hand, ericalean species from across the order share recurrent combinations of character states: for instance, 75.1% of the 381 sampled species share structurally bisexual flowers with a differentiated perianth, a single whorl of sepals, and a single whorl of five petals (electronic supplementary material, table S2). Many of these common traits are likely plesiomorphic and evolutionarily constrained within the order. On the other hand, some floral traits, such as petal union and stamen and integument numbers, are highly variable in Ericales (see electronic supplementary material, table S2), although they have traditionally been considered stable within major angiosperm clades [35]. These two conflicting patterns are the main reasons for the pre-molecular phylogenetic placement of Ericales' taxa in more than 10 different angiosperm orders [44,49]. Within each suprafamilial clade, our analysis showed that there are two main patterns of space occupation by families. The balsaminoid, styracoid, sarracenioid and ericoid clades are all morphologically homogeneous, with the families overlapping (electronic supplementary material, table S1). For instance, the sarracenioids (Sarraceniaceae, Roridulaceae, and Actinidiaceae) are all characterized by, e.g. free petals to which stamens are not (or only basally) attached. The recovery of sarracenioids based on morphology is consistent with their phylogenetic relationships and illustrates a case where floral traits typically have a higher diagnostic value than vegetative traits [66]. In contrast with the pattern described just above, families of the polemonioid, primuloid and Pentaphylacaceae-Sladeniaceae clades occupy distinct regions of the morphospace (electronic supplementary material, table S2). For instance, in the polemonioids, flowers of Polemoniaceae and Fouquieriaceae significantly differ (see also [45]): Fouquieriaceae flowers have a free calyx and more than five stamens arranged in two whorls and free from the corolla, whereas Polemoniaceae flowers generally have a united calyx and a single whorl of five stamens that are more or less fused with the corolla. Consequently, Polemoniaceae and Fouquieriaceae were placed in various different orders before being recovered in molecular phylogenies [35,67]. A discussion about Ericales families whose placement in the phylogeny is not resolved is given in the electronic supplementary material. In Ericales, the most variable family, Lecythidaceae, is a pantropical family of trees (340 species). The least variable family, Mitrastemonaceae, is a root-parasitic family composed of only one Asian and one Central American species (it is thus the smallest family of Ericales). The positive correlation between disparity and the species richness of families can be explained by the fact that family size is highly variable in Ericales, ranging from two species in Mitrastemonaceae to 4010 in Ericaceae. Although disparity measured as mean pairwise distances is generally robust against differences in group sizes [60], a group containing a thousand times more species than another is very likely to be more diverse. Additionally, reproductive isolation is often due to differences in floral traits [68,69], hence a correlation between floral disparity and species number is, in many cases, to be expected. In a study on the disparity of Neotropical pollen morphotypes, Mander [70] found a similar positive correlation between disparity and family size, with some exceptions (e.g. Poaceae and Papilionideae). Such exceptions to the overall patterns are also revealed in our study and may be illustrated by comparing, for instance, Sapotaceae and Lecythidaceae. Although Sapotaceae are more than three times as speciose as Lecythidaceae, their flowers are only half as morphologically diverse (figure 1). These two families are both constituted mainly of tropical trees, and have similar stem ages [33,42]. Other factors might thus explain their contrasted disparity, such as different diversification events, polyploidization events, species distribution or ecology. For instance, the higher diversity of Lecythidaceae may be linked to specialized floral adaptations towards different functional groups of pollinators. Lecythidaceae have evolved highly elaborate and specialized pollination mechanisms involving bees, beetles and bats [71–73], whereas Sapotaceae appear to be characterized by more generalist insect or sometimes bat pollination systems [42,74]. Only four ericalean families together (Lecythidaceae, Sapotaceae, Primulaceae and Ericaceae) contribute 50% of the Ericales disparity. Some families, like Ericaceae and Pentaphylacaceae, display high partial disparity because they are widely spread in the space; i.e. are themselves highly variable. Alternatively, some families, like Sapotaceae and Balsaminaceae, display high partial disparity because they are distributed in the periphery of the space, i.e. they increase the overall disparity by adding new traits or trait combinations. Such traits are, e.g. two whorls of sepals in Sapotaceae, and distally united filaments, nectar spur and zygomorphic flowers in Balsaminaceae (see electronic supplementary material, table S2). Zygomorphy and the presence of nectar spurs might explain the peculiar pattern found in Balsaminaceae: the family displays very low morphological disparity in the investigated organizational floral traits (figure 1b), in spite of the fact that it is extremely speciose. Balsaminaceae is composed of two genera: Impatiens (1100 spp.) and Hydrocera (1 sp.). The high taxonomic diversity in Impatiens may result from rapid radiation during the Pliocene and Pleistocene, triggered by climatic fluctuations resulting in refuge areas [75]. This recent diversification might explain the low structural disparity in Balsaminaceae (although perianth shape and colour, not coded here, are highly variable [76]). Zygomorphy and nectar spurs are often considered to be key innovations associated with high diversification rates: they are often linked to speciation through increased specialization in pollination, e.g. through precise pollen placement on pollinator bodies [77,78], and spur length filtering for the most efficient pollinators [79]. The lack of correlation between disparity and the age of families is not surprising, as disparity does not necessarily steadily increase over time [80]. The nine fossils we included showed low levels of disparity and contributed little to the total disparity of Ericales. These fossils are not considered to be closely related to each other [8,40], and have been tentatively assigned to a number of different families, albeit all belonging to a group consisting of ericoids, sarracenioids, styracoids and Pentaphylacaceae-Sladeniaceae (the assigned positions of the fossils, according to original papers, are highlighted in figure 1a). They are close in age (72–94 Ma) to the Ericales' initial diversification (crown age: ca 104 Ma [33]) and the features they share (e.g. bisexual flowers, pentamerous and actinomorphic whorls of sepals and petals, free stamens, and superior ovaries) could thus represent plesiomorphies for the order. Compared with other types of fossils (e.g. permineralized fossils or compression/impression fossils), charcoalified fossils often show excellent preservation of morphological and anatomical features [8]. The availability of such fossils thus offers opportunities for future work on the extinct disparity of flowers, also at the level of the angiosperms as a whole. Our results indicate that, in Ericales, morphological variation differs considerably among the flowers' three functional modules, probably because of the different selective regimes they are submitted to. The perianth varied significantly less than the rest of the flower, as expected. Contrastingly, the androecium varied marginally more than the rest of the flower, and the gynoecium showed neither less nor more variation than the rest of the flower, although we expected it to show more variation, due to its complex organization. In our dataset, most species are characterized by actinomorphic, whorled, and pentamerous flowers (like most core eudicots [5,10]). In general, the perianth is structurally less complex than the reproductive floral organs [27] and the characters we could code for were mostly related to organ number and arrangement (Bauplan; [81]). In the perianth, these characteristics are likely spatially and functionally constrained during development and anthesis and mostly stable at higher taxonomic ranks, with one explanation being that the number of perianth organs and their arrangement depends on meristem size during early development [10]. Stamen number, however, appears to be much less constrained, ranging from two to several hundred in our dataset (see electronic supplementary material, table S2). Large stamen numbers can easily be accommodated even on a relatively small floral base as the filament bases are generally small [7]. Polystemony (i.e. flowers with more stamens than perianth organs) has apparently evolved along several separate lineages in Ericales [35]. In addition to the diversity in stamen numbers, ontogenetic patterns of androecium development and anthetic stamen arrangement are also particularly diverse in Ericales, including complex ring primordia with multiple stamen whorls and stamens arranged in fascicles [82]. With the exception of early diverging angiosperms, there are probably only few other groups of angiosperms with such labile and diverse patterns of stamen numbers and arrangement as the Ericales. It seems likely that this lability in the androecium has played a major role during the evolutionary history and diversification of the Ericales and has allowed the group to explore new evolutionary paths in connection with different functional groups of pollinators. Finally, the gynoecium, with its multiple functions and complex morphogenesis, is often considered the most complex module of the flower [7,32]. All Ericales in our dataset were syncarpous. Syncarpy is relatively stable in angiosperms [81] and is believed to have many advantages [83]. For instance, it allows for the presence of a centralized canal for pollen germination that allows a single pollen load on a stigma to potentially reach all the ovules of the same flower [83]. It has been proposed that syncarpy allows for higher levels of synorganization both within the gynoecium and also between the gynoecium and other floral organs [32], like in Balsaminaceae and Tetrameristaceae, where the syncarpous gynoecium is highly synorganized with the androecium [46]. Once syncarpy has evolved, it is therefore likely to remain stable, and it is a factor decreasing disparity in the gynoecium. However, other gynoecial traits such as ovary position, type of placentation, the number of ovules per carpel and the number of integuments are remarkably variable across Ericales (electronic supplementary material, table S2). These contrasting trends of variation might explain the lack of signal for more or less variation of the gynoecium. Overall, variability occurs less at the level of floral organization (e.g. organ number, organ arrangement), than at the level of floral construction (architecture, mechanical properties) and mode (traits like organ shape and colour). Floral mode often directly concerns interactions with pollinators and is generally highly variable, even at low taxonomic ranks [81,84]. Such traits could not be included in our analysis, as most of them would not have been applicable throughout the order. Our dataset is more representative of floral organization and construction, with some exceptions: the most variable characters in the perianth are the union of sepals and of petals (electronic supplementary material, figure S5), typically linked to pollination, allowing for the formation of corolla tubes and a canalized access to floral rewards [85], channelling pollinator movements so that they touch stamens and stigmas [77]. The dataset supporting this article has been uploaded as part of the supplementary material (electronic supplementary material, table S2). M.C., F.J. and J.S. designed the research; M.C., S.L., M.v.B., F.J., H.S. and J.S. generated the dataset; M.C. and S.G. analysed the data; M.C., S.L. and J.S. wrote the paper. In addition, all authors contributed to the writing of the manuscript and gave final approval for publication. We have no competing interests. This work was supported by the Austrian Science Fund (FWF P 250077-B16). We thank S. Sontag for helping with literature search, T. Palme for creating the three-dimensional pdf figure, D. Breedlove, P. Herendeen, D. L. Nickrent, T. Rodd and A. Weissenhofer for pictures, W. L. Crepet and K. C. Nixon for the authorization to reproduce pictures, and U. Schachner for proofreading. We thank A. Dellinger, E. Reyes, Y. Städler, B. Kezzim, A. van Holt and L. Carrive for data entry. In addition, we are thankful to L. Mander and an anonymous reviewer for their critical and helpful comments and suggestions. FootnotesElectronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3717514. © 2017 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited. References
Page 6Explaining the evolution of species’ ranges is fundamental to understanding how biodiversity is distributed and maintained [1–3]. Species' ranges are influenced by biotic (e.g. predation, parasitism and competition), and abiotic factors (e.g. climate) [4,5]. Yet, we still do not fully know how species’ geographical ranges evolve and what factors fuel range expansions [6,7]. Generally, species’ ranges are limited by the inability of populations at the range edge to adapt to environmental pressures before going extinct [6,8,9]. Range expansions often result in smaller population sizes at the range periphery and decreased genetic diversity [10–12]. Thus, peripheral populations’ adaptive potential is low and their risk of extinction high. Adaptive evolution that prevents extinction could occur in such populations via new mutations or gene flow [8,13]. However, the waiting time for adaptive mutations is potentially too long to rescue edge populations, and gene flow from conspecifics will most likely consist of alleles from the range centre [14,15], which may be poorly adapted to the range periphery [8,16]. Whether such gene flow has positive or negative effects will differ depending on the novelty of the habitat into which expansion occurs. When expansion occurs into relatively non-novel habitat, gene flow from other conspecific populations inhabiting similar environments can provide an increase in genetic diversity or adaptive alleles to foster local adaptation in peripheral populations [8,17]. By contrast, when expansion occurs into novel habitat, theoretical and empirical work has shown that gene flow from the centre of the range can have an opposite effect, generating an influx of maladaptive alleles that prevent local adaptation in peripheral populations [8,16]. An alternative to conspecific gene flow as a source of genetic variation is hybridization with a resident species. Although hybridization is often deleterious [18,19], it is sometimes beneficial [20,21]. In such cases, introgression of heterospecific alleles may provide populations at the range edge with a source of genetic variation [22,23], including the transfer of specific adaptive alleles from one species to another [20,24–28]. This can result in rapid adaptation by peripheral populations, allowing for further expansion into the novel habitat [29,30]. Whether hybridization plays an important role in range expansion remains an open question, especially in animals, as most tests of the hypothesis have been in plants [21,29,31]. Yet, evaluating hybridization's role in the evolution of species’ ranges is important for ascertaining hybridization's role in the origins and distribution of biodiversity. Indeed, understanding the relationship between hybridization and range expansion is increasingly important for practical reasons as evidence shows that global change is altering the distribution of animal and plants species around the world [32–35] and hybridization events could become more common as a result [36]. We addressed the above issues with two goals for this study. Using a population genetic approach, we: (i) ascertained whether encountering a novel environment might limit range expansion as theory predicts; and (ii) evaluated the potential role of hybridization in expansion into a novel habitat. To achieve these goals, we used Plains spadefoot toads, Spea bombifrons, as a model system. Spea bombifrons occupy a wide range throughout the southwestern and central United States (figure 1) and are thought to be ancestral to the central plains region [37]. After the most recent glacial retreat, S. bombifrons appears to have expanded its range northward [37] through grassland habitat similar to the ancestral region, with further northern expansion taking place in current populations [38,39]. By contrast, S. bombifrons also may have expanded their range southwestward into an entirely different biome: the desert [37]. Museum collections record S. bombifrons in the Southwest USA in the late 1800s, so this expansion is not contemporary. However, in some populations in Arizona, the relative abundance of S. bombifrons increased within the last 30 years [40]. Figure 1. Map showing species’ ranges and sampling locations of S. bombifrons and S. multiplicata. Central Oklahoma is the likely origin of the S. bombifrons range and grey circles represent S. bombifrons only sampling sites. White circles represent sites with both species sampled. Inset shows breeding S. bombifrons (top) and S. multiplicata (bottom) males. Sampling site key located in table 1. (Photos by David W. Pfennig; map by Travis Taggart.)
Spadefoot toads breed, and their tadpoles develop, in ephemeral ponds that potentially dry before the tadpoles successfully metamorphose. This putative southwestward expansion of S. bombifrons is therefore striking because a limiting environmental factor for these amphibians is ponds that last long enough for tadpole metamorphosis. Indeed, a congener, S. multiplicata (Mexican spadefoot toad), that is ancestral to the desert region has shorter developmental times that enable tadpoles to more likely metamorphose before their desert ponds rapidly dry [41]. Where S. multiplicata and S. bombifrons co-occur, they potentially hybridize and produce viable offspring. Female hybrids can backcross to both parent species (hybrid males are sterile; [42,43]), thereby generating introgression between the two species [42,44]. Critically, hybrid tadpoles develop faster than pure S. bombifrons tadpoles, resulting in a fitness benefit for the expanding species to hybridize in a dry, desert environment [45]. In fact, S. bombifrons females that occur in sympatry with S. multiplicata have evolved facultative mate preference where they prefer conspecifics when breeding in deep, long-lasting ponds, but switch their preference to S. multiplicata in shallow, ephemeral ponds [45]. Consequently, S. bombifrons females primarily contribute to the production of F1 hybrids and the incidence of hybridization increases with decreasing pond size [46]. Because S. bombifrons appears to have undergone two distinct range expansions, and because they hybridize with a resident species in the context of one of those expansions, S. bombifrons is ideally suited to address the issues raised above. We specifically compared the genetic effects of range expansion into a novel, desert environment with expansion into a non-novel grassland environment, and ascertained whether hybridization has potentially facilitated the expansion of S. bombifrons into the southwestern USA. To do so, we used microsatellites and population genetic analyses to: (i) investigate patterns of genetic diversity and population structure within S. bombifrons; and (ii) evaluate patterns of introgression between S. bombifrons and S. multiplicata. Because the desert provides a novel habitat, we expected populations there to suffer genetic effects from bottlenecks and reduced gene flow compared to populations in the non-novel grassland habitat. Our findings support these predictions. Our results also suggest that hybridization has enhanced genetic variation in populations of southwestern S. bombifrons; this hybridization might have enabled their expansion into novel habitat. We obtained 217 samples of S. bombifrons from 21 locations across the USA through collection efforts and museum samples (figure 1). Locations ranged from a grassland environment in Nebraska to a desert environment in Arizona and included one location in Arizona that did not overlap locally with S. multiplicata (i.e. allotopy) (table 1). Additionally, we obtained 93 S. multiplicata samples from three sympatric locations in Texas and Arizona (electronic supplementary material, table S1). Genotype data for the Arizona S. multiplicata individuals were previously reported in [37]. Adult specimens of S. bombifrons from Arizona sympatry and both S. bombifrons and S. multiplicata from Texas sympatry were used to ensure accurate species identification (see figure 1 inset photos). Museum sample IDs are provided in electronic supplementary material, table S2. We genotyped each sample using 10 polymorphic microsatellite markers that were previously shown to not be in linkage disequilibrium (electronic supplementary material, table S3; method details in electronic supplemental materials) [47–49]. We used the software Arlequin v 3.5.1.2 [50] to calculate observed and expected heterozygosity for each location (electronic supplementary material, table S4). Deviation from Hardy–Weinberg equilibrium was calculated with an exact test contrasting observed and expected heterozygosity in Arlequin using a Markov chain with a chain length of 1 000 000 and 100 000 dememorization steps. We then corrected for multiple testing using a sequential Bonferroni correction at α = 0.05 for each locus in each population. All of our loci were in Hardy–Weinberg equilibrium in at least 70% of our sampling locations (electronic supplementary material, table S4). To understand the impact of habitat novelty on range expansion, we first examined population structure of S. bombifrons across non-novel grassland and novel desert environments using the software STRUCTURE v. 2.3.3 [51]. For the STRUCTURE analysis we also included samples of S. multiplicata individuals from seven locations. Because we included S. multiplicata in this analysis, we used only the eight microsatellite loci that could be amplified in both species (electronic supplementary material, table S3). We implemented 100 000 burn-ins followed by 200 000 Markov chain Monte Carlo runs. We also used an admixture model with uncorrelated allele frequencies to avoid the risk of overestimating the number of populations and the LOCPRIOR model to provide the software with collection information for each toad to ensure the detection of subtle population structure. We started simulations with K values of 1–28, to reflect the 28 sampling locations (table 1; electronic supplementary material, table S1). For each K, we ran 10 simulations to check for consistency between runs, and used the log likelihood [51] and delta K method [52] to determine the most likely number of genetic populations (electronic supplementary material, figure S1). To confirm our results, we used all 10 loci and calculated FST and RST statistics [53,54] to measure genetic differentiation between S. bombifrons populations. Permutation tests (using 10 000 permutations) implemented in ARLEQUIN v. 3.5.1.2 [50] were used to determine whether pairwise FST and RST values were significantly different from 0. We also performed an analysis of molecular variance (AMOVA) and calculated p-values based on permutation tests (using 1000 permutations) in ARLEQUIN for the northern, central, and desert regions to examine differences in the level of population structure across the range. Finally, using RST statistics we performed a Principal Coordinate Analysis using the cmdscale function in R (v. 3.0.1). To understand if habitat novelty results in different colonization mechanisms and demographic effects, we used Poptools [55] to standardize sample sizes to seven individuals per collection site before comparing levels of genetic diversity (using the value 1-Qinter, the inter-individual diversity within populations), which were measured using Genepop v. 4.1.0 [56]. Locations with fewer than seven samples were excluded resulting in 17 sites used in the analysis. A sample size of seven was chosen to optimize statistical power with number of sites. We also calculated allelic richness using ADZE-1.0 [57], which uses a rarefaction approach to account for unequal sample sizes. For this analysis, we set the minimum sample size per site to seven individuals with no missing data across all loci. By setting the minimum samples size as seven individuals with no missing data, we optimized the number of sites used (12 total for this analysis) while minimizing bias that could result from including sites with too few individuals. We compared values between grassland and desert regions using a Welch two sample t-test implemented in R v. 3.1.2. To examine possible deleterious effects of a novel environment, we tested for recent population bottlenecks in S. bombifrons populations with at least 10 individuals samples using a Wilcoxon test [58] for heterozygosity excess across loci and a two-phase mutation model in the software Bottleneck [59]. We accounted for possible null alleles and determined corrected allele frequencies using the Brookfield 1 estimator [60] implemented in Micro-checker, v. 2.2.3 [61]. We also performed this analysis with uncorrected frequencies; however results were not qualitatively different so we only report results based on corrected frequencies. If novel habitats limit expansion, we expected signatures of bottlenecks in the desert, but not grassland, S. bombifrons populations. To calculate the amount of hybridization across the range, we calculated gene flow between the species in Texas and Arizona using likelihood ratio tests implemented in the coalescent-based software package MIGRATE-N 3.2 [62]. For this analysis we used the Brownian motion approximation to the ladder (‘stepwise’ or ‘one-step’) mutation model and Bayesian inference with multiple heating chains to jointly estimate parameters with three replicates [63,64]. This also allowed us to determine the directionality of gene flow. We identified specific loci showing patterns of introgression by comparing hybrids with both parental species using FST. Hybrids were identified as admixed individuals (more than 10% assignment to heterospecific populations) based on inferred ancestry by STRUCTURE. Using the FST calculations, we could detect when admixed individuals were more genetically similar to the heterospecific, indicating introgression at that locus. We further identified signatures of introgression at locus SpeaC7 with alleles primarily found in the desert toad, S. multiplicata appearing in S. bombifrons individuals. We resampled without replacement to obtain population sizes of 10 and examined changes in frequencies of the putative heterospecific allele across the S. bombifrons range. Allele frequency values at sampling sites were used to generate an allele frequency surface map by inverse distance weighted (IDW) interpolator in ArcGIS v. 10.4.1 (ESRI, Redlands, California, USA). IDW estimates values by averaging nearby data points, with closer points carrying more influence. Across grassland populations (table 1 and figure 1), STRUCTURE analysis revealed a clinal pattern of increasing membership to the ‘green’ group with increasing distance northward into Nebraska from Oklahoma, the putative centre of the range (figure 2). This is consistent with a northward range expansion as suggested by previous observational and genetic data [37–39]. Despite evidence of a range expansion, S. bombifrons populations showed stable levels of heterozygosity (F14,135 = 0.54, p = 0.91; figure 3a), genetic diversity (F10,99 = 1.26, p = 0.27; figure 3b) and allelic richness (F7,72 = 0.50, p = 0.83; figure 3c) across the grassland portion of their range from Texas to Nebraska (figure 3). These latter findings contrast with previous empirical and theoretical work in other systems showing that an expanding species will exhibit decreasing genetic diversity [11,12,65–68]. We further found that FST and RST values were low across much of the northern and central portions of the range (electronic supplementary material, table S5 and figure S2) of S. bombifrons. An AMOVA analysis using RST confirmed that differences among sites did not account for a significant amount of observed variation in the northern (Kansas, Nebraska, Colorado: 3.4% variation, p = 0.17) or central (Oklahoma, Texas: 5.18% variation, p = 0.06) populations. Thus, where S. bombifrons has expanded its range across similar habitat, it maintains stable levels of genetic diversity and allelic richness, likely via gene flow among other grassland populations. Figure 2. Structure plot showing that K (number of distinct populations) = 5 for S. bombifrons and S. multiplicata. The northern range of S. bombifrons shows increasing differentiation from the central range whereas the most southwestern portion of the range is highly differentiated. Individual toads are indicated by vertical bars and colour denotes population membership.
Figure 3. Genetic trends across S. bombifrons range. Sampling sites are ordered from north to southwest. (a) High levels of heterozygosity (Ho) are maintained throughout the range. (b) Genetic diversity (using the value 1-Qinter, the inter-individual diversity within populations) was also maintained. (c) Allelic richness declines in desert populations.
Contrary to the non-novel grassland populations, allelic richness decreased among desert populations in New Mexico and Arizona (t4.77 = 6.90, p < 0.01; figure 3c). Additionally, we found the desert-inhabiting populations to be highly genetically differentiated compared to grassland populations in our FST, RST and STRUCTURE analyses (figure 2; electronic supplementary material, table S5 and figure S2). Not only are the southwestern populations significantly differentiated from the northern and central populations, they are significantly differentiated from one another. An AMOVA analysis using RST confirmed that differences among sites accounted for increased population structure in the desert populations (New Mexico, Arizona; 9.33% variation; p < 0.001). Additionally, we found evidence for bottlenecks in two Arizona populations (Zent, p < 0.01 and Wilcox, p < 0.01). Outside of Arizona, we did not detect evidence for bottlenecks. The decline in allelic richness, high population structure, and evidence for bottlenecks is consistent with the desert habitat limiting gene flow and restricting movement among conspecific populations. Nevertheless, genetic diversity and heterozygosity appear to be maintained among desert populations of S. bombifrons as neither measure was significantly different from values found in grassland populations (genetic diversity: t4.66 = 2.54, p = 0.06; heterozygosity: t6.73 = 1.28, p = 0.24; however the marginal p-value for diversity could reflect insufficient power to detect a difference). To investigate if hybridization might have enhanced genetic diversity in S. bombifrons in the novel desert habitat, we surveyed S. multiplicata from two areas of sympatry—Texas and Arizona. We found that outlier S. bombifrons individuals located in Texas were genetically similar to S. multiplicata, indicative of hybridization (figure 2). Additionally, we found outlier S. multiplicata individuals in Arizona appearing genetically similar to Arizona S. bombifrons, again pointing towards hybridization. Rather than equal gene flow between the species, we found asymmetrical gene flow with the recipient species differing based on sympatric location. Migrate-n confirmed these findings, indicating a higher level of gene flow from resident S. multiplicata to S. bombifrons in Texas (S. multiplicata → S. bombifrons 5.43 immigrants/generations; S. bombifrons → S. multiplicata 1.40 immigrants/generation) but a higher level of gene flow from invading S. bombifrons to native S. multiplicata in Arizona (S. bombifrons → S. multiplicata 4.09 immigrants/generation; S. multiplicata → S. bombifrons 1.89 immigrants/generation). Additionally, individual examination of the markers revealed introgression across multiple loci (electronic supplementary material, table S6). For example, in sympatry, admixed individuals were genetically indistinguishable from heterospecifics at multiple loci, while showing significant differentiation from their conspecific population. For locus SpeaC7 in particular, a S. multiplicata allele was not found in high frequencies in central or northern S. bombifrons populations outside of sympatry, but was maintained at a relatively high frequency in S. bombifrons throughout the desert habitat (figure 4). Given that we detected introgression across multiple loci with only a handful of microsatellite markers, it is possible that hybridization has introduced a larger amount of genetic variation than observed here. Figure 4. Introgression and maintenance of heterospecific allele in S. bombifrons. ArcGIS frequency surface map of heterospecific (S. multiplicata) allele at locus SpeaC7 over the range of S. bombifrons based on observed population frequency at collection sites (black pentagons). See Methods for details. This heterospecific allele first appears in Texas sympatry and is maintained throughout the desert region. It is present even in the most westward S. bombifrons desert population, which is allotopic.
We used Plains spadefoot toads, S. bombifrons, to examine the population genetic effects of habitat novelty during range expansion and to evaluate how hybridization with a resident species can impact genetic variation during expansion into novel habitat. In contrast to populations in the ancestral grassland environment, we found that the novel desert environment was associated with reduced gene flow and recent population bottlenecks in S. bombifrons. Although these factors led to a reduction in allelic richness, hybridization with the resident S. multiplicata appears to have led to the transfer of genetic variation in the novel desert habitat. Such transfer of genetic variation via hybridization might have facilitated the range expansion of S. bombifrons into a novel habitat. Generally, species are expected to evolve expanded ranges when edge populations adapt to local conditions and become sources of dispersers [6,8,9,69]. Adaptability of peripheral populations therefore sets the limits of a species’ range. A key factor that limits adaptability is genetic diversity: in the absence of genetic diversity, populations are unable to evolve in response to local selective pressures [70,71]. Ironically, a common signature of range expansion is reduced genetic diversity because edge populations are often the result of serial founder events or suffer population crashes (and, concomitantly, genetic bottlenecks) [10–12]. Although dispersal and the resulting gene flow among conspecific populations can reintroduce genetic variation into peripheral populations [72,73], such gene flow can inhibit adaptation if alleles from the range centre are maladaptive at the range edge [8,16]. Moreover, dispersal might be limited across novel habitat [74–76]. Therefore, the novelty of the habitat into which expansion occurs might critically impact the adaptability of peripheral populations. Our results illustrate this dynamic between range expansion and habitat novelty. In the northward range expansion into a relatively non-novel grassland habitat, allelic richness and genetic diversity levels remain high among S. bombifrons populations (figure 3). Additionally, FST and RST values throughout the grassland regions are relatively low, suggesting ample movement of toads throughout the grassland range, which maintains a high level of genetic diversity and a low level of population differentiation (electronic supplementary material, table S5). Conversely, the southwestward range expansion by S. bombifrons into novel desert habitat showed a different pattern. Desert populations were highly genetically differentiated not only in comparison with the rest of the range (figure 2), but also between closely located populations in the desert itself (electronic supplementary material, table S5). Such strong population differentiation was likely influenced by bottleneck events. The differences we observe between the northern and southwestern range expansions by S. bombifrons highlight how expansion into novel versus non-novel habitats can generate variation in gene flow and population genetic patterns. Indeed when a species expands into a similar habitat, gene flow among populations is more likely to result in maintenance of allelic richness and adaptive potential [8,17]. By contrast, our results indicate that a novel habitat might restrict movement, so that populations are less likely to receive migrants (and genetic rescue) and are more likely subject to population crashes and extinction. The novel habitat thus generates a negative genetic impact, which—without some counterforce—could limit peripheral population adaptability. One such counterforce is hybridization with a resident species that is locally adapted [21,30]. In spadefoots, hybridization with desert-adapted S. multiplicata may be one way in which S. bombifrons maintains genetic variation in the novel habitat. Patterns of introgression are consistent with this, as is our finding that desert populations of S. bombifrons contain high levels of genetic diversity despite evidence of recent bottlenecks. Specifically, we examined two areas of sympatry to evaluate the effects, if any, of hybridization on S. bombifrons desert expansion. Interestingly, we found that directionality of introgression differed between these regions (figure 2). In Texas, introgression was from the resident S. multiplicata into the expanding S. bombifrons, whereas the opposite was true in Arizona. Although theory predicts that the rarer, expanding species should receive massive introgression from the resident species, the female driven hybridization and hybrid male sterility seen in S. bombifrons may result in the inverse pattern in Arizona (sensu [77]). Additionally, a study examining hybrid mate choice in this system found that hybrid females had no mate preference between parental species [78]. Thus, the relative abundance of parental species might drive patterns of backcrossing and introgression in a given region [78]. Indeed in Texas, S. bombifrons are more common, resulting in a higher likelihood for hybrid females to mate with S. bombifrons, thereby moving S. multiplicata alleles into the S. bombifrons population. In Arizona, by contrast, there are relatively fewer S. bombifrons, so hybrid female behaviour might contribute to introgression of S. bombifrons alleles into the S. multiplicata population. Whether hybrid behaviour contributes to the patterns of introgression we observed requires further study. Nevertheless, our data suggest that S. bombifrons has received S. multiplicata alleles in Texas, where they initially encountered S. multiplicata. Furthermore, we found that S. bombifrons has not only received heterospecific alleles, but has maintained an introgressed allele following their southwestward expansion deeper into the desert region (figure 4). This allele is maintained despite evidence of bottlenecks and low gene flow among Arizona populations of S. bombifrons (but see below), and despite introgression primarily from S. bombifrons into S. multiplicata in Arizona. We also observed a relatively high frequency of this heterospecific allele in an allotopic S. bombifrons population at the western edge of the species’ range where S. multiplicata is absent. This latter result emphasizes that ongoing hybridization is not necessary to maintain this genetic variation. Why this allele persists is not clear, but one explanation is that it is linked to a functional locus under selection. Such a pattern would be expected if S. bombifrons acquired adaptive alleles from S. multiplicata that enabled them to expand into the desert habitat. Although this explanation is speculative at this point, we can conclude that the allele is not being purged from the populations as would be expected if it were associated with reduced fitness hybrids. The exception to this pattern of heterospecific allele maintenance was the sympatric Zent population in Arizona in which the heterospecific allele is absent. However, Zent also shows one of the strongest signatures of a recent bottleneck and is a newly discovered (and possibly newly established) site with very low population size (N = approx. 12–20 adults of both species in recent samplings). The allele may therefore have recently been lost through genetic drift. Although it is possible that this putatively heterospecific allele was a result of convergent evolution or shared ancestry, rather than introgression, both possibilities seem unlikely. Texas S. multiplicata samples have a high frequency of this allele (approx. 50% in some populations), making it likely to be shared during hybridization. Additionally, we do not see significant frequencies of this allele in any S. bombifrons populations north of the Texas sympatric zone. Given the spatial pattern of the allele frequency (figure 4), and that hybridization between these two species occurs [46], the most parsimonious explanation for its presence in S. bombifrons is introgression. As S. bombifrons expanded into the novel desert environment of the southwestern US, the receipt of S. multiplicata alleles in Texas could have provided S. bombifrons with adaptive genetic variation that enabled them to colonize the novel habitat and further expand southwestward. Future work examining differential levels of introgression across the genome and its adaptive significance, if any, are underway to evaluate this possibility. Regardless, this study indicates that hybridization with S. multiplicata has altered the population genetics of S. bombifrons. Our findings suggest that hybridization with a resident species may be a way in which expanding species can maintain levels of genetic diversity in a novel habitat, which could enable further expansion. Given shifting species’ ranges [32–35,79] and the likelihood that hybridization will become increasingly common [36], the need to evaluate hybridization's role in range expansion is more pressing now than ever [30]. Data are provided as electronic supplementary material. K.S.P. and A.A.P. conceived of the project and its design. K.S.P. collected samples in addition to those provided by the museums acknowledged below. A.A.P. and R.G. genotyped samples and analysed data. A.M.R. provided microsatellite data for S. multiplicata. A.A.P. wrote the paper in collaboration with K.S.P. R.G. and A.M.R. approved and edited the final versions. The authors declare no competing interests. This work was supported by a grant from the National Science Foundation (IOS-1555520) to K.S.P.; a UNC Tony and Elizabeth Long Research Award to R.G.; and grant K12GM000678 from the Training, Workforce Development and Diversity division of the National Institute of General Medical Sciences (NIGMS), National Institutes of Health (NIH) to A.A.P. We thank the following museums for S. bombifrons samples: Museum of Vertebrate Zoology, University of California, Berkeley; Sam Noble Oklahoma Museum of Natural History; University of Kansas Herpetology Tissue Collection; and Sternberg Museum of Natural History, Fort Hays State University (FHSM). We are also grateful to David Pfennig and Audrey Kelly for comments on the paper. FootnotesElectronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3723952. References
Page 7The current species richness of a group of organisms results from the diversification process that has occurred throughout its evolution. In plants, a variety of intrinsic and extrinsic factors affect the diversification process [1]. Among those factors, changes in climatic conditions [2], the colonization of new geographical areas [3] or the evolution of particular traits might create new possibilities for species diversification [4,5]. In angiosperms, traits such as biotic pollination, floral symmetry and nectar spurs, which are all related to specialized pollination and the ability to generate reproductive isolation, have been proposed as key innovations due to their positive effects on diversification [6,7]. The role of specialized biotic pollination in the diversification of angiosperms is a long-standing question [8], but the mechanisms that led to the apparent association between pollination and species richness are still rather unexplored [9]. One hypothesis is that diversification in angiosperms has been enhanced by the effect of pollinator specialization on reproductive isolation. Spatial and temporal differences in the availability of the most effective pollinator across the species range could produce pollinator shifts, floral divergence, reproductive isolation and, ultimately, speciation in plants [10]. Evidence for pollinator-shift effects in plant speciation have been found for Costus [11], Gladiolus [12] and Lapeirousia [13], and a review of available species-level phylogenies estimated that around 25% of the divergence events could be associated with pollinator shifts in angiosperms [14]. Although these results suggest that frequent pollination shifts have occurred during the speciation events in angiosperms, a large proportion of these events could still occur within specific pollination systems. Indeed, an alternative hypothesis proposes that diversification rates in angiosperms increase with specialization on certain guilds of pollinators, rather than with pollinator shifts per se [12]. For example, vertebrate pollination, and in particular pollination mediated by birds, is associated with plant species richness in various clades [15,16]. The evaluation of the role of ornithophily in the diversification of the whole Gesneriaceae family has recently indicated distinct patterns between Old and New World lineages [17]. The evolution of bird pollination (specifically hummingbird pollination) was associated with an increase in diversification rates in the New World, while no influence was detected for the lineages in the Old World. However, a necessary step to further understand the effects of hummingbird pollination in the New World plant diversity is to evaluate whether the frequent shifts among pollinator groups or the specialization on hummingbird pollination is influencing plant diversification [18]. Surprisingly, the relative contribution of these two processes remains unexplored. The aim of this study is to evaluate the tempo of evolution of functional groups of pollinators, in particular hummingbirds, and their impact on the diversification rates of the Neotropical lineage of the family Gesneriaceae, hereafter referred to as Gesnerioideae, which is the lineage of this family with hummingbird interactions. Specifically, by expanding the most recent phylogenetic sampling by 129 species, we conducted an accurate evaluation of the pollination syndromes, their evolution and their impact in diversification rate shifts for the subfamily. The Gesnerioideae is a clade of herbaceous plants, shrubs or more rarely small trees. It contains 75 genera and over 1200 species found exclusively in the Neotropics, with the exception of few Southwest Pacific taxa in the tribe Coronanthereae [19,20]. Molecular dating and biogeographic reconstructions have estimated an origin of the Gesnerioideae in South America during the Early Oligocene, with a rapid range expansion into most Neotropical regions [21]. The species in this subfamily exhibit a large diversity of floral morphology associated with repeated adaptations to different pollinators, such as hummingbirds, bees and bats [22–27] (figure 1). Therefore, this clade is particularly interesting to test the mode and tempo with which plant–pollinator interactions have evolved and how they have influenced species diversification. Figure 1. (a) Examples of plant–pollinator interactions in the Gesnerioideae (images 1–7); photo information in the electronic supplementary material, table S2. Photo credits: images 1, 2, 4, 5, 7 by I. SanMartin-Gajardo; image 3 by a. Weber; image 6 by L. Freitas. (b) Discriminant analysis conducted for 118 plant species and nine floral traits. Numbers refer to species shown in (a). Electronic supplementary material, table S2, provides the morphological data and the source of pollinator observations for each species. (c) Bayesian common ancestor (CA) phylogenetic reconstruction showing one stochastic mapping of pollination syndromes. White and gray boxes correspond to taxonomic tribes and subtribes. Names following classification by [20]: N, Napeantheae; C, Coronanthereae; T, Titanotricheae; Sp, Sphaerorrhizinae; O, outgroups. Colours on branches correspond to pollination syndromes: blue, insect; red, hummingbird. Trait states: blue, insects; red, hummingbirds; green, bats. Grey concentric circles have 10 Myr span.
Here, we reconstructed one of the largest species-level phylogenies for a group of Neotropical plants based on four DNA loci and a wide sampling of Gesnerioideae species to test for temporal variations and trait-dependent rates of diversification at a continental scale. We specifically investigated (i) the number and timing of transitions among pollination syndromes in the subfamily, and (ii) the temporal match between the evolution of hummingbird pollination and hummingbird diversification in South America [28]. Additionally, we tested (iii) whether the evolution of hummingbird pollination has contributed to the elevated Gesnerioideae diversity in Neotropics, and (iv) whether diversification rates were associated with recurrent shifts of pollinators or the observed species richness was driven specifically by hummingbird-mediated pollination. Addressing these questions in such a large and diverse group of plants will contribute to a better understanding of how ecological factors shape current patterns of species richness in the Neotropics [15]. Our taxonomic sampling consisted of 583 species representing all the 75 recognized genera in Gesnerioideae and about 50% of the species in the subfamily [20]. The sampling of each tribe and outgroups is detailed in the electronic supplementary material, methods S1. A total of 475 sequences were amplified from field samples for this study (see ‘Data accessibility’ below) and merged to available Genbank sequences. Sequences were aligned using MAFFT v. 7 [29]) and all sites were scored for accuracy of the alignment using Guidance [30]. We identified the best substitution model for each DNA region using the Akaike information criterion (AIC) as implemented in the phymltest function in R (ape package [31]; see details in the electronic supplementary material, methods S1). Relationships among species were reconstructed by Bayesian inference using MrBayes v. 3.2 [32]. Data partitions and details of Bayesian inference are described in the electronic supplementary material, methods S1. Divergence times were estimated using a relaxed clock model with uncorrelated lognormal prior distribution for the rates of substitution and a birth–death prior for the age of each node as implemented in BEAST v. 1.7.0 [33]. Secondary calibration was performed by imposing priors for the divergence times for the clade containing all Gesneriaceae (including Sanango racemosum and members of the Didymocarpoideae family). MCMC settings and tree sampling are described in the electronic supplementary material, methods S1. The predictability of pollination syndromes is largely debated [34,35]. However, a recent meta-analysis supported the concept of pollination syndromes, especially for tropical plants [36], and encouraged the use of floral characters as a proxy for pollination interactions in macro-evolutionary studies [16,17]. In Gesnerioideae, several studies combining field observations and multivariate analyses of morphometric data have demonstrated that suites of floral traits could predict specialized pollination by hummingbirds, bees and bats in Drymonia [27], Gesnerieae [25], Nematanthus and Codonanthe [37], and Sinningieae [24]. To further test the validity of pollination syndromes, we assessed thoroughly the correlation among floral traits and functional groups of pollinators among the species of Gesnerioideae with documented pollination systems. An extended bibliographic search was conducted to identify all species with published information about their pollinators (electronic supplementary material, table S2). Flowers of these species were characterized using nine morphological traits reflecting their variation in size, shape and colour (electronic supplementary material, table S2). Trait values were derived from published morphometric datasets, monographic revisions and our own measurements of flowers collected in the field or in living collections, or from scaled images available on J.L.C.'s website (www.gesneriads.ua.edu). Among these traits, the degree of corolla constriction (i.e. tubular versus bell-shaped corolla) and the presence of pouched or urn-shaped corolla have been identified as key traits to discriminate hummingbird- from bee- and bat-pollinated flowers in different groups of Gesnerioidae [24,25,27,37]. Experimental results have also demonstrated the role of flower constriction and anther exsertion in improving the morphological fit between hummingbirds and flowers and/or in deterring less efficient pollinators like bees [38]. We used a discriminant analysis to maximize the differences in each trait among functional groups of pollinators (i.e. hummingbirds, bats, insects and generalists), and to estimate their predictability for the identification of pollination/shape associations. We used the lda and predict functions from the R package MASS [39]. The most discriminant floral traits were then used to predict the functional groups of pollinators for the species included in the phylogeny that lack direct observation of pollinators. In all subsequent analyses requiring binary states (see below), the bat-pollinated species (8 out of 590 species) were merged into the hummingbird pollination syndrome category. We based this choice on the fact that (i) hummingbirds and nectarivorous bats are both vertebrates with hovering ability, (ii) certain bat-pollinated species are generalists (pollinated also by hummingbirds during late afternoon and at dawn [25]), and (iii) according to a three-state stochastic mapping analysis, most of the bat-pollinated species in Gesnerioideae evolved recently from hummingbird-pollinated species [40] (see the electronic supplementary material, figure S5). The study of trait evolution has largely improved by considering evolutionary time in the modelling of a trait change [41], and by including the species diversification process itself in binary-state speciation and extinction (BiSSE) models [42]. Here, we incorporate most of these improvements by jointly modelling the evolution of pollination syndromes and trait-dependent diversification rates (binary-state trait). For this, we used estimates of transition rates between hummingbird and insect pollination syndromes from the BiSSE model that decomposes the evolutionary process into state-specific speciation and extinction rates and two transition rates. We performed an ancestral state reconstruction, accounting for the influence of diversification, using the asr function from the R package diversitree [43] to estimate the marginal probability of each state at each node. This function is only available for the BiSSE model, and not the other extensions of this model, such as the cladogenetic state change speciation and extinction (ClaSSE) model used in the diversification analysis. The temporal assessment of insect and hummingbird-adapted flowers was done by mapping changes in pollination syndromes across the Gesnerioideae phylogenetic tree. We incorporated the BiSSE estimates of ancestral states into the stochastic mapping (modifying the simmap function in the R package phytools [44]; script available from ‘Data accessibility’ below) and ran 200 reconstructions on independent trees. For each stochastic mapping, we divided branch lengths into time bins of 1 Myr and recorded the number of transitions from and to hummingbird pollination syndrome in each bin. We reported the time bin at which 95% of the stochastic mappings have at least one transition event as the onset time for each type of transition. We performed an additional three-state stochastic mapping without considering trait-dependent diversification to explore the evolution among hummingbird, bat and insect pollination syndromes (see the electronic supplementary material, methods S2). First, we tested whether diversification rates were constant or varied through time using the R package TreePar [45]. The model settings are described in the electronic supplementary material, methods S3. Second, we tested a range of trait-dependent diversification models to assess correlations between evolution of pollination syndromes in Gesnerioideae and changes in speciation and extinction rates. Those models included the BiSSE [42] and ClaSSE [46] classes. These models allow us to distinguish whether diversification rates are associated with any particular pollination syndrome (i.e. BiSSE) or whether they change in response to shifts in pollination syndromes (switches between insect and hummingbird syndromes, i.e. ClaSSE). We compared the BiSSE and ClaSSE models with a recently proposed trait-independent model where unobserved states, which are independent of our pollination syndromes, account for differences in the diversification process (called CID2 [47]). A binary trait was used to represent the pollination syndromes (insect as state 0; hummingbird as state 1). Sampling fraction was accounted in all trait-dependent and CID2 models as 0.53 and 0.48 for insect and hummingbird pollination, respectively. We estimated different speciation rates in BiSSE to obtain the parameters λ0 and λ1 for the speciation associated with insect and hummingbird pollination syndromes, respectively. For the ClaSSE model, we denoted the speciation within pollination syndromes (λ000, λ111), the speciation associated with a switch in trait for one of the descendant species, namely from insect to hummingbird pollination syndrome (λ001, λ101), and the speciation rates associated with switch in both descendant species (λ011, λ100). Each model also included two state-specific extinction rates (μ0, μ1) and two transition rates (q01, q10). We compared eight BiSSE and 13 ClaSSE models using maximum-likelihood estimates in the R package diversitree [43]. The best model was selected based on AICc, and we estimated the posterior distributions of each parameter for the best model in a Bayesian framework [5]. Priors and MCMC parameters are described in the electronic supplementary material, methods S3. Methods associating traits and diversification should be taken with caution [48]; these issues were minimized in our dataset and analyses (see the electronic supplementary material, methods S4), rejecting that rates may vary over the tree or through time [49]. The best models of molecular evolution were GTR + Γ and GTR + Γ + I for the nuclear and chloroplast DNA partition, respectively (electronic supplementary material, figure S1). The MrBayes and BEAST analyses resulted in congruent topologies. Our phylogenetic reconstruction (figure 1c) constitutes one of the largest species-level phylogenetic analysis for Neotropical plants. The topology corroborates the formal previous classifications [20], namely that Gesnerioideae comprises five tribes and 12 subtribes (posterior probabilities > 0.99; electronic supplementary material, figure S4). Relationships among tribes had a high support (posterior probabilities > 0.99), except that Titanotricheae, Napeantheae and Beslerieae (composed of the genera Besleria, Gasteranthus, Reldia, Cremosperma, Shuaria, Anetanthus and Tylopsacas) formed a clade (PP = 0.508 in the BEAST MCC tree) sister to the rest of the Gesnerioideae. The tribe Coronanthereae was sister to the Gesnerieae in agreement with prior results [19]. Five highly supported clades were resolved in the Gesnerieae, corresponding to the subtribes Gesneriinae, Gloxiniinae, Columneinae, Sphaerorrhizinae and Ligeriinae. Generic and infrageneric relationships largely agree with previous phylogenetic results obtained for these lineages [37,50–53]. Out of 74 genera of Gesnerioideae, eight appeared non-monophyletic and are still in need of further taxonomical revision (Achimenes, Diastema, Gesneria, Mandirola, Paliavana, Phinaea, Sinningia and Vanhouttea). Overall, 118 species with documented pollination systems were recorded from the literature (electronic supplementary material, table S2). Among them, 82 species were pollinated by hummingbirds, 19 species pollinated by bees, three species pollinated by other insects (butterfly, diptera and moth), and seven species pollinated by bats (electronic supplementary material, table S2). Seven other species are pollinated by a mix of nocturnal and diurnal visitors (e.g. hummingbird, bat and moth). These generalist species of Gesnerioideae have so far been recorded only on the Caribbean islands in pollinator-depauperate environments [26,54]. The discriminant analyses explained a large proportion of the floral trait variability (axes 1 and 2 with a 76.35% and 22.34% of explained variance, respectively). Linear discriminant axis 1 had a positive loading for corolla tube shape and lateral compression, and a negative loading for corolla length and the corolla width at mouth. Linear discriminant axis 2 had a positive loading for most of the traits except corolla length and tube shape (electronic supplementary material, table S3 and figure S2). The predictability of each group of functional pollinators was high (hummingbirds = 0.974, insects = 0.954, bats = 1.00, generalists = 0.66), and their separation in the morphological space was clear (figure 1b). Only five species of 118 (approx. 5%) have a group predictability lower than 0.8; these are one bee-pollinated (S. villosa), three generalists (G. viridiflora, R. leucomallon, R. vernicosum) and one hummingbird-pollinated (P. sericiflora, a species with flower morphology related to the bat syndrome but effectively pollinated by hummingbirds [22]). The standardized coefficients of each trait determine the contribution of the respective trait to the discriminant function among the groups. Based on these values (electronic supplementary material, table S3), we selected tube shape and lobe symmetry (electronic supplementary material, figures S2 and S3) as a proxy to assign pollination syndrome for species in the phylogeny whose pollination biology is unstudied in the field. Using this approach, and the information listed in the electronic supplementary material, table S2, we inferred 351 species with hummingbird pollination syndrome, eight species pollinated by bats and 231 species with insect pollination syndrome among the 590 taxa included in our phylogenetic tree (electronic supplementary material, table S1). The BiSSE estimates of transition rates between pollination syndromes indicated a median rate from insect to hummingbird pollination syndrome of 0.009, and from hummingbird to insect pollination syndrome of 0.044. Our stochastic mapping showed that pollination syndromes evolved on average from insect to hummingbird 31.50 (±10.07) times. Transitions to hummingbird pollination syndromes first occurred around 18.5 Ma and then increased in frequency over time (figure 2a). These transitions were reconstructed at or near the crown of large clades of Gesnerioideae, such as Besleria, Ligeriineae, Gesneriinae, Gloxiniinae and Columnea + Glossoloma genera (figure 1c). Reversions from hummingbird to insect pollination were highly frequent (on average 76.50 ± 18.06 times). These reversions to insect pollination started around 12.5 Ma and were mainly reconstructed on terminal branches or within clades including few species (figures 1c and 2b). Our three-state reconstruction treating bat and hummingbird syndromes separately showed that transitions to bat-pollinated flowers (all bat-pollinated species in Gesnerioideae have been observed in the field) have occurred at least seven times, since around 9.5 Ma, and mainly from hummingbird-adapted flowers (electronic supplementary material, figures S4 and S5). Figure 2. (a,b) Estimated number of transitions through time for pollination syndromes. Numbers below the pictograms correspond to the mean total number of transitions between the states and the standard deviation. Stars denote the starting point in time where at least one transition is recorded in 95% of the reconstructions. Grey bar is the age of the most recent CA of extant hummingbirds (20.3–24.7 Ma [28]). (c) State-dependent speciation rate estimates from ClaSSE model. λ = speciation rate (specific parameters in Diversification analysis section), 0 = insect, 1 = hummingbird. Colours and pictograms correspond to the binary pollination syndrome: blue for insect, and red for hummingbird states.
Our analyses of temporal shifts during the diversification of Gesnerioideae detected a single shift in diversification rate (p-value < 0.001; electronic supplementary material, figure S6) that most probably occurred around 18.5 Ma (95% confidence interval = 5.0–25.5 Ma). The mean net diversification rates were 0.067 and 0.177 Ma−1 for the periods before and after the shift, respectively. The comparison of models for trait-dependent diversification indicates that the model with the best AICc is the ClaSSE with all different speciation rates (though λ011 and λ101 constrained to be zero), extinction rates constrained to be zero and equal transition rates (table 1; see complete electronic supplementary material, table S7). We clearly rejected a trait-independent process shaping the diversification of the Gesnerioideae (CID2 null model with a δAICc = 25.0 from the best model). The MCMC parameter estimates suggested that speciation rates within pollination syndromes were higher than those associated with shifts between them (figure 2c). Furthermore, species within the hummingbird pollination syndrome have at least a twofold higher rate (mean λ111 = 0.252 Myr−1, 95% HPD = 0.193–0.314) than species within the insect pollination syndrome (mean λ000 = 0.102 Myr−1, 95% HPD = 0.071–0.133). All speciation rates associated with a shift in pollination syndrome, regardless of the direction of the shifts, are lower and close to 0.01, and thus an order of magnitude lower than the rates within pollination syndrome (λ001 = 0.006, λ011 = 0.003, λ101 = 0.004, λ100 = 0.010; figure 2c). Posterior distributions of extinction rate showed a higher extinction rate for hummingbird pollination syndrome species (mean μ0 = 0.015, 95% HPD = 0.000–0.041, and μ1 = 0.027, 95% HPD = 0.000–0.077). Transition rates between pollination syndrome states supported a higher rate of reversals to insect pollination syndrome (mean q01 = 0.006, HPD = 0.000–0.009, and mean q10 = 0.023, HPD = 0.000–0.040), and were of similar magnitude to the rates estimated by BiSSE (electronic supplementary material, table S4).
We found that diversification results are robust to the misidentification of functional groups of pollinators at the tips of the phylogenetic tree. First, the difference in rates of speciation between the two pollination syndromes (λ000 and λ111) is persistent if we remove the species-rich genus Columnea, which includes exclusively hummingbird pollination syndrome species (electronic supplementary material, figure S7b), and the few bat-pollinated species (electronic supplementary material, table S6). Second, the test for possible misidentification of functional groups of pollinators indicated that our estimation of speciation and extinction rates are extremely robust to up to 10% of equivocal states (for both insect and hummingbird) and that even 15% of misidentification leads to qualitatively similar results (electronic supplementary material, figures S8–S10). Finally, the simulations of traits, whose evolution is independent from the diversification process, showed that the estimated speciation rates within pollination syndromes (i.e. λ000 and λ111) are equal, under the null hypothesis, as well as the extinction and transition rates (electronic supplementary material, figure S7a). The effect of hummingbird pollination syndrome on diversification that we detected is thus not likely to be due to a particular shape of the phylogeny that could lead to a false detection of an association between traits and speciation [49]. We showed that hummingbird pollination probably played a role in the diversification dynamics of Gesnerioideae in the Neotropics. Two lines of evidence support this result. First, the diversification of this subfamily increased substantially around 20 Ma. This period corresponds closely to the dating for the common ancestor of hummingbirds into South America (22.4 Ma, 95% HPD: 20.3–24.7 Ma in [28]; 24–25 Ma in [55]) and the first appearance of plant species with hummingbird pollination syndrome in the Gesnerioideae (figure 2a). Second, we clearly show that species with hummingbird pollination syndrome have higher rates of speciation compared with species with insect pollination syndrome. On the other hand, we did not find high speciation rates associated with transitions between pollinator syndromes. These results indicate that species richness in this plant group has been driven by speciation within hummingbird-pollinated lineages, without involving shifts among functional groups of pollinators, contradicting the classical pollinator-shift hypothesis [14]. Our study suggests that Gesnerioideae was ancestrally pollinated by insects and that at least 31 transitions to hummingbird and bat pollination syndromes occurred during its evolution (figure 1c; electronic supplementary material, figure S4). The repeated evolution of hummingbird pollination syndrome in independent Gesnerioideae lineages centred in different geographical areas, such as the Brazilian Atlantic forest, Andes, Caribbean islands and Central America [21], is indicative of the success of this ecological interaction in multiple biomes of the Neotropics. We found also frequent state reversals from hummingbird to insect pollination syndrome providing evidence against the hypothesis that the evolution of hummingbird pollination could act as a dead end from which reversals to insect pollination are rare or no longer possible [14,56–58]. The convergent evolution of distinct floral morphologies and pigmentation in relatively short periods of time, as well as the reversibility of this system, are striking and encourage the investigation of whether a shared molecular basis might control these phenotypic transitions [59]. Flowers corresponding to the hummingbird pollination syndrome appeared in Gesnerioideae around 18.5 Ma (figure 2a), when hummingbirds were already present and diversifying in South America [28,55]. This rather early origin of hummingbird flowers, and the inferred south American origin of the Gesnerioideae [21], indicate that this plant group could be one of the oldest to have established interactions with the first hummingbirds living in the tropical regions of South America, unlike more recent hummingbird-adapted plant lineages [15,16,55,60]. This old interaction might have given the Gesnerioideae species enough time for the evolution of hummingbird pollination in separate lineages, and a substantial amount of transitions back to insect pollination. Our finding of a preferred trait-dependent model indicates an effect of pollination syndromes in the diversification of the Gesnerioideae. A twofold increase in speciation rates for species with a hummingbird pollination compared with insect pollination syndrome suggests that floral morphologies associated with hummingbird pollinators may promote mechanisms that lead to the generation of new species. By contrast, speciation rates associated with shifts in pollination syndromes (i.e. between insect and hummingbird pollination syndromes) were between 20 and 80 times lower than those within pollination syndromes (electronic supplementary material, table S5). These results indicate that the classical pollinator-shift hypothesis driving plant speciation does not explain the Gesnerioideae diversification. Instead, species richness in this plant group has been driven by speciation within hummingbird-pollinated lineages, without involving shifts among functional groups of pollinators. Why hummingbird pollination promotes plant speciation remains unclear, but we propose plausible mechanisms as a starting point for future research. First, the evolution of tubular or gullet-like flowers characterizing most hummingbird flowers may have directly accelerated speciation by promoting specialized relationships with the different categories of hummingbird species and the evolution of rapid prezygotic reproductive barriers [61,62]. Hummingbirds vary dramatically in bill size and shape, and these characteristics largely match with the flower morphology of the species they feed on [63–65]. Although most plant species are visited by several hummingbird species, specialization in the plant–pollinator network was found to increase at low and medium elevations [66], and in species-rich communities in which closely related hummingbirds visited distinct sets of flowering species [67]. Thus, floral specialization in hummingbird-pollinated species is frequently more specialized than initially assumed, which may result in greater potential for pollen segregation and, eventually, speciation. Second, flower specialization and specific pollen placement on the hummingbird body may limit interspecific pollen transfer among species sharing pollinators, thereby increasing the number of plant species that can co-occur in the same community [37,68,69]. It has been suggested that this process could decrease extinction rates [9], but also potentially increase the carrying capacity of hummingbird-pollinated lineages per unit of area, a factor that can limit the decline of diversification rates over time [70]. Finally, beyond floral specialization, bird pollination could have a direct impact on gene flow and the geography of plant speciation. Compared with insects, bird pollination increases the efficiency of pollen transfer and deposition [71], potentially affecting the connectivity among natural populations even in the context of a patchy distribution of suitable habitat [72]. Such improved transfer of pollen over long distances could enable the maintenance of isolated population providing the precondition for allopatric speciation. This has been suggested for the Andean species of Passiflora [73] and other hummingbird-pollinated lineages in the Gesneriaceae like Columnea and Dircaea [24,74]. Finally, and consistent with the last argument, hummingbird pollination is considered to be more efficient than insect pollination in Neotropical cloud forests at middle to high elevations, because insects are less active in cool, foggy, and wet conditions [75,76]. Altogether these patterns suggest that hummingbird-pollinated species could have more opportunities to colonize a wider geographical and climatic range, and to establish complex plant–pollinator interactions compared with insect-pollinated lineages. Although these hypotheses remain largely untested, such factors could potentially trigger diversification in angiosperm. Testing these hypotheses requires more complete morphological characterization of the plant species and plant–pollinator interaction data, to better understand how biotic interactions have shaped biodiversity and macro-evolutionary patterns in the Neotropical region. Our study provides new insights into one of the most intriguing factors influencing the diversification of Neotropical plant lineages, namely the impact of hummingbird pollination. However, additional ecological factors should not be excluded to conduct holistic examinations as encouraged by recent studies [15,16]. The evolution of epiphytism and different growth forms, the colonization of new biomes, and the geological history have potentially influenced the diversification and distribution of Gesnerioideae [17]. Currently, factors are usually evaluated independently, making it difficult to test their joint effects [77], and they should be taken with precaution, for instance, in traits with very few evolutionary transitions (such as epiphytism in the species-rich Columneinae clade). An optimal methodology should consider multiple factors simultaneously and allow a particular combination of those (i.e. colonization of a new area, with an in situ new trait state) to affect the diversification process (as discussed in [77,78]). Such approaches are, however, not yet available, and there is a current need to develop thoughtful tests for modelling simultaneously the success of multiple ecological interactions, considering the caveats of the methods and data, while integrating global biodiversity patterns. We identified a strong and positive effect of hummingbird pollination syndrome on the process of species diversification in the subfamily Gesnerioideae. This effect was probably triggered by the repeated acquisition of hummingbird pollination when this pollination niche became available in South America during the Early Miocene. Plants with a hummingbird pollination syndrome have a twofold increase in the rate of speciation, suggesting a positive effect of hummingbird pollination on the establishment of reproductive isolation. Our findings complement the global understanding of the diversification processes leading to the exceptional diversity of flowering plants in the Neotropics, and provide new directions towards further testing the role played by plant–pollinator relationships in the build-up of plant diversity. The data and additional details of the approaches used in this paper are available in the electronic supplementary material and at Dryad Digital Repository [79]. R script available from www2.unil.ch/phylo/files/software/make.simmap.BiSSE.R. M.L.S.-S., N.S. and M.P. planned and designed the research and performed the analyses. M.P. and J.L.C. conducted fieldwork and gathered data. M.L.S.-S., J.R., J.L.C., N.S. and M.P. wrote the manuscript. All authors gave final approval for publication. The authors have declared that no competing interests exist. This study was funded by the Ville de Genève, the Faculty of Biology and Medicine at the University of Lausanne and the Swiss National Science foundation (grant no. CRSII3_147630). We thank Alain Chautems, Lianka Cairampoma, Andrea O. De Araujo, Valquiria Ferreira Dutra, Mauro Peixoto, Angela Cano, Harri Lorenzi, Gabriel E. Ferreira, Juvenal E. Batista, Mireya Correa, Marina Wolowski and the staff of the Smithsonian Tropical Research Institute for their contribution to the fieldwork; Régine Niba and Fadil Avdija for their contribution to the laboratory work; Yvonne Menneret, Marina Magnette, Matthieu Grillet, Alexandre Chappuis, Bertrand Guigon and Vincent Goldschmid for the propagation and maintenance of the Gesneriaceae collection at Geneva. Collection permits were granted by the CNPq in Brazil (CMC 038/03) and in the ANAM in Panama (SC/P-43–10). Images in figure 1a were kindly provided by Ivonne SanMartin-Gajardo, Anton Weber and Leandro Freitas. We thank Daniele Silvestro for the helpful discussions and support for the analyses, and Catherine Graham and two anonymous reviewers for their advice on the manuscript. All analyses were run at the High-performance Computing Center (Vital-IT) from the Swiss Institute of Bioinformatics. FootnotesElectronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3718624. References
Page 8Mutation rates, like most organismal phenotypes, are subject to natural selection [1]. Alleles that increase the genomic mutation rate can not only be subject to direct selection but also to indirect selection on the fitness effects of mutations at other loci. A rich body of evolutionary theory predicts that in asexual populations, mutator alleles with little or no direct fitness effect can rise to fixation when they are linked to rare beneficial mutations. The probability of mutator hitchhiking is related to the supply and magnitude of such mutations, as well as the population size [2–6]. Experimental studies with Escherichia coli [7–12] and Saccharomyces cerevisiae [13–15] have provided evidence for this phenomenon, and clinical isolates of asexual pathogenic microbes have been shown to contain elevated frequencies of mutator strains [16–20], suggesting that mutator alleles have a propensity to rise in frequency during invasion of a new ecological niche. In contrast with asexual populations, mutator hitchhiking is not predicted to occur in sexual organisms. According to simulations and theoretical models, sex and recombination will erode linkage, separating beneficial mutations from the mutator alleles that caused them and preventing an increase in frequency of the mutators [5,21–25]. This prediction has garnered experimental support in S. cerevisiae [26]. However, some population genetic models have shown that under certain restrictive circumstances, fluctuating environments could select for mutator alleles in sexual populations [27–29]. Sexual microbes in the real world are mostly eukaryotes with facultatively sexual life cycles that undergo infrequent sex between periods of clonal growth [30]. The degree of linkage within the genome, as well as the strength and direction of selection and environmental change, are usually unknown. Thus, it is unclear whether natural populations of sexual microbes are expected to contain considerable frequencies of mutator strains, and to our knowledge none have been reported. The budding yeast S. cerevisiae has been collected from locations and ecological settings around the world [31–34]. Its life cycle can include asexual growth in haploid and diploid forms, as well as both outcrossing and extreme inbreeding (mating-type switching and intra-ascus mating) [35]. Outcrossing rates are inferred to be extremely low in natural populations, with most growth occurring clonally in the diploid phase [36,37]. Saccharomyces cerevisiae is also known to be an opportunistic pathogen [38,39]. Pathogenic isolates have been shown to have high levels of heterozygosity [40], suggesting that sexual outcrossing events prior to colonization of a new ecological niche may be associated with opportunistic pathogenicity. The global population of S. cerevisiae contains two alleles in the genes MLH1 (G761D) and PMS1 (R818 K) [41] that, when engineered together in a laboratory strain background, increase the mutation rate 20–400× above wild-type [42]. Mlh1p and Pms1p form a heterodimer that plays a major role in the process of DNA mismatch repair, a system for identifying and correcting mutations that occur during replication and recombination [43]. The mutation rate-increasing substitutions are found within the domains that mediate the Mlh1p–Pms1p interaction [44], although the degree of any disruption in the interaction is unclear. We will refer to the MLH1 and PMS1 alleles described above as the ‘mutator’ or ‘incompatible’ allele combination, though we show in this study that these alleles do not always cause a strong mutator phenotype when found together. Until recently, no natural isolates had been found carrying both alleles, in line with theoretical expectations of the absence of mutators in sexual populations. However, strains recently sequenced by Strope et al. [34] include a cluster of clinical isolates that carry both mutator alleles. This cluster of strains could represent a rare case of a facultatively sexual microbe with a naturally occurring mutator phenotype, and one that may be associated with invasion into a new ecological niche (the human body). Here we report direct estimates of genomic mutation rates in the clinical background containing the two incompatible alleles. Furthermore, we investigate the phylogenetic history of the mutator alleles and examine levels of genomic variation in both mutator and non-mutator genetic backgrounds. We find that these isolates have a mildly elevated mutation rate (approximately four times greater than closely related non-mutator strains), significantly lower than previously reported for incompatible combinations of MLH1 and PMS1 alleles assayed in a single laboratory genetic background [42]. Genomic analysis of mutator and non-mutator strains provided no evidence of a historically elevated mutation rate. We conclude that the mutational effect of the incompatible mismatch repair alleles is dampened by background genetic modifiers, and discuss possible explanations for our observations. Fifteen strains were chosen from the 100-genomes yeast panel [34]: four that contained both of the alleles necessary to confer a mutator phenotype (MLH1D761 and PMS1K818), six that contained only the PMS1K818 allele and five that contained only the MLH1D761 allele. Two MLH1D761-PMS1K818 strains that were reported to have a mutation rate approximately 100× wild-type [41], EAY1370 and EAY1363 (S288c and SK1 genetic backgrounds, respectively; generously provided by Eric Alani), were used as a control. All strains assayed were heterothallic haploid MATa (i.e. derived from the original diploid background), as the mutator alleles and the reporter for the mutation rate assay are recessive. In order for the control strains to be used in the fluctuation assays, the URA3 gene was restored via lithium acetate transformation [45] using strain YJM128 as the PCR template. The URA3 protein sequence in all strains used in this study was identical. Strains were Sanger sequenced at the PMS1 and MLH1 loci with the following primers to verify allele status: MLH1for-GCAGGTGAGATCATAATATCCC; MLH1rev-GGGCATACACTTTCAAATGAAACAC; PMS1for-CAGATAAACGATATAGATGTTCATCG; PMS1rev-CCTTCGAAAATGAGCTCCAATCA. Strains are listed in table 1.
Strains were grown in YPD [46]; spontaneous ura-mutants were selected on synthetic complete agar plates [46] supplemented with 1 mg ml−1 5-fluoro-orotic acid (5FOA) and 60 mg ml−1 uracil. Cultures were also treated with deflocculation buffer (20 mM sodium citrate, 5 mM EDTA), as strains YJM1639, YJM1851, YJM1747 and YJM1707 flocculated. The modified Jones protocol was used to estimate mutation rates [47]; following the procedure of Raynes et al. [26], spontaneous ura- mutants were assayed. All strains were grown overnight in YPD from freezer stocks; a new tube of 10 ml of YPD was inoculated with approximately 1000 cells and grown for 48 h. Five replicate cultures of 30 ml YPD were then inoculated with approximately 100–500 cells and grown for 48 h. The cultures were centrifuged, resuspended in 5 ml deflocculation buffer, washed and resuspended in 5 ml or 10 ml of water. For each replicate, 100 ul of an appropriate dilution was plated on YPD to estimate the population size and 300 ul was plated on 5FOA plates to select for spontaneous ura- mutants; for control strains, 10 ul + 90 ul water was also plated on 5FOA plates. Mutation rates were estimated using Mutation Rate Calculator software (provided by P. D. Sniegowski and P. J. Gerrish). This assay was replicated four times, with at least one control mutator strain included in two of the replicates. R v. 3.2.2 [48] and the lme4 package [49] were used to implement mixed models to test for differences in mutation rate. All models used the logarithm of mutation rate as the response variable. Assay was modelled as a random effect. Likelihood ratio tests were used to test for an assay effect in models containing a random strain effect and fixed effect for mutator alleles. To test for differences between groups of strains (engineered mutators, natural mutators and non-mutators), a random strain effect and a fixed effect for group differences were included. Coefficient estimates from these linear models quantified the fold differences between groups. For each of the 100 yeast genomes studied by Strope et al. [34] (except strain M22, which has a large amount of incomplete sequence), at least 10 kb of DNA sequence was extracted surrounding the mutator alleles, amino acid 761 of MLH1 and 818 of PMS1. We used a method developed by Gabriel et al. [50] and implemented in the program Haploview [51] to partition the sequences into blocks that did not show strong evidence of historical recombination. The sequence block containing the site of the mutator allele substitution in each gene was examined, and haplotype networks were constructed using PopART (http://popart.otago.ac.nz) with the median joining network option [52]. Complete genomes of the 15 strains whose mutation rates were measured were obtained from Strope et al. [34]. The genome of S. paradoxus was obtained from Scannell et al. [53]. These 16 complete genomes were aligned with mugsy v. 1r2.2 [54] using default options. MAF files were manipulated using maftools v. 0.1 [55]. To clarify comparisons, only alignment blocks where sequences from all 16 taxa aligned were considered. The majority of sequence fell into this subset of the full alignment—for the 12.16 Mb S288c genome, 11.75 Mb (97%) was present in the complete alignment and 11.24 Mb (92%) was present in alignment blocks containing sequence from all 16 taxa. In order to obtain a reliable set of polymorphic sites, the following criteria were used to filter further: (i) ends of alignment blocks were trimmed if the sequence of any species consisted of only gaps; (ii) only biallelic single nucleotide variants were considered, not insertions–deletions (indels) or multi-allelic sites; and (iii) for a particular site, any strains with a gap or N at that site were ignored. Strope et al. [34] sequenced the genomes of 93 S. cerevisiae strains derived from isolates collected globally with a particular focus on clinical isolates. These strains were haploid or homozygous diploid segregants of the isolates themselves. Surprisingly, four of the sequenced strains carried incompatible alleles in the mismatch repair genes MLH1 (G761D) and PMS1 (R818 K) that are thought to confer a dramatic increase in mutation rate when they co-occur [41]. While both alleles are known to be segregating in natural populations (24% and 12%, for MLH1D761 and PMS1K818, respectively), it was previously postulated that the alleles would not be found together in nature because the combination would be selected against due to the accumulation of deleterious mutations [41]. Indeed, a previous investigation into a diverse panel of yeast isolates uncovered no strains with the combination [42]. While this manuscript was in review, Bui et al. [56] published an analysis of MLH1 and PMS1 across an even larger panel of 1010 natural isolates and found only 19 isolates carrying the alleles together in either heterozygous or homozygous form, confirming that it is unusual to find this ‘mutator’ allele combination. The four strains recently discovered by Strope et al. [34] to be carrying the mutator allele combination are all clinical isolates classified as admixed ‘mosaics’, which do not show pure ancestry from a single S. cerevisiae population [34]. The isolates from which these strains derived were collected from geographically disparate locations in California (n = 3) and North Carolina (n = 1). To understand the history of mutator alleles at MLH1 and PMS1, we analysed the haplotype structure of regions surrounding these loci in the 100-genomes yeast panel [34]. We partitioned this sequence into blocks showing little evidence of historical recombination, and found that each derived allele appears to have arisen once (on a single haplotype; figure 1) and spread widely within S. cerevisiae, as evidenced by its presence among diverse wine/European and mosaic strains. This observation is expected due to the small or non-existent fitness effect of each mutator allele when found with a wild-type allele at the opposing locus [41]. Figure 1. Haplotype networks constructed using variation surrounding mutator alleles. (a) MLH1 and (b) PMS1 networks. In both networks, nodes indicate individual haplotypes and are sized proportional to the number of strains carrying that haplotype. Lines connect closely related haplotypes, with notches indicating the number of mutations that separate adjacent haplotypes. Boxed strain names (or boxed numbers near large nodes) indicate strains assayed in this study. White circles indicate haplotypes carrying the derived allele (the mutator combination consists of the derived allele at both loci) and grey circles indicate haplotypes carrying the ancestral allele. Both networks are derived from sequence data of 99 strains, including all strains assayed in this study (Material and methods). Strain labels indicate diploid parent name (to facilitate comparison to genomes sequenced in [34]). (Online version in colour.)
We assayed mutation rate in the four strains carrying the two incompatible mismatch repair alleles, 11 related strains carrying compatible combinations of alleles, and two laboratory strains engineered to carry the incompatible alleles and previously shown to exhibit high mutation rate (control strains). We refer to these strains as natural mutator strains, non-mutator strains and engineered mutator strains, respectively. Although the strains we studied were derived from diploid isolates that could be heterozygous for the incompatible alleles, and therefore not exhibit a mutator phenotype, these derived strains have direct relevance to S. cerevisiae population biology: a single round of sporulation followed by intra-ascus mating or mother–daughter mating, the most common forms of mating in Saccharomyces yeasts [36], would produce a diploid yeast cell homozygous for both mutator alleles. At least one of the natural clinical isolates, YJM523, the diploid isolate from which YJM555 (isogenic to YJM1785) was derived, was shown to be homozygous for both mutator alleles [56], and should therefore exhibit the mutator phenotype in the clinical setting from which it was isolated. We measured mutation rate multiple times (median 4, range 3–5) for each strain, except the two previously assayed engineered mutator strains (n = 1 or 2) and a single non-mutator strain (YJM1876; n = 1). Given the large number of mutation rate measurements conducted in this study, we carried out the mutation rate assay in multiple batches and attempted to include all strains in each batch. Experimental variation led to significant variation in mutation rate estimates between batches ( Figure 2. Mutation rates of strains carrying mutator and non-mutator allele combinations. (a) For each fluctuation assay, clinical strains with only one of the mutator alleles were averaged to provide the baseline mutation rate; the mutation rate for each strain was divided by this baseline. Shapes (see legend) show data from individual assays and diamonds indicate average for the strain over all assays. (i) Strains sequenced in [34] and assayed in this study. (ii) Comparison of mutation rates of engineered mutator strains from [41] to 100 genomes natural mutator strains. Note the difference in y-axis scale between the panels. (b) Estimates of mutation rate by strain. Shapes (see legend) show measurements from individual assays. Boxes show 25%, 50% and 75% quantiles. In both panels, strain labels indicate either diploid parent name (to facilitate comparison to genomes sequenced in [34]) or genetic background (for engineered mutator strains in [41]). See table 1, for more details on each strain. (Online version in colour.)
Although the mutator allele combination did not appear to give rise to substantially higher mutation rates in the natural mutator strains, this may not reflect the historical influence of this pair of alleles. In particular, an initial elevation of genomic mutation rate when the mutator alleles were first combined on a single genetic background could be dampened by the subsequent appearance of a modifier of mutation rate [25]. This scenario could manifest as an increase in genome-wide mutations in only those strains carrying the mutator combination (due to the historically higher mutation rate in those backgrounds) despite the present similar overall genomic mutation rates for both mutators and non-mutators. To search for evidence of this pattern, we examined the number of derived mutations in each strain. We first aligned full genomes of the assayed strains to the genome of S. paradoxus, the closest extant relative of S. cerevisiae. Next, we examined sites that were polymorphic in S. cerevisiae and for which we could determine the derived allele using the S. paradoxus sequence. We found that the number of derived alleles in natural mutator strains was not higher than the number in non-mutator strains (natural mutator mean 44 892; non-mutator mean 45 793). The same pattern held when we restricted our focus to the natural mutator strains and their two closest relatives at the whole-genome level (YJM1775 and YJM1787) in order to focus only on mutations that occurred along branches separating the natural mutators from their close relatives (natural mutator mean 30 218; non-mutator mean 30 427). Nevertheless, it is difficult to falsify the hypothesis that background modifiers of mutation rate arose after the incompatible mutator alleles were brought together in a single genetic background, for at least two reasons. First, the mutation rate modifier allele(s) could have arisen shortly after the mutator alleles were brought together, resulting in a very small excess of mutations due to a short period of elevated mutation rate. Second, the complex population history of diverse yeast isolates, along with an unknown natural life cycle (i.e. clonal growth versus sexual reproduction), could obscure historical patterns of sequence divergence within and among lineages. We have shown that the mutation rate in natural mutator strains is slightly elevated (approximately 4×) compared with their closest non-mutator relatives. Why, then, do we not observe any enrichment of derived mutations in natural mutator strains? There are two possibilities. First, we have few strains to compare that are closely related to the natural mutator strains. A larger number of evolutionary lineages could provide power to detect subtle trends in patterns of mutation over time. Second, the mutator alleles may have been present in heterozygous form, and thus not exhibiting the mutator phenotype, for a significant portion of their evolutionary history. The lack of a strong effect of the mutator alleles on mutation rate in the clinical isolate backgrounds that we tested suggests a role for epistatic interactions with other loci. Previous work investigating this incompatibility utilized natural mutator alleles that contained additional single nucleotide polymorphism (SNP) variants within MLH1 and PMS1 [42]. These results demonstrated that intragenic modifiers could modulate the mutation rate over a 20-fold range, although in all cases the mutation rate was still significantly elevated (at least 24× [42]). By contrast, our results suggest the presence of modifiers that greatly suppress the elevated mutation rate phenotype. Indeed, the slight elevation in mutation rate that we observed among natural mutator strains matches roughly with variation in mutation rate among different backgrounds with compatible allele combinations (up to sixfold increase relative to S288c) [42]. We used sequence data [34] to search for possible intragenic modifiers in the set of strains we assayed. Demogines et al. [42] found two alleles, PMS1F165C and MLH1L271P, that modified mutation rate in strains with the incompatible mutator allele combination. The PMS1F165C allele that increases mutation rate was not found among the strains we examined. However, the MLH1L271P allele that decreases mutation rate threefold to fourfold was polymorphic among the natural mutator strains. Strain YJM1785 (isogenic to YJM555) carries the reference T allele at position 595 697 on chromosome XIII, while the other three natural mutators carry the C allele that dampens the increased mutation rate due to the incompatible alleles. We observed a modestly lower mutation rate (approximately twofold) among the natural mutators carrying this ‘protective’ allele, in line with previous observations [42]. In addition to these variants, there were 18 SNP and indel variants in MLH1 and 40 in PMS1 that were polymorphic in the full panel of 100 Genomes strains we studied, of which 11 and 13 result in non-synonymous substitutions, respectively. Of these, only two SNPs in MLH1 and one SNP in PMS1 were polymorphic among the natural mutator strains. However, all had the same allelic pattern as the MLH1L271P allele (YJM1785 carried one allele and the other natural mutators carried the other), preventing us from determining whether any of these additional non-synonymous polymorphisms also modulate mutation rate in the incompatible background. We expanded our focus to the full genome to search for candidate loci that could contain mutation-rate-lowering modifiers of the incompatible allele combination. First, we focused on a cluster of mutator and non-mutator strains (YJM1775, YJM1751 and YJM1753—isogenic to YJM320, YJM541 and YJM554, respectively) that are among the most closely related among 100 genomes strains. We compared YJM1775 along with S288c and Sigma1278b (non-mutators that may not contain compensatory mutations) versus YJM1751 and YJM1753 (natural mutators that presumably do contain compensatory mutations). However, there were no polymorphic sites that had one allele shared by the mutators and another by the non-mutators. Next, we expanded our focus to search for sites where all natural mutators shared one allele but their two closest non-mutator relatives at the whole-genome level (YJM1775 and YJM1787—isogenic to YJM320 and YJM681, respectively) shared another allele. Unfortunately, we found over 1000 sites that fit these criteria. We used the STRING database [57] to filter these sites and found only one potential polymorphism of interest, in the MutS homologue MSH5. However, this polymorphism is synonymous, and a predicted interaction between MLH1 and MSH5 is probably due to their shared role in promoting crossing over [58,59], rather than DNA mismatch repair. Overall, our inability to detect candidates may be because these strains are sufficiently divergent that a few large-effect compensatory mutations are obscured by the large amount of variation at other loci, or because there are many small-effect modifiers of mutation rate scattered throughout the genome. Moreover, these modifiers may have already been present in the genomic background on which the mutator alleles were first brought together. The complex machinery that is required to replicate DNA with high fidelity [60] makes it possible for many loci to influence the mutation rate. Given the challenges of making precise measurements of the mutation rate phenotype, uncovering modifiers of small effect using QTL approaches is not feasible. In this study, we characterized the mutation rate in strains derived from natural isolates carrying a pair of incompatible alleles in mismatch repair genes that were previously shown to lead to a much higher genomic mutation rate when combined in one background [41]. The existence of naturally derived strains carrying this pair of alleles is unexpected given the prediction, from previous theoretical and experimental work, that natural populations of sexual microbes should not contain considerable frequencies of individuals carrying mutator alleles. We show that the mutation rate of laboratory strains engineered to carry the incompatible MLH1D761 and PMS1K818 mutator alleles is high, as demonstrated previously [41], but that these alleles only mildly elevate mutation rate in strains derived from natural clinical isolates. Thus, the effect of the incompatible alleles on mutation rate is modulated by genomic background. The observation that all natural mutator strains identified by Strope et al. [34] are derived from clinical isolates raises the possibility that the human body, an unusual environment for S. cerevisiae, poses unique selective pressures that could select for the mutator combination hitchhiking along with advantageous alleles during a period of clonal growth. Our data do not strongly support this contention given the relatively slight elevation of mutation rate in natural mutator strains. However, the haploid natural mutator strains we tested were derived from diploid clinical isolates that were probably highly heterozygous [40]. Thus, it remains a possibility, albeit unlikely, that within the population of S. cerevisiae present in a single patient, these alleles could be present at moderate frequencies along with other mutation rate modifiers such that rare offspring containing the mutator alleles and no suppressors could contribute adaptive mutations that are quickly detached from their high-mutation-rate background via recombination. Given the possible link between mutators and opportunistic pathogens, the question of the effect of mutator alleles in pathogenic strains of a facultatively sexual microbe is of particular interest. Our results demonstrate that a cluster of clinical yeast strains with putative mutator alleles are in fact not mutators. More broadly, it is intriguing that hundreds of isolates of the facultatively sexual yeast S. cerevisiae—from around the globe and from various ecological niches—have been sequenced, and the only isolates found to contain the incompatible mutator combination do not actually have a strongly elevated mutation rate. It is worth noting that mutator alleles are more likely to be maintained in facultatively sexual organisms than in obligately sexual organisms due to periods of clonal growth where such alleles cannot be decoupled from beneficial mutations. Overall, our results provide rare and compelling evidence that supports evolutionary theory suggesting that sexual organisms are unlikely to harbour alleles that increase their genomic mutation rate. Mutation rate data, genome sequence data and code to reproduce statistical and genomic analyses is available at http://dx.doi.org/10.5061/dryad.125p3 [61]. H.A.M. conceived of the study. B.M. and H.A.M. carried out assays to estimate mutation rate. D.A.S. and H.A.M. analysed the data. D.A.S. conducted bioinformatic analyses of genomic variation. H.A.M. and D.A.S. wrote the paper. P.M.M. provided strains and overall advice on the project. All authors gave final approval for publication. The authors have no competing interests. This work was supported by a College of William and Mary Faculty Research Grant (H.A.M.), NIH F32 GM110997 (D.A.S.) and NIH R01 GM098287 (P.M.M.). We thank John McCusker and Eric Alani for strains, and P. D. Sniegowski for helpful comments on the manuscript. FootnotesReferences
Page 9Infectious diseases have impacted human populations throughout history. While studies of contemporary diseases benefit from modern methods which allow rapid collection and dissemination of information about the diseases' effects on populations, studying diseases in the past is more challenging. Some historical disease events, such as the plague that struck Europe in the fourteenth century, are relatively well understood, but in other cases little or conflicting information is available about the diseases and their impacts. One such event is the series of epidemics caused by colonization of the Americas by Europeans. By some accounts, disease killed more than 90% of the native population and caused widespread social disruption, but other reports suggest smaller impacts (see varying estimates of population mortality rates in Dobyns [1] and Crosby [2]). Despite uncertainty about the disease-related effects of European colonization of the Americas, that interaction has been studied much more extensively than have the migrations into Africa, in particular southern Africa. In fact, Europeans were not the first immigrant group to settle in southern Africa. There is evidence for at least two earlier within-Africa population movements. One may be associated with the introduction of pastoralism to southern Africa around 2 000 years ago [3–7]. Admixture studies find a small fraction of east African pastoralist ancestry in southern African pastoralist Khoekhoe populations, indicating that the cultural practice of pastoralism may have been transported to southern Africa by a relatively small number of east African individuals who assimilated into the local populations [3,8]. Despite the low number of immigrants participating in this migration, the introduction of pastoralism likely resulted in a large increase in disease burden, particularly zoonotic disease. A later migration of Bantu-speaking farmers from west and central Africa (e.g. [9]) arrived in the south around 1 200 years ago [10]. This migration was a larger-scale movement of people and resulted in the many Bantu-speaking groups found in southern Africa today (e.g. [8,9,11]). This later migration that involved many individuals and the new cultural practice of farming is likely to be a better candidate for a large-scale effect on disease burden. While many of the aboriginal southern African San populations remained mobile hunter–gatherers, they may have been exposed to diseases associated with sedentary or herding lifestyles through interactions with immigrant groups and local groups that adopted those modes of subsistence. More recently, European colonists began arriving in southern Africa around 1650. They first settled close to the southern coast where they primarily came into contact with indigenous groups living close to the African south coast (Khoekhoe herders and Tuu-speaking San groups, most probably ancestral to, e.g. ‡Khomani and Karretjie groups). This interaction resulted in disease epidemics, including several documented smallpox epidemics in the 1700s that killed up to 90% of the Cape Khoekhoe groups [12]. However, while some effects of diseases introduced during European colonization are better understood than those due to earlier within-Africa movements, questions remain even about that period. Studies of both the within-Africa and European migrations are hindered by the lack of pre-arrival indigenous population records, which makes it impossible to estimate the impacts of introduced diseases using traditional measures of mortality and morbidity. Genetic analyses, however, offer another account of population history and enable us to find signatures of past events that we could not otherwise measure. Episodes of natural selection, as would occur during epidemics of introduced infectious diseases, are expected to leave signatures in the genome [13] such as extended lengths of haplotype homozygosity (measured with long-range haplotype scores such as iHS [14] and XP-EHH [15]), differentiation between populations (e.g. FST), and change along a specific lineage in a three-way population differentiation comparison (the population branch statistic, PBS [16]). These different statistics would capture signals of natural selection at different time points, from very recent (perhaps a couple of hundred years) to far back in time (beyond human emergence) [13]. We compared genome-wide population-genetic data for signals of disease-related selection in the ‡Khomani and the Ju|’hoansi peoples of southern Africa. The indigenous peoples of southern Africa are the San people (hunter–gatherers) and the closely related Khoekhoe people (pastoralists) who belong to a common branch of the human lineage that diverged more than 100 000 years ago from all other modern humans, thus representing the earliest diversification event among modern humans [8,17,18]. The ‡Khomani is a San group that historically resided in the southern Kalahari region of southern Africa while the Ju|’hoansi, another San group, historically resided in the northwestern part of southern Africa (figure 1). These two populations are estimated to share common ancestry 35 000 years ago [8]. Their different geographical locations have resulted in disparate levels of contact with outside groups entering southern Africa within the last 2 000 years. The Ju|’hoansi population has been isolated throughout its history and has had low levels of contact and gene flow with outside groups, whereas the ‡Khomani population has experienced much more contact and gene flow with both immigrants practising farming and the local indigenous groups that adopted pastoralism [3,8,19,20]. Figure 1. Sampling locations of the Ju|’hoansi (red), a population with a history of isolation, and the ‡Khomani (blue), a population with abundant contact with Khoekhoe pastoralists, Bantu-speaking farmers, and European colonists. (Online version in colour.)
We used two methods to compare signals of selection in these populations. Since any statistic designed to detect selection will also—at least to some extent—be affected by demography (such as bottlenecks, admixture, and expansions), we first explicitly contrasted summary statistic values of genetic tests for selection in immune genes versus all genes in the two populations. As demography affects genetic variation at all (autosomal) genes in a population equally, this approach allowed us to control for possible biases due to demography. Next, we examined regions of the genome that were in the top fraction for three selection statistics (iHS, PBS, and FST) for enrichment of immune genes. We combined three summary statistics to reduce the number of false positives [21]. This dual approach allowed us to examine signals of selection on immune genes throughout the genome as well as in regions with the strongest indicators of selection, resulting in a fuller understanding of infectious disease-related selection than using a single method. We first tested whether there was a difference in selection signals in the full set of immune SNPs compared to the full set of (non-immune) genic SNPs (using a Student's t-test and/or a Mann–Whitney U-test), and whether such a difference varied between the Ju|’hoansi and the ‡Khomani. Using a weighted block jackknife approach, we examined to what extent a signal among SNPs was due to a few genomic regions with tightly linked SNPs (see below). The full list of genes (80 922 genes) overlapped with 642 560 SNPs in our dataset. The immune gene list (855 genes) contained 33 578 SNPs (5.2% of the full list). In both the Ju|’hoansi and the ‡Khomani, iHS values were significantly greater for immune genes compared with all genes (figure 2, Ju|’hoansi p = 0.026; ‡Khomani p =5.1 × 10−6, Mann–Whitney U-test). Figure 2. iHS results (a) quantile-quantile (qq) plot for SNPs in immune genes (x-axis) versus SNPs in all genes (y-axis) in the Ju|'hoansi (mean |iHS| = 0.791 among 13 126 SNPs in immune genes and mean |iHS| = 0.778 among 415 266 SNPs in other genes, p = 0.0263 based on Mann–Whitney U-test), (b) distribution of |iHS| at SNPs in immune genes (red) and at SNPs in all genes (black) in the Ju|'hoansi, (c) qq plot for SNPs in immune genes (x-axis) versus SNPs in all genes (y-axis) in the ‡Khomani (mean |iHS| = 0.805 among 14 422 SNPs in immune genes and mean |iHS| = 0.778 among 458 852 SNPs in other genes, p = 5.08 × 10−6 based on Mann–Whitney U-test), (d) distribution of |iHS| at SNPs in immune genes (red) and at SNPs in all genes (black) in the ‡Khomani. (Online version in colour.)
To test whether frequency changes were different between SNPs in immune genes and SNPs in other genes in the two populations, we used the PBS statistic that produces a three-way population topology proportional to differentiation among groups (we used the Herero, a Bantu-speaking group, as an outgroup). The length of the Ju|’hoansi branch based on SNPs in immune genes is slightly shorter (mean = 0.029) than the branch based on SNPs in all genes (mean = 0.030, Mann–Whitney U-test p-value = 0.010). By contrast, the ‡Khomani branch based on SNPs in immune genes is much longer (mean = 0.015) than the branch length based on SNPs in all genes (mean = 0.00017, Mann–Whitney U-test p-value = 0.0023) indicating stronger selection in immune genes than all genes in the ‡Khomani (figure 3a,b,e,f). Figure 3. PBS and BLD analysis for Ju|'hoansi versus ‡Khomani (p-values based on Mann–Whitney U-test) (a) qq plot of PBS values for SNPs in immune genes (n = 19 737, x-axis) and SNPs in all genes (n = 627 769, y-axis) in Ju|'hoansi (mean PBS = 0.0289 among SNPs in immune genes and mean PBS = 0.0300 among SNPs in other genes, p = 0.00997), (b) qq plot of PBS values for SNPs in immune genes (x-axis) and SNPs in all genes (y-axis) in ‡Khomani (mean PBS = 0.00147 among SNPs in immune genes and mean PBS = 0.00017 among SNPs in other genes, p = 0.0023), (c) qq-plot for BLD values of SNPs in immune genes (x-axis) and SNPs in all genes (y-axis) in Ju|'hoansi versus ‡Khomani (mean BLD = 0.0274 among SNPs in immune genes and mean BLD = 0.0298 among SNPs in other genes, p = 0.00527), (d) distribution of BLD at SNPs in immune genes (red) and at SNPs in all genes (black), (e) tree representation of PBS values for SNPs in immune genes, and (f) tree representation of PBS values for SNPs in other genes. (Online version in colour.)
In addition to testing whether values are different for SNPs in immune genes and SNPs in other genes for each population separately, we tested whether the distribution of the difference between Ju|’hoansi and ‡Khomani was different for SNPs in immune genes and for SNPs in other genes. By contrasting FST-based population branch lengths for different categories of SNPs, we computed the difference between PBS(Ju|’hoansi) and PBS(‡Khomani) (i.e. Branch Length Difference, BLD) for SNPs in immune genes and SNPs in other genes. These two distributions were then compared. BLD was significantly smaller (Mann–Whitney U-test p = 0.0053) among immune SNPs (mean = 0.027) than among genic SNPs (mean = 0.030) (figure 3c,d). BLD being smaller among immune SNPs is consistent with a relatively longer ‡Khomani branch at immune SNPs than at other genic SNPs. That both means are positive indicates a longer Ju|’hoansi branch, probably an effect of a larger proportion of Bantu-speaking admixture (Bantu-speaking Herero, used as an outgroup) in ‡Khomani than in Ju|’hoansi, consistent with previous observations [8]. Using the East African Maasai as an outgroup instead of the Herero in the PBS analysis led to very similar results (data not shown). The XP-EHH statistic captures differences between pairs of populations in extended haplotype homozygosity that signals local adaptation. XP-EHH values were significantly (t-test p = 4.7 × 10−9) more negative at SNPs located within immune genes (mean = −0.052) compared with SNPs in all genes (mean = −0.012; figure 4). Negative XP-EHH values for the Ju|’hoansi and ‡Khomani comparison corresponds to longer haplotypes in ‡Khomani relative to haplotypes in Ju|’hoansi, consistent with relatively stronger selection on immune genes in ‡Khomani than in Ju|’hoansi. To verify that this result is not due to bias caused by comparing a small set of genes to a much larger set of genes, we replaced immune SNPs with a random set of SNPs from the full set of genic SNPs and then performed the same XP-EHH analysis contrasting Ju|’hoansi and ‡Khomani. We repeated these 100 times and only four of these were significant at p < 0.05 (close to the expected five out of 100 under the null model of no difference) with a minimum p-value of 0.0029 (compared to the observed p-value of 4.7 × 10−9). Figure 4. XP-EHH analysis for Ju|'hoansi versus ‡Khomani (a) qq-plot for SNPs in immune genes (x-axis) versus SNPs in all genes (y-axis); (mean XP-EHH = −0.0523 among 20 662 SNPs in immune genes and mean XP-EHH = −0.0122 among 648 000 SNPs in other genes, p = 4.7 × 10−9 based on Student's t-test) and (b) distribution at SNPs in immune genes (red) and at SNPs in all genes (black). (Online version in colour.)
To test the effect of populations used for our comparisons, we also performed the iHS, XP-EHH, PBS, and BLD analyses using the Karretjie [8], another San population exposed to migrant populations, in place of the ‡Khomani. Results from this comparison corroborate the results presented here: the Ju|’hoansi appear to be less affected by selection on immune genes than the Karretjie (electronic supplementary material, figures S4–S6). It is important to point out that we assess significance using statistical tests that assume independence among SNPs. A weighted block jackknife analysis to study the effect of linkage shows that the only statistic that has non-overlapping 95% CIs between immune and genic SNPs is |iHS| in ‡Khomani (see the electronic supplementary material, figure S1), which points to a very distinct difference between signals of recent selection on immune and genic SNPs. It also suggests that there are specific regions driving the signal. To investigate this closer, we identified blocks (of 5 Mbp) that were driving the difference between the immune and genic SNPs for each statistic (see Material and methods). The results are shown in the electronic supplementary material, table S1. In addition to examining differences in summary values of selection statistics at immune and all genes, we also identified genomic regions with the strongest signals of selection, as indicated by high values of three separate selection statistics, and evaluated these regions for enrichment of immune genes. We first selected the 10 most significant iHS windows in each population and compared these windows with the same genomic location in the other population. We visually examined the 20 pairs of windows to determine whether there was evidence of selection also in the other population. Of the top 10 iHS windows in the Ju|’hoansi, six showed evidence of selection in the ‡Khomani as well, indicating that most of the genome-regions with the strong selection affect both groups. Of these four windows with Ju|’hoansi-specific evidence of selection, two contained no genes and the other two contained no immune genes. Because we were searching for regions with high values for all three summary statistics that contained immune genes, and no iHS windows in the Ju|’hoansi passed the iHS selection step, none were considered further in our analysis (table 1). By contrast, nine of the top 10 iHS windows in the ‡Khomani showed evidence of selection unique to that population. Of these nine windows, seven contained genes with immune function and were considered for FST analysis.
We used a 99th percentile FST cut-off of 0.238 (mean FST = 0.0184) to determine extremely differentiated Ju|’hoansi–‡Khomani SNPs. Of the seven regions in the ‡Khomani selected during the iHS step, five contained SNPs with FST values above this cut-off. The number of significant SNPs ranged from one to 11 per window (21 in total). We then calculated PBS values for the 21 SNPs selected via the combined iHS and FST steps. Because all of the regions being examined were chosen due to selection in the ‡Khomani, we focused on SNPs for which the ‡Khomani (not the Ju|’hoansi or Herero) had the long branch (a cut-off of 2.5 times the second-longest branch). This resulted in four windows containing eight SNPs that showed strong evidence of selection and differentiation due to adaptation along the ‡Khomani lineage. Varying the iHS, FST, and PBS cut-offs did not qualitatively change our results (table 2).
Three of the four regions selected by this three-test process contain SNPs with distinct signs of selection along the ‡Khomani branch in close proximity to immune genes (table 3). SNPs in the fourth region were 600 kb away from the nearest immune gene (HSPD1). Owing to the distance, this region was excluded from further analyses. The eight immune genes in the three retained regions were two members of the Fc-receptor-like cluster (FCRL4 and FCRL5), located on chromosome 1 around 157.5 Mb (electronic supplementary material, figure S7a); the Butyrophilin family (BTN2A1, 2A2, 3A1, 3A2, and 3A3, located in the extended MHC on chromosome 6, around 26.4 Mb (electronic supplementary material, figure S7b); and PRSS16, also in the extended MHC region on chromosome 6 around 27.5 Mb (electronic supplementary material, figure S7c).
The indigenous populations of southern Africa experienced different levels of interactions and exposure to groups migrating into the region based on their historical locations. Their past demographic histories coupled with varying degrees of interactions with external groups and concomitant exposure to their unfamiliar diseases, may explain signatures of selection at immune genes. Here, we have used the Ju|’hoansi as a representative group for a population with minimal exposure to incoming farmer/herder cultures in the past 2 000 years and the ‡Khomani as a representative of an exposed population in the same time period. We find several lines of genetic evidence in accord with our hypothesis that the ‡Khomani underwent stronger selection on immune function than did the Ju|’hoansi. Using the framework based on extended lengths of haplotype homozygosity [14,15], we find that iHS values are significantly higher for immune genes than for all genes in both the Ju|’hoansi and ‡Khomani (figure 2). This result indicates that the immune system may have been a target of selection in both populations over a long period of their history in southern Africa. However, XP-EHH values for immune gene regions and genic regions show smaller values at SNPs in immune genes in the Ju|’hoansi and ‡Khomani (figure 4), indicating that while selection on the immune system may have occurred in both populations, it has likely to have had a stronger effect in the ‡Khomani. The divergent selection pressure between the populations is further demonstrated by PBS and BLD analyses. Branch lengths as estimated by PBS are significantly longer for immune genes than for all genes in the ‡Khomani, with the converse true in the Ju|’hoansi (figure 3). BLD values are significantly smaller at immune genes than at all genes suggesting historically stronger directional selection at immune genes in the ‡Khomani than in the Ju|’hoansi. While the direction of selection is not indicated by XP-EHH and BLD for the comparison of immune genes versus all genes (and it could in principle be explained by stronger directional selection at non-immune genes in Ju|’hoansi than in ‡Khomani), the iHS and PBS analyses show that it is immune genes that adapted faster. The large difference in power between the framework based on extended lengths of haplotype homozygosity (which had much smaller p-values) and the FST-based framework suggests that linkage disequilibrium patterns, rather than frequency differences, contain most of the information, perhaps indicating more recent selection (e.g. [15]). A more careful inspection of the XP-EHH and BLD distributions (figures 3 and 4) suggests that while the significant result for BLD is based on a few outliers, the XP-EHH result is due to a general left-skew of the XP-EHH values. That the selection tests give different results may also be due to the fact that they capture different aspects of selection. iHS is best at detecting recent selective sweeps that have not gone to fixation [14], while PBS and XP-EHH measure more ancient events [15,22]. To precisely determine the age of the selective events that gave rise to the genomic signals is difficult. The Ju|’hoansi and the ‡Khomani diverged around 35 000 years ago [8], which provides an upper limit for the time of selection. The signals could be a result of selection during any or all of the three major population movements we are aware of, or could be a result of earlier unknown migration events. Alternatively, the signals could be due to adaptation to disease exposure in general as a result of repeated exposure to immigrants and unfamiliar diseases. We note that although we contrast the statistics on a genome-wide scale, the qq-plots suggest that the signals are driven by relatively few regions. In fact, a more conservative weighted block jackknife analysis to identify the genomic regions driving the difference between immune and genic SNPs (electronic supplementary material, table S1) shows that the test statistics are generally driven by distinct regions of the genome. The MHC region, for example, drives |iHS| signals, but not XP-EHH, PBS, and BLD (electronic supplementary material, table S1). These regions are also distinct from those identified using stringent cut-offs for the three-test statistics. This effect may be caused by the fact that the test statistics are sensitive to selection events of different time frames. This observation also suggests that the difference between the Ju|’hoansi and ‡Khomani is not due to a single event, rather it is the result of a combination of events (such as the greater exposure of the ‡Khomani to both Bantu-speaking and European migrants or to other factors that may have introduced differences prior to their exposure to more recent immigrant populations). One potential concern is that not only are the ‡Khomani more likely to have suffered from a higher disease burden than the Ju|’hoansi, they also have a larger proportion of genomic material of Bantu-speaker ancestry. Since the farming Bantu-speaker populations likely also experienced an increased disease exposure as a consequence of their change in subsistence mode, it is possible that the difference in selection signals between San populations at immune genes compared to all genes merely reflects different levels of admixture. However, this scenario predicts a stronger difference between immune genes and all genes in the Bantu-speaking population compared to the San populations. We do not find that effect with either iHS (see the electronic supplementary material, figure S2) or PBS (electronic supplementary material, figure S3) for the Bantu-speaking population suggesting that such a scenario does not explain the difference in disease adaptation between the ‡Khomani and the Ju|’hoansi. Finally, pathogen load has been shown to be correlated to climate, specifically to precipitation and temperature [23]. At least for recent climate data, we could not detect any difference with respect to precipitation and temperature between the geographical areas of the Ju|’hoansi and the ‡Khomani (electronic supplementary material, table S2). However, to fully investigate this possibility, climatic data over long timescales is required, and such a test would also require a strong assumption of geographically stable populations. Our second approach, using extreme values of the selection statistics to locate regions of the genome with strong indications of selection unique to each population, also indicates that selection on immune genes has been a stronger force in the ‡Khomani than the Ju|’hoansi. While there were seven potentially immune SNPs in the ‡Khomani that were highly significant for all three tests, no SNPs in the Ju|’hoansi met these criteria. The eight genes located near significant SNPs have a range of roles in the immune system, some more well defined than others. PRSS16 on chromosome 6 encodes a thymus-specific serine protease involved in MHC class II antigen presentation to T cells during positive selection [24], and it shows a dramatic signal of selection in ‡Khomani (electronic supplementary material, figure S7c). Another selected region on chromosome 6 contains the butyrophilin (BTN) genes, including BTN2A1, 2A2, 2A3, 3A1, and 3A2 (electronic supplementary material, figure S7b). BTN family members are structurally similar to B7 co-stimulators and those whose function has been investigated are inhibitory co-stimulators with immunosuppressive function [25]. A third region with extreme signals of selection in the ‡Khomani contains two members of a family of Fc receptor-like genes, FCRL4 and FCRL5 (electronic supplementary material, figure S7a) on chromosome 1. These are B-cell membrane receptor proteins with both inhibitory and stimulatory signalling subunits [26]. Several of the genes identified as putative targets of selection have inhibitory functions, but further investigation will be required to know whether the variants selected in the ‡Khomani lead to up- or downregulation of these genes, i.e. whether selection favoured increased or decreased immune response. The type of immune response that would be beneficial depends on the diseases to which a population is exposed. While certain diseases are more efficiently fought with an increased immune response, others, including some influenza pandemics [27,28] and SARS [29], cause damage via over-activation of the immune system. The roles of these genes in response to specific infectious diseases are still unknown. One disease known to have affected indigenous southern African populations is the repeated epidemics of smallpox during European colonization [12], which had severe impact on the affected populations. Studies of the vaccinia virus, the closely related poxvirus from which the smallpox vaccine was derived [30], have shown that one way poxviruses evade the immune response is by blocking signalling pathways, particularly those activating the Toll-like receptor [31,32] and complement pathways [33,34], two important components of the innate immune system. Several of the immune system genes identified in this study are involved in signalling in the adaptive immune system, which could potentially compensate for the downregulation of the innate system caused by smallpox. A recent study identified two recently emerged alleles in the ‡Khomani that also affect signalling, indicating this may have been a common target of selection [35]. This study also implicates a very high diversity at the KIR2DL1 gene (a gene that interacts HLA-C in the major histocompatibility complex) in ‡Khomani with some variants originating in this population and subsequently transmitted to neighbouring populations. None of the statistics that we employ suggests selection (or population-specific selection) at this locus however (data not shown). As functional analyses of more immune system genes become available, it may be possible to make more definitive links between genes apparently under selection and their roles in various diseases, both those known to have affected the indigenous populations, such as smallpox and influenza during European colonization, as well as unknown diseases that may have affected the populations during the earlier migrations. It is possible that the genetic variants under selection in the immune genes in ‡Khomani were introduced via admixture (i.e. adaptive introgression) from Bantu-speaking populations, such as the Herero. There is clear evidence of genetic material from Bantu-speaking populations in the ‡Khomani (e.g. [8]). Although we cannot rule out adaptive introgression for single immune regions, we note that (i) the differentiation (FST) between the ‡Khomani and the Herero is very similar for immune genes and other genes (0.0538 versus 0.0528, p > 0.05, t-test) and (ii) contrasting the frequency of the immune SNPs likely to be under selection (table 3) across worldwide populations (electronic supplementary material, figure S8), suggest that the frequency changes have occurred specifically in ‡Khomani (six of seven SNPs showed substantial change, in contrast to one for the Ju|’hoansi and none for the Herero). We hypothesize that the majority adaptive immune gene-variants in ‡Khomani are endogenous. However, regardless of the origin (introgressed or not) of these variants, they are under stronger selection in the ‡Khomani than in the Ju|’hoansi. It is important to note that many genes not associated with immune response are likely to have experienced different selective pressures in the ‡Khomani and the Ju|’hoansi. This would decrease our power to detect a difference between immune and all genes even if there are population-specific differences in selection pressure on the immune system. Additionally, the simultaneous use of three stringent criteria for candidate regions of selection in our second approach potentially omits true signals of selection and differentiation (but is important for avoiding false-positive results). That we find signals of differential selection at immune genes in the ‡Khomani compared to the Ju|’hoansi using two conservative approaches increases our confidence in the results. Our results indicate that selective pressure on immune genes has been strong for indigenous southern African populations, but also that it was a considerably stronger force in the ‡Khomani than the more isolated Ju|’hoansi. The regions with the strongest signals of selection in the Ju|’hoansi contained no immune genes while there were at least eight immune genes in the regions with the strongest signals of selection in the ‡Khomani, supporting the theory of less selective pressure on Ju|’hoansi immune system genes. Our findings suggest that rapid adaptation of immune function can result from contact with external groups and their unfamiliar diseases. We used 2 286 795 high-quality filtered SNPs typed by the Illumina Omni 2.5M SNP array ([8], data available at [36]). Related and exceptionally admixed individuals were removed from the analysis, resulting in sample sizes of 17 Ju|’hoansi individuals and 17 ‡Khomani individuals (for further details of sampling and processing of the data, see [8]). Sampling locations for the Ju|’hoansi and ‡Khomani are indicated in figure 1. We first compared summary statistics for SNPs in immune system genes and SNPs in all genes to examine differences in selection between the two in each population. We created a list of SNPs in all genes using the hg19 gene list (http://genome.ucsc.edu/cgi-bin/hgTables?command=start). We assembled the start and end positions of each gene and combined any overlapping intervals. SNPs in the resulting intervals were selected from the full SNP dataset, yielding the list of all genic SNPs. Defining a set of immune genes is not straightforward as the immune system is complicated and involved in many interactions. We examined several lists of immune genes and chose the Immunome Database [37,38], because its clear inclusion/exclusion criteria define a core immune gene set. Genes included in the Immunome must have a specifically immune function, or if a part of another system, pathway, or interaction, the gene must have a clear role in immune processes [37]. This allowed us to detect selection on immune function while avoiding potential confounding effects of selection on non-immune roles of genes with broader functions. Of the 893 genes in the Immunome, 38 were either on sex chromosomes or only on certain haplotypes of autosomes and so were excluded from the analysis, resulting in a list of 855 immune genes. SNPs in these genes were selected from the full SNP dataset, yielding the list of immune SNPs. We calculated four summary statistics for these two sets of SNPs. We computed the integrated haplotype score (iHS [14], calculated following Pickrell et al. [19]) and the population-specific branch length (PBS [16]) for the Ju|’hoansi and ‡Khomani separately. The relative iHS (XP-EHH [15]) and the BLD were calculated between the Ju|’hoansi and ‡Khomani. Both PBS and BLD rely on a rescaling of pairwise FST values according to the relationship T = − ln(1 − FST). We used Weir & Cockerham [39] to calculate the pairwise FST values between Ju|’hoansi, ‡Khomani, and Herero (a Bantu-speaking population used as the outgroup, n = 8). No conditioning on SNPs being polymorphic in any of the populations was performed. The BLD between Ju|’hoansi and ‡Khomani was calculated as T(Ju|’hoansi,Herero) – T(‡Khomani,Herero). Note that BLD between population ‘pop1’ and population ‘pop2’ using ‘pop0’ as the outgroup, by construction, equals PBS(pop1,pop2,pop0)-PBS(pop2,pop1,pop0). For |iHS|, PBS, and BLD, that are not normally distributed, we used the Mann–Whitney U-test to assess whether the distribution of the statistics were different for SNPs in immune genes and SNPs in all genes. For XP-EHH (which closely follows a normal distribution), we assessed statistical significance using the Student's t-test. To determine what genomic regions were driving the difference between immune SNPs and genic SNPs, we followed Busing et al. [40] in performing a weighted block jackknife analysis. We divided the genome into 5 Mbp blocks, removed each individually, and re-calculated the p-value of the difference between immune and genic SNPs. We ordered the blocks according to how much their removal increased the p-value (lowered the significance). The dataset was then decimated by removing first the top block on this list, then the second block, and so on. For each additional block that was removed, the p-value for the difference between immune and genic SNPs was calculated. This was repeated until the difference between the two categories of SNPs was no longer significant (p > 0.05). In this way, we identified the 5 Mb blocks driving the difference between immune and genic SNPs. To select genes with strong signatures of selection, we combined three summary statistics, examining only genes within 100 kb of SNPs that had high values for all three. Selected SNPs belonged to the top 10 iHS windows, top 1% FST values, and had a PBS branch length ratio above 2.5 (the branch length for a given SNP was 2.5 times as long in the population of interest as in the second-longest population). The top iHS windows were assessed by calculating p-values for non-overlapping 200 kb windows [19]. Adjacent windows with p-values below 0.01 were merged and assigned the lowest p-value among the merged windows. The cut-offs for each test were applied separately to the full dataset. More stringent cut-offs were also investigated for each summary statistic to examine the impact of stringency on the results (table 2). A high value for all three summary statistics indicated that selection had acted on the corresponding genomic region in only one of the populations due to the differentiation required to generate significant FSTand PBS values. SNPs with high values for all three tests were checked for proximity (±100 kb) to immune genes. Gene functions were investigated using GeneCards [41,42], the UCSC Genome Browser [43,44], and literature searches. Genes were considered immune-related if there was strong evidence for a functional role in immune processes. As we were interested in selection due to infectious disease, genes involved strictly in autoimmune or tumour-related disease were not included. This study investigates published data, see Schlebusch et al. [8] for ethical description. This study investigates published data [8] for a data description and access. The data can be accessed at: http://jakobssonlab.iob.uu.se/data/. C.M.S. and M.J. conceived the project; K.A.O., P.Sj., and P.Sk. analysed data; K.A.O., P.Sj., C.M.S., and M.J. wrote the paper; C.M.S. and H.S. conducted community engagement and fieldwork; H.S. contributed the samples. All authors read and approved the manuscript. The authors declare no competing interests. This project was supported by grants from the Erasmus Mundus Programme (to K.O.), the Wenner-Gren Foundations (to C.M.S.), the South African Medical Research Council (to H.S.), Göran Gustafsson foundation, Knut and Alice Wallenberg Foundation and the Swedish Research Council (to M.J.). We are grateful to all subjects who participated in our research and would like to thank members of the Working Group of Indigenous Minorities in Southern Africa (WIMSA) and the South African San Council for facilitating research among the San people. Computations were performed at the Swedish National Infrastructure for Computing (SNIC-UPPMAX). FootnotesElectronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3721801. References
Page 10Habitat destruction is one of the major drivers of biodiversity loss worldwide [1,2]. While habitat loss has obvious immediate and high-impact ecological consequences, habitat degradation, in contrast, has slower, more subtle effects that are more difficult to detect [3]. Coral reefs are ecosystems that are at particular risk from habitat degradation. In these ecosystems, the health of corals are of prime importance as they represent ecosystem engineers, providing habitat to hundreds of animal and plant species [4]. Recent climatic changes, operating through an increased frequency of severe storms and ocean warming, have threatened the health and resilience of these ecosystems [5,6]. In fact, the Great Barrier Reef, the world's largest coral reef system, has recently experienced a period of ocean warming that may leave a tract of 1000 km of coral reefs experiencing 50–90% coral death [7]. While traditionally, biodiversity loss has been assessed through species extinctions, a few have argued that a missed component that often precedes those species extinctions are the alterations of ecological interactions in which these species are engaged [8,9]. Hence, studying changes in the way species interact in degraded coral reef ecosystems could provide insights into the resilience of the community in the face of environmental change. Predation is a major force shaping communities, and has been ascribed a fundamental role in the promotion and maintenance of biodiversity. Due to the highly variable nature of predation, both in space and time, prey have evolved numerous ways to decrease their risk of capture. These adaptations include behavioural, morphological or life-history changes. Predation pressure, for instance, dictates where individuals forage, set up a territory and with whom they mate [10]. Some prey have predator-induced morphologies, such as protective spines or helmets that help reduce their rate of predation [11]. Others show these defensive morphologies from birth [12]. Prey with complex life histories can sometimes alter the timing of their transition from one stage to the next based on predation risk in either stage. For instance, predators capitalizing on eggs can induce prey to hatch earlier than those that are not exposed to predators [13]. Conversely, prey detecting predation risk in the next life stage can delay their transition in order to reach larger sizes before entering the next stage, thus increasing their chance of surviving [14]. Many more examples of phenotypic plasticity exist in response to predation [15]. Such alterations in prey defences have cascading effects, in the form of trait-mediated indirect interactions (TMII). It has been suggested that TMII are more regulatory in prey populations than traditional consumptive, density-mediated interactions [16]. Most TMII are inducible and expressed in a threat-sensitive manner, that is, they are expressed with a ‘strength’ that matches that of the risk perceived. Thus, in order to know when and how much to invest in antipredator defences, prey need to assess their risk of predation using cues from their environment. In aquatic ecosystems, most prey rely on visual and chemical cues to assess risk [17]. Because visual cues are often limited by light availability and by highly complex habitats like coral reefs or kelp beds, and can be manipulated by predators via crypsis, it is not surprising that many aquatic prey have a strong reliance on chemical information to inform them about risk [18]. One of the most common classes of chemicals that aquatic prey use are injured conspecific cues, often referred to as alarm cues in fishes. These cues are present in the skin or tissues of conspecifics and are only released into the water column via mechanical damage to the skin or tissue, as would typically occur during a predator attack. As such, they represent a highly reliable indicator of risk, and are known to elicit immediate and dramatic antipredator responses in nearby conspecifics. These responses are highly conserved and documented in a wide variety of taxa, including corals, molluscs, crustaceans, fishes and larval amphibians [19]. The widespread use of these cues in aquatic systems illustrates the critical role they play for the survival and maintenance of populations. Indeed, these cues have been shown to elicit most trait-mediated indirect interactions discussed above, and many more, such as facilitating learned predator recognition [20]. Not surprisingly, the presence of these cues has been linked to increased prey survival during staged predator–prey encounters [21–23]. As such, these cues are considered a major source of information for prey. Our present study aimed to assess the effect of coral degradation on the ability of coral reef fishes to detect and respond to injured conspecific cues. Previous work suggests that the Ambon damsel, Pomacentrus amboinensis, fails to respond to injured conspecific cues when the cues pass over a patch of degraded coral [24,25]. Recent research also suggests that this species is also unable to learn the identity of novel predators using chemical alarm cues, but a congeneric specialist of dead coral habitats was still able to use information contained within the alarm cues to identify threats [26]. This important ecological difference between closely related species begs the question of how widespread the negative effect of coral degradation on the use of chemical information is to coral reef fishes. Specifically, our first experiment investigated how widespread this phenomenon was, by testing six common and closely related damselfish species, sampling the species from a variety of habitats. Two species, Pomacentrus moluccensis and Chromis viridis, are found on healthy live corals (live coral associates). Two species, P. chrysurus and P. nagasakiensis, are commonly found on coral rubble (dead coral associates), while our last two species, P. amboinensis and P. wardi, are found on mixed habitat types (mixed associates). Each species was tested in both a live and dead coral environment for their response to their species’ injured conspecific cues or a heterospecific control. Predictions from our previous studies were that the alarm cue response of fish that are coral obligates may be most affected by coral degradation, while those more typically associated with dead and degraded habitats may have evolved a mechanism to circumvent the problem. A second experiment was performed to try and tease out the mechanism behind the results of experiment 1, specifically to test whether the lack of response of P. amboinensis in degraded coral water was due to the inactivation of the cues in that environment, or whether it was due to the inability of P. amboinensis to detect the cues via sensory interference. Newly settlement-stage juvenile damselfish (five Pomacentrus species and one Chromis sp.—see electronic supplementary material for more details) were collected overnight using light traps moored in open water around Lizard Island (14′40° S, 145′28° E), in the northern Great Barrier Reef, Australia, in November 2015. The juveniles, sorted by species, were held in 20-l flow-through holding tanks and fed three-times a day with brine shrimp (Artemia nauplii). Apogonids (cardinalfish) were caught on the fringing reef using hand nets and fed fish pellets daily. They were used as heterospecific control (see below). The first experiment consisted of exposing six common species of damselfish juveniles to their injured conspecific cues or a heterospecific control (controlling for the smell of any fish; apogonid) in seawater flowing past either live or dead-degraded coral. The experiment followed a 6 × 2 × 2 completely randomized design. The second experiment investigated possible mechanisms responsible for the loss of response of fish to alarm cues in degraded environments. We chose P. amboinensis and P. nagasakiensis juveniles for this experiment, as the former is affected while the latter appears unaffected by water that has been in contact with dead-degraded coral. The two species were maintained in two habitats (live or dead coral water), and were exposed to each other's injured cues or apogonid cues in a 2 × 2 × 3 completely randomized design. We predicted that if the absence of response of P. amboinensis is mediated via a deactivation of its alarm cues (hypothesis 1), then neither species should respond to P. amboinensis cues, while they should both respond to P. nagasakiensis cues. If, on the other hand, P. amboinensis cannot respond to its alarm cue due to sensory interference (hypothesis 2), then we predicted that P. amboinensis should not respond to the alarm cues from a closely related species, while P. nagasakiensis should respond to both cues. The protocol used to test the fish was identical for both experiments. Groups of 10 juveniles were placed into 12-l plastic exposure tanks, which had flowing seawater from a header tank containing either live or dead coral. The header tank consisted of a 15-l Amundsen bucket containing either a piece (approx. 60 cm in circumference) of healthy, live Pocillopora damicornis, a hard bushy coral commonly found at our field site, or an equal sized piece of dead-degraded coral that was encrusted with algae. The header tanks were equipped with an airstone, and had constantly flowing fresh seawater at a rate of 1 l min−1 (one tank turnover every 12 min). The header tank was plumbed in such a way that allowed the overflow to enter the exposure tanks. Both coral types were changed daily. The fish were kept in the exposure tank for 48 h before the test phase. Following the exposure phase, fish were moved individually into 5-l plastic tanks, equipped with a sand substrate, a moulded plastic replica of branched coral (15 cm high) serving as shelter, and an air stone, to which was attached a 1.5 m long injection hose. A 4 × 4 cm grid was drawn on the tank to facilitate data collection. Each test tank received flow-through water from a header tank containing live or dead coral, as described above. The difference was that the flow-through from the header tank was divided among five testing tanks. Each test tank thus received water at a rate of approximately 1 l/5 min (one tank turnover every 25 min). The fish were left to acclimate overnight and were tested the following day. The bioassay followed established protocols [19] and is described in details in the electronic supplementary material. In short, the behaviours of each fish (number of feeding strikes and line crossed, as measures of feeding and activity) were observed for 3 min before and after the introduction of a stimulus (5 ml of alarm cues or apogonid cues). Reductions in feeding and activity are both well-established antipredator responses. We tested 244 fish (n = 10–11/treatment) in experiment 1 and 148 fish (n = 12–13/treatment) in experiment 2 (see electronic supplementary material for size). The observer was blind to the treatment and the order of treatments was randomized. Given that feeding and activity are not independent variables, the two were analysed simultaneously using a MANOVA approach. Pre-stimulus data were first analysed to ensure there was no difference among treatment groups prior to stimulus injection. Pre- and post-stimulus data were then used to calculate a per cent change in behaviour [(post-pre)/pre] and the resulting variables were used in subsequent analyses. For experiment 1, both analyses (one for prestimulus baseline, one for behavioural change) were carried out using a three-way MANOVA, testing the effects of species, habitat type (dead versus live coral) and cue type (heterospecific versus conspecific cues) on behavioural responses. Subsequent two-way MANOVAs were used to explore possible interactions. For experiment 2, both analyses were performed using a three-way MANOVA, testing the effect of species (P. amboinensis versus P. nagasakiensis), coral type (live versus dead) and cue type (P. amboinensis, P. nagasakiensis or apogonid control). Subsequent two-way MANOVAs and Tukey HSD post-hoc comparisons were performed to explore interactions. For all tests, data met parametric assumptions. The only factor explaining differences in pre-stimulus values was species (Pillai's Trace: F10,440 = 3.1, p = 0.001), indicating that fish from the same species exposed to different coral waters or given different cues did not differ in their baseline activity levels. No other factor was significant (all p > 0.4). Change in behaviour was influenced by a three-way interaction among species, cue and coral (Pillai's Trace: F10,440 = 3.6, p < 0.001, figure 1). Splitting the analysis between the two coral types revealed that, in live coral, all fishes responded to conspecific cues with a significant antipredator response (cue: Pillai's Trace: F2,105 = 152.4, p < 0.001). We failed to find an effect of species (Pillai's Trace: F10,218 = 1.3, p = 0.3) or an interaction between cue and species (Pillai's Trace: F10,218 = 1.3, p = 0.2), indicating that all species responded similarly to their respective alarm cues. On dead coral, however, a significant species by cue interaction (Pillai's Trace: F10,222 = 3.3, p = 0.001) indicated that species differed in their responses to alarm cues. Species found on live coral failed to respond to their alarm cues in dead coral (P. moluccensis: F2,18 = 1.3, p = 0.3; Chromis: F2,18 = 0.7, p = 0.5). Dead-degraded associates, on the other hand, maintained their response to alarm cues (P. chrysurus: F2,17 = 54, p < 0.001; P. nagasakiensis: F2,18 = 25, p < 0.001). Interestingly, species living in mixed habitats showed mixed responses, with P. amboinensis failing to respond to alarm cues (F2,17 = 0.2, p = 0.8), and P. wardi displaying a full antipredator response to the alarm cues (F2,18 = 5.7, p = 0.012). Figure 1. Mean (±s.e.) proportion change in the number of feeding strikes (top panel) and line crosses (bottom panel) for damselfish species that are typically associated with live coral only (live coral associates), degraded-dead coral only (dead coral associates) or species that settle in both types of habitat (mixed associates). Fish were maintained in water from either live or dead coral and exposed to cues from heterospecific apogonid (empty bars) or cues from injured conspecifics (solid bars) (n = 10–11/treatment).
None of the treatment groups differed in their pre-stimulus baseline (Pillai's Trace: all p > 0.4). Change in behaviour was affected by an interaction between coral and test cue (Pillai's Trace: F4,272 = 12.1, p < 0.001, figure 2). Splitting the analysis by coral revealed that, in live coral, the responses of the fish were affected by an interaction between species and test cue (F4,136 = 2.8, p = 0.028). Specifically, both species displayed a significant antipredator response to the Pomacentrid alarm cues compared to the apogonid control (Tukey post-hoc comparisons: P. amboinensis versus apogonid: p < 0.001, P. nagasakiensis versus apogonid: p < 0.001 for both variables). However, each species responded to their own cues with a stronger intensity than to the one of the close relative (2 × 2 MANOVA: species × cue interaction: F2,45 = 4.3, p = 0.02). In dead coral, however, the pattern was different. Fish behaviour was affected by the type of cue they received (F4,136 = 20.1, p < 0.001), but there was no species by cue interaction (F4,136 = 0.7, p = 0.6). Both species responded with a significant antipredator response to P. nagasakiensis cues compared to apogonid cues (Tukey post-hoc comparisons: p < 0.001 for both variables), but failed to show a statistically significant response to P. amboinensis cues (p = 0.08 and p = 0.8 for feeding and activity respectively). For P. nagasakiensis cues, once again, the response from conspecifics was stronger than that of close relatives (p = 0.032). Figure 2. Mean (±s.e.) proportion change in the number of feeding strikes (a) and line crosses (b) for P. amboinensis or P. nagasakiensis tested in live coral water (top panels) or dead-degraded coral water (bottom panels). The fish were exposed to cues from distantly related apogonid (control, empty bars), cues from P. amboinensis (light grey bars) or cues from P. nagasakiensis (dark grey bars).
Coral degradation had dramatically different effects on the efficacy of alarm cues among closely related species. Of the six species tested, half maintained their response to injured conspecific cues in degraded corals, while the other half completely lost their ability to respond to the cues in the degraded habitat. This was a striking result, because although the composition of the active substance in the alarm cues is still unknown (and likely different for all species, since we do not see taxa-wide responses to a single cue), a number of studies found that these compounds were highly conserved among closely related species. For instance, several species of salmonid trout from a few genera can respond to each other's cues, although the strength of the response decreases with increased phylogenetic distance [27]. Similar results are found in other species, including damselfish [28]. Our findings emphasize that the interaction between the background olfactory landscape and chemical alarm cues is species specific and can differ between closely related fish. Results suggest for P. amboinensis that the lack of response in a degraded environment may stem from a deactivation of the active component of their alarm cue. Indeed, while the expected cross-species response is intact in live coral environments, neither P. amboinensis nor P. nagasakiensis can respond to P. amboinensis cues in degraded coral. Interestingly, they can both respond to P. nagasakiensis cues in that same environment. That result suggests that P. amboinensis alarm cues are modified by the chemistry of water from degraded corals, while the same water does not affect P. nagasakiensis cues. We speculate that a chemical group nearby the active site of P. amboinensis’ cue either changes its conformation or binds with a water-borne compound, which blocks access to the active site, rendering the cue inactive. Another potential explanation for our results would be that the responses to species-specific alarm cues are mediated by species-specific receptors in the olfactory rosette, and that degraded coral water contains a compound that would block the receptors for P. amboinensis alarm cues in the rosette of both species. While technically possible, the information we have to date with regards to olfactory perception and neurobiology [29], the multi-compound nature of the alarm cues [30] and the principle of parsimony makes this alternative explanation less likely in our opinion. Exploring both these suppositions would require some advances in the field of vertebrate predation-related chemical ecology. The chemistry of these interactions remains sadly understudied [19,31]. Based on the previous findings, one of two scenarios, not necessarily mutually exclusive, may explain the pattern of responses we observed. First, the pattern of response follows that of the species’ habitat. Although we cannot test this hypothesis rigorously, our limited sample size (n = 4 species) provides preliminary evidence that habitat may be a good predictor of the impact of coral degradation on cue use. Both species typically associated with live coral lost the ability to respond appropriately to injured cues in a degraded habitat, while both species typically associated with rubble and dead coral maintained the appropriate cue response. This pattern was also found for P. coelestis, a dead coral associate [26]. Hence, the different sensitivity to degraded coral habitat could stem from local adaptation to microhabitat conditions, a hypothesis already present in the literature [32,33]. Rubble has always been a part of coral reef ecosystems. When corals die, their exoskeletons break apart and form rubble-dominated microhabitats, until new corals recruit and take over. Species that live in those habitats may have selected the habitats due to the combined benefits from lower competition and their unique ability to detect alarm cues, an ability that was inherently present or was selected for by predation-mediated natural selection. The second scenario that could explain the pattern of response is phylogeny. Two relatively recent studies have defined the phylogenetic relationship among Pomacentridae [34,35]. Both of them have relationships among four of our species, but neither of those have tested P. wardi. From these two papers, we can make some general groupings: Chromis viridis is the most distantly related, P. moluccensis and P. amboinensis are sister species, and P. nagasakiensis and P. chrysurus are also closely related to each other. This phylogenetic pattern also matches our response patterns, with P. nagasakiensis and P. chrysurus maintaining their response to injured cues in degraded habitat, while P. moluccensis and P. amboinensis both lost their responses in the degraded habitat. Interestingly, according to Cooper et al. [34], P. coelestis is closely associated with P. chrysurus, and we see concordance in the response pattern of the two species in degraded habitats. It is difficult to conclude anything for the other species. Following the principle of parsimony, the change seen from a phylogenetic point of view may in turn explain the ecological segregation of the species based on their ability to use predation-related cues in degraded habitats. For the species that lost their response to alarm cues, the ecological consequences are likely significant, with a potential decrease in all alarm cue-mediated indirect effects. The immediate effect of alarm cues is to warn nearby conspecifics of a recent predation attack. The increase in vigilance results in an immediate increase in survival over the next several minutes to hours [21]. However, alarm cues also facilitate learning and other lasting effects including investment in morphological defences [36,37]. Without these cues, these species will likely be much more vulnerable to predation. Many coral reef species, including our damselfishes, have a bipartite life history where pelagic larvae recruit to the reefs after 10–25 days and settle to become benthic juveniles. This transition is linked to a predation-mediated population bottleneck whereby 60–90% of juveniles are consumed within the first 2 days post-settlement [38]. There is immense selection for prey that can use alarm cues to reduce risk of predation. The loss of these cues by some members of the community will have far-reaching consequences for restructuring the community. For instance, the cross-species responses seen in our second experiment may indicate benefit for some species to associate with other species that can provide them with valuable public information regarding predation risk, such as would happen during cross-species social learning among guild members [39,40]. The present study provides a viable mechanism that explains the relatively rapid loss of species from systems where hard corals have died, despite the maintenance of topographic complexity for years after death. It provides a link between the expansion of dead-coral-dominated landscapes and their rapidly altered communities, such as those seen in the Caribbean [41]. A common pattern seen in many ecosystems is that generalist species that are able to survive in modified habitats have a competitive edge over specialists in the face of habitat change [42–44] and these species make up the new, modified community in altered environments. Our results provide evidence that some coral reef fish species are functionally more generalist than others, as demonstrated by their ability to use predation-mediated cues in both pristine and degraded coral environments. As such, we predict that these species will make up a higher proportion of the fish community in the reefs of the future, and that those that cannot adapt may slowly disappear. All work followed James Cook University Animal Ethics Protocols A2080 and A2005. The data have been uploaded as electronic supplementary material. M.C.O.F., D.P.C. and M.I.M. designed the project; M.I.M. and B.J.M.A. collected the fish and the corals; M.C.O.F. and D.P.C. collected the data; M.C.O.F. analysed the data and wrote the first draft; all authors contributed to the final version. We have no competing interests. Funding for this work was provided by the ARC Center of Excellence of Coral Reefs (M.I.M.) and the Natural Sciences and Engineering Research Council of Canada (M.C.O.F., D.P.C.). We thank the staff of Lizard Island Research Station. FootnotesElectronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3726595. © 2017 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited. References
Page 11An escape response allows prey to evade predators with rapid locomotion [1]. Because of its potential to directly affect survivorship, natural selection may favour animals that can execute an escape response with high locomotor performance. Indeed, the physiology and mechanics of locomotion features many traits that appear to be adaptations for a rapid response and fast motion. Escape responses are controlled by large-diameter command neurons (e.g. the giant axon of squid [2]), which often recruit specialized muscles (e.g. the axial musculature of fish [3]), which may animate an appendage that is dedicated to escape behaviour (e.g. the uropods of crayfish [4]). Prey may direct this escape in an optimal direction [5], or may alternatively benefit from moving randomly so as not to be predictable [6,7]. However, it does not necessarily follow that any enhancement in speed or variation in heading will have a positive effect on a prey's survival. Fish predators commonly approach their prey at a relatively slow speed [8,9] and this could permit escape by a prey operating below its maximal capacity. The aim of this study was to test whether improvements in kinematics related to locomotor performance and sensing affect prey survival by examining predator–prey interactions in zebrafish (Danio rerio, Hamilton 1922). We addressed this aim with a novel approach that combines experimentation with mathematical modelling. Our methodology was developed to meet the challenges to understanding the coupled dynamics of predators and prey. Owing to the sensing of both animals, the actions of the prey may (or may not) be a response to the predator, which may (or may not) be a response to prior motion by the prey. Regression analysis is generally insensitive to such interdependency, yet may succeed in resolving dominant features of successful prey [10] or predators [11]. It has additionally been helpful to study behavioural responses to an artificial predator or prey that is experimentally controlled and therefore not coupled to the actions of the animal [11–15]. An alternative approach has attempted to formulate a behavioural algorithm of one animal by considering their responses to the measured kinematics of the other [16]. Our present approach similarly included measurements of predator–prey kinematics, but these were used as a basis for an agent-based probabilistic model that calculated the trajectories of both animals from a series of behavioural actions. A series of simulations by this model allowed for a predictive consideration of the effects of differences in behaviour on prey survival. Zebrafish provide an excellent experimental subject for studying the mechanisms of predator–prey interactions. The larval stage of this species serves as a model for the neurophysiological [17–19] and biomechanical [20,21] basis of behaviour. Predator–prey interactions may be experimentally replicated in the lab, where adults and juvenile zebrafish strike at larval zebrafish with suction feeding and the larvae respond with a fast-start escape response [22]. Therefore, the interactions between zebrafish of different life-history stages replicate the principle predator and prey behaviours that characterize a broad diversity of piscivorous interactions [5,10]. When approaching an evasive prey, adult zebrafish move much slower than their maximum speed [22], which is common among suction-feeding fishes [8,9]. A slow approach presumably allows greater control over the direction and timing of the suction feeding, which is limited to a brief duration over a small region in front of the mouth [23,24]. The prey, by contrast, respond with an explosive escape response with speed that generally exceeds that of the predator. As suggested by prior experiments [25] and theory [5], the relative speed and size of predator and prey has the potential to affect the strategy of both animals. Therefore, we performed experiments with juvenile (three- to four-months old, 2.0 ± 0.4 cm, N = 19) and adult (more than or equal to nine-months old, mean ± 1 s.d. = 3.4 ± 0.5 cm, N = 19) predators with a nearly twofold difference in body length to examine the effects of scale. All experiments were conducted using wild-type (AB line) zebrafish larvae (5–7 days post fertilization) as prey (D. rerio). The fish were placed in a hemispherical aquarium (Ø = 8.5 cm) composed of white acrylic, which served as a translucent diffuser of IR illumination (940 nm), provided by three lamps (CM-IR200-940, CMVision, Houston, TX, USA), positioned below (figure 1a). These lamps provided high-intensity illumination that was invisible to the fish [26], while visible illumination at low intensity was provided by overhead fluorescent lights. High-speed video cameras (1000 fps at 1024 × 1024 pixels, FASTCAM Mini UX50, Precision Photron Inc., San Diego, CA, USA) were each fitted with a 55 mm lens (f/2.8 Micro Nikkon AIS, Nikon Inc., Melville, NY, USA) and positioned at a distance that permitted a view of the entire aquarium, where they recorded both predator and prey. Predation experiments were performed by recording the swimming of one predator and one prey fish in the aquarium (figure 1a). This began by placing the fish on opposite sides of a partition. After acclimation (15 min), we lifted the partition and recorded the fish until the predator successfully ingested the prey. Figure 1. Kinematic measurements and mathematical modelling for studying predator–prey interactions in zebrafish. (a) Three high-speed video cameras recorded one larval prey and one predator fish (adult or juvenile) that were placed in a hemispherical aquarium. A representative video frame shows an adult in close proximity to the prey. In the inset, orange markers denote the locations of morphological landmarks for the eyes of the predator (+) and the swim bladder in the prey (open circle). (b) A state diagram illustrates the major components of the agent-based probabilistic model used to simulate the interactions between predators and prey (see table 1 for symbol definitions and parameter values). Each fish behaves according to an algorithm specific to a particular behavioural state and the probability of transitioning between states is determined by random-number generators with probability distributions matching our measurements (figure 2). Predators (in red) operate between tracking (T) and striking (S) states and prey are either resting (R) or escaping (E). The outcome of a strike is determined by the capture probability (C, equation (2.2)). Simulations of this model were performed with a Monte Carlo method to generate probability distributions of prey survival (see Material and methods for details).
Figure 2. Descriptive statistics of swimming kinematics. (a–e) Measurements of the probability density (circles) for the kinematics of the prey (a–c) and predator (d–e) for experiments that included either a juvenile (purple) or adult (green) predator. Points denote binned measurements for probability density with a sample size determined using the Freedman–Diaconis rule, which yielded N ∼ 17 measurements per point (see table 1 for total sample sizes). Measurements of escape angle (a) were not significantly different between the two types of experiments and these measurements were therefore combined. For each parameter, we performed a nonlinear least-squares fit to the measurements for a log-normal probability density function (equation (2.1)). The log-mean and log-standard deviation values from these fits (table 1) were consequently used to describe the probability of events in our mathematical model (figure 1b). (f) The capture probability was measured as a function of distance between the predator and prey and we similarly performed a curve fit to approximate this relationship (equation (2.2)) for our model.
Our video recordings were used to perform measurements of three-dimensional kinematics. We calibrated the cameras using direct-linear transform (DLT) using ‘Digitizing Tools’ software in MATLAB (2015a, MathWorks, Natick, MA, USA) [27]. We used the position of the predator's two eyes to calculate a mean position that approximated the buccal cavity (figure 1a). The posterior margin of the swim bladder was found on the prey's body, which approximates the centre of mass [28] and the heading was measured by matching an ellipsoid to the body. We acquired these landmarks at five key events in each interaction between predator and prey: (i) the initiation of a predator's approach towards the prey, the (ii) opening and (iii) closing of the predator's mouth during a strike, the (iv) initiation and (v) completion of the prey's escape response. Descriptive statistics were used to characterize the probability of actions by the predator and prey. We recorded the predator-specific parameters of the strike distance (s), the distance from the prey at which a strike was initiated, and the strike duration (τ), the period between the opening and closing of the mouth during suction feeding. For the prey, we found the reaction distance (l), the distance from the predator at which the escape response was initiated and the escape angle (θ), the change in heading from the resting orientation to the escape path. The escape duration (η) included the C-start and subsequent undulatory swimming, until the larva ceased moving. The frequency distribution for each of these parameters was found to be well approximated by the following log-normal probability density function: 2.1 where x is a particular behavioural parameter (s, τ, l, θ or η), μ is the log mean and σ is the log standard deviation. We determined best-fit values for μ and σ for each behavioural parameter by maximum likelihood to binned measured values. The bin size was determined by the Freedman–Diaconis rule, which yielded a number of samples per bin of around 17 measurements. We used the the non-parametric two-sample Kolmogorov–Smirnov test (i.e. KS-test) [29] for statistical comparisons between parameters because they failed to conform to normal distributions.The probability that the strike of a zebrafish predator is successful depends critically on the distance between the mouth of the predator and the prey [22]. We therefore binned measurements of capture probability (C) and approximated its dependency on distance with the following sigmoidal function: 2.2 where d is the strike distance, d0 is the decay distance and r is the decay rate. This function additionally approximates the spatial variation in fluid forces that act on prey when subject to suction feeding [30]. The best-fit values for d0 and r were determined by least squares.An agent-based probabilistic model was developed to simulate the conditions of our experiments. This model predicted the two-dimensional motion of a predator and prey [31] according to algorithms that were specific to the behavioural state of each animal (figure 1b). The predator's states were tracking (T) and striking (S) and the prey's were resting (R) and escaping (E), consistent with prior observations [13,22]. The duration of states and probability of transitioning between states were determined by random-number generation that conformed to the log-normal probability distributions (equation (2.1)) and range of measured values. Prey capture was similarly found with a probability that depended on the proximity to the prey distance at mid-gape (equation (2.2)). A simulation was terminated if a strike was successful, otherwise the predator reverted to the tracking state (figure 1b). Single values for speeds and latencies were used for all simulations (table 1), determined by trial and error to replicate measured survivorship and approximate prior measurements [22,32]. The prey escape angle was directed with respect to the right or left side of the body with a probability of moving contralateral to the predator (table 1) as previously measured [13]. Predator and prey kinematics were calculated at a fixed time step (5 ms) and simulations began with random positioning of the prey within one aquarium diameter (Ø = 8.5 cm) of the predator. Simulations were scripted in MATLAB to calculate the motion of both agents and their behavioural states, which consequently determined the number of evasions before prey capture. This model represents a Markov chain that treated the predator and prey's actions as probabilistic, but each outcome also depended on the determinism of their kinematics. Each behavioural state operated by a distinct set of rules. In the tracking state, the predator moved at a fixed approach speed with a direction always headed towards the prey with a time delay (figure 1b). The prey transitioned from rest into the escaping state when the predator moved within the reaction distance with a latency [33]. The prey escaped along a straight path with speed that varied as a single saw-toothed pulse, with the maximum value (the peak of the sawtooth) attained at 20% of the escape duration. This well-characterized prey speed, as determined by a frame-by-frame analysis of escape swimming for 12 larvae. This model simplified many aspects of the complexity of predator–prey interactions. It assumed two-dimensional kinematics that were not bounded by a barrier. We tested the model's predictions by comparing its predictions of survivorship for 1000 simulations against our measurements (two-sample Kolmogorov–Smirnov). Survivorship was defined as the number of individuals surviving a particular number of strikes, divided by the initial population. We tested whether any single probability distribution was sufficient to predict observed patterns of survivorship by drawing parameter values from a single parameter distribution while leaving all other parameters fixed. These simulations failed to match the measured patterns of survivorship that was achieved by drawing from distributions for all five possible parameters (electronic supplementary material, figure S1), which suggests that variation in all parameters is necessary for a predictive model. We performed an analysis of the model to evaluate the parameters that had the greatest effect on prey survival. This was achieved by running a Monte Carlo series where individual parameters were varied in increments of 10% between −90% and 100% of their original mean value. For parameters described by a probability distribution, the log-mean parameter, μ, was adjusted to create the desired per cent change in the mean of the distribution with a fixed σ. The range of possible random values for each distribution was also adjusted to retain the same cumulative probability range. The effect of these manipulations was assessed by comparing survivorship against the model's prediction without any parameter variation using a Kruskal–Wallis test. The behaviour of both predator and prey was similar whether the predators were juvenile or adult zebrafish. Prey responded in both cases with indistinguishable differences in escape angle (KS-test: p = 0.86, N = 164). Prey reacted at a mean distance to juvenile predators ( Although they exhibited similar behaviour, adult predators were more effective than juveniles. Juveniles did not succeed in capturing prey beyond a distance of 3.2 mm (N = 91), whereas adults were successful at a maximum distance that was about three times greater (10.4 mm, N = 77) and showed a decay distance that was greater by about the same factor (table 1, figure 2f). We found that this result was unaltered by excluding the furthest point measured for capture distance among adult predators (figure 2f, KS-test: p = 0.28). We tested whether the superior feeding of adults was due to juveniles approaching the prey with inferior accuracy by measuring the bearing, the prey's radial position relative to the predator's heading. The bearing when prey initiated an escape was not significantly different (KS-test: p = 0.15) between juveniles (N = 91) or adults (N = 77). At mid-gape, however, prey succeeded in evading juvenile predators to the extent that the median bearing (13.1°), was less than half that of adult predators (30.0°), which was a significant difference (KS-test: p < 0.01). Therefore, adult predators were more capable of adjusting their strike towards escaping prey than juvenile predators. This ability was included in our model through the measured probability distributions of capture success (figure 2f). The model predicted kinematics and prey survivorship that were similar to our measurements. The temporal sequence of events in the model offered a reasonable approximation of the kinematics of live predator–prey interactions (figure 3a,b). Measured survivorship decreased monotonically from the first strike and extended to as many as 20 strikes, with a slower decline in survivorship for juvenile predators (figure 3c). The survivorship predicted by the model was statistically indistinguishable for both adult (KS-test: p = 0.93, N = 73) and juvenile (KS-test: p = 0.86, N = 91) predators. Furthermore, all trends from the parameter analysis of the pursuit–evasion model were similar between the adults (figure 4) and juveniles (electronic supplementary material, figure S2). Figure 3. Comparison between experimental measurements and modelling. (a) Trajectories of predator and prey from a representative experiment (i) and simulation (ii). The position of predator and prey that correspond to particular time points are shown with connecting arrows. (b) Ethograms for these trajectories illustrate the temporal changes in the predator's swimming and strike (i), which are respectively modelled by the tracking (T) and striking (S) states (ii). The prey's behaviour while motionless and during escape (i) were respectively modelled as resting (R) and escape (E) modes (ii). For both ethograms, the distance (d) between predator and prey are shown. Particular moments in the trajectories are highlighted with vertical lines that correspond with the same-coloured arrows in (a). (c) The probability that a prey survives over a particular number of strikes is shown for adult (i) and juvenile (ii) predators for experiments (dark grey) and simulations (light grey).
Figure 4. Parameter analysis of the mathematical model to examine the effects on escape probability. (a) We individually varied the mean of a parameter among simulations by manipulating its probability distribution (figure 2, see table 1 for parameter definitions and values). Each point represents the survival probability of prey among 1000 simulations and filled circles denote a significant difference (KS-test: p < 0.05) from the observed probability. Simulations that varied in escape angle (θ) differed by a interval of 0.127 rad. (b) Variation in escape probability was examined with respect to both escape speed and reaction distance. All simulations shown here used an adult predator, although similar results were obtained with a juvenile predator (electronic supplementary material, figure S2).
Our parameter analysis revealed that escape speed and reaction distance were the only parameters that showed any noteworthy effect on survival. Changes in escape duration, escape direction and escape angle yielded statistically insignificant or otherwise small changes in escape probability (figure 4a). An increase in escape speed similarly had a negligible effect on survival, but survival probability did decline when speed was reduced by 50% or more. By contrast, an increase in reaction distance caused escape probability to elevate by as much as 16% and a decrease of at least 30% had a dramatically adverse effect on survival (figure 4a). A two-dimensional parameter analysis of reaction distance and escape speed (figure 4b) showed little evidence for an interactive effect between these two parameters. We found that the survival of larval fish does not increase by escaping at a faster speed or by varying direction, but only by responding from a greater distance. These results were derived from an agent-based probabilistic model (figure 1b) that calculated the trajectories of predator and prey and the outcome of predatory strikes (figure 3a,b). This model successfully replicated prey survivorship (figure 3c) by simulating behavioural actions that matched our measurements (figure 2). Our analysis of its predictions suggests that prey survival in larval fishes may only be enhanced by increasing the sensitivity of predator detection. The survival of prey depends largely on the actions of the predator. In contrast with the explosive speed of an escape response [20], adult zebrafish tend to approach their prey substantially slower than their capacity, often by braking [32]. We found that the approach speed amounted to less than one-third the maximum speed of escaping larvae (table 1), which is consistent with previous measurements [22]. This approach relates strategically to the mechanics of feeding. The suction feeding of fishes succeeds in capturing prey in only a small region around the mouth over a duration of merely tens of milliseconds [34–36]. A slow approach is common among suction-feeding fishes and is likely a means of enhancing strike accuracy [8,9]. This style of predation is seen in many species of fish [37]. Furthermore, our data suggest that adult and juvenile zebrafish are more likely to capture when approaching larval zebrafish with a slower approach (electronic supplementary material, figure S4). Therefore, the limited range of suction feeding may constrain some predators to a slow approach while offering prey an opportunity to escape [24]. Despite this strategic advantage for prey, adult zebrafish captured prey on the first strike more than one-quarter of the time and rarely needed more than three strikes to be successful (figure 3c). The effectiveness of an escape has previously been considered by classic pursuit–evasion models of fish predation. This theory resolves how the direction of an escape affects the distance between predator and prey [5,31], generally by modelling a single encounter with the assumption that both animals move with a fixed heading and speed. A recently developed version of this model suggests that animals like zebrafish operate in a ‘slow-predator’ strategic domain [38], where the predator moves more slowly than the prey. In this domain, no optimal escape angle exists and prey may evade predators with a broad range of effective escape directions. A faster escape speed serves only to modestly expand the range of effective escape angles. Consistent with these predictions, our model found a monotonic decrease in survival as we reduced escape speed below half of the observed value (figure 4a). We additionally found only modest differences in survival between experiments using adult and juvenile predators (figures 2 and 3c), despite a nearly twofold difference in body size and speed. These predictions would not hold for cases where the predator is faster than the prey. In such a strategic domain, an optimal escape angle arises for the prey and failure to move in that direction is predicted to adversely affect survival [5,38]. Ram-feeding fishes strike at prey while swimming at a relatively high speed and may thereby place prey at a strategic disadvantage. Success in ram feeding may, in-turn, require superior coordination in directing and timing a strike [11]. Ram feeding in this way shows greater similarity in strategy to flying predators such as birds [15], bats [16] and insects [39]. Prey may benefit by escaping in a direction that conforms to an optimal value [5], by being unpredictable [6] or escaping along a trajectory with a small radius of curvature [7,40]. The reaction distance has broad strategic significance. Classic pursuit–evasion models support the simple notion that prey are more evasive if they start from further away [5,31,38]. This principle is consistent with evolutionary models that contrast the fitness benefit of responding from a distance against its potential costs [41,42]. For example, escape responses that are initiated unnecessarily may be energetically expensive, prohibit foraging or succeed in revealing cryptic prey [43]. Responding from a great distance may even be inferior on purely strategic grounds. A prey that is slower than a predator, but capable of executing a tight turn, may be more evasive when initiating this manoeuvre at the final moments of a predatory strike, rather than providing the opportunity for the predator to adjust course [7]. Therefore, a greater reaction distance offers a clear strategic benefit in zebrafish (figure 4), but may not be universally advantageous. The primacy of reaction distance underscores the strategic importance of predator detection. Responding to a predator from afar depends on the sensitivity of receptor organs and the capacity of the nervous system to rapidly recognize a threatening cue and trigger an escape response. As in invertebrate zooplankton [14] and insect prey of spiders [44,45], zebrafish larvae use flow sensing to detect the bow-wave of flow generated by an approaching predator [22]. Flow-sensing may be augmented by olfactory cues [46], though zebrafish do not acquire a sensitivity to the alarm pheromone Schreckstoff until a later stage of growth (more than 48 dpf) [46]. Using the lateral line system, zebrafish larvae may respond to flow up to a distance of 1.3 cm [13] ahead of a gliding zebrafish adult. This range encompasses many of the responses that we recorded (figure 2c), which supports a role of the lateral line in our experiments. This would agree with previous experiments that have shown flow sensing to be necessary for survival in zebrafish larvae [22]. Flow sensing offers the capacity to trigger an escape with a very brief (less than 10 ms) latency [47]. By contrast, a looming visual stimulus succeeds in stimulating an escape response in zebrafish larvae [48], but the demands for visual processing necessitate a latency that is at least 20-fold longer than for flow sensing [49]. Nonetheless, the greatest reaction distances that we observed were outside the range of the lateral line and therefore likely generated by the visual appearance of the predator. The visual system consequently offers prey fish the means to enhance survival by responding to a predatory threat from a distance (figure 4). We found that zebrafish larvae operate in a slow-predator strategic domain when preyed upon by adults and juveniles of the same species. As a consequence, increasing the speed or varying the direction of an escape response has a negligible effect on survival (figure 4). Survival may instead be enhanced by initiating the escape from greater distance by rapidly identifying the predator as a threatening visual stimulus. These findings offer valuable insight into the key strategic factors that govern predator–prey interactions in a diversity of similar fishes and other animals that operate with similar strategy. Data related to this study may be accessed through the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.mr720 [50]. The study was designed in collaboration between A.N. and M.J.M.H. A.N. and C.N. performed all experiments and kinematic analysis. The agent-based probabilistic model was created by A.N., with guidance from M.J.M.H. The manuscript was written collaboratively by A.N. and M.J.M.H. We declare we have no competing interests. This research was supported by grants to M.J.M.H. from the National Science Foundation (IOS-1354842) and the Office of Naval Research (N00014-15-1-2249). A. Soto helped with the model. The manuscript was benefitted from feedback from A. McKee, A. Carrillo, A. Peterson, M. Byron, B. Free, D. Paley, R. Blob and two anonymous reviewers. FootnotesElectronic supplementary material is available online at http://dx.doi.org/10.6084/m9.figshare.c.3729922. References
Page 12Three-dimensional (3D) digital morphological data are commonly employed by palaeontologists and biologists in research. In palaeontology and anthropology, the widespread application of tomography (especially X-ray computed tomography, CT), laser and structured light scanning, and photogrammetry has revolutionized the study of morphology [1–4]. In biology, optical microscopy, magnetic resonance imaging (MRI) and contrast-enhanced CT are important tools for investigating soft-tissue anatomy [5–10]. The revolution brought about by these technologies has increased the amount and detail of anatomical information recovered from fossil and living organisms, transforming the nature of scientific enquiry in related fields. The resulting datasets are often reconstructed and presented as 3D digital models, which are themselves sometimes used in downstream analyses, including geometric morphometrics [11,12], finite element analysis (FEA) [13], multibody dynamics analysis (MDA) [14] and computational fluid dynamics (CFD) [15], thereby facilitating quantitative tests of functional and evolutionary hypotheses [3]. These types of studies have yielded important advances in our understanding of the anatomy of living and fossil organisms (e.g. [10,16,17]), as well as fundamental aspects of their biology, from feeding mode [18–20] to mobility [21,22], development [23,24] and physiology [25–27], as well as developments in taxonomic practice [28,29]. Barriers to data sharing and access to specimens can be eroded because data exist as digital files that can be easily copied and readily distributed, allowing simultaneous analysis by multiple researchers [30]. These attributes should also enhance the verifiability and reproducibility of studies, facilitating the reuse of data and metadata, more in-depth interrogation of any given dataset, and broader-scale comparative analyses through the assembly of large datasets of multiple specimens or taxa. However, authors of studies involving 3D digital datasets of biological and palaeontological specimens often do not publish their supporting data, meaning that results and conclusions cannot easily be verified or replicated, and that this potentially valuable source of novel data cannot be further explored [30]. Ultimately, digital data collected but unpublished are likely to be lost to science [2,28]. This also represents a substantial waste of financial and other resources, and places vulnerable original specimens at greater risk of damage or loss, as the same specimens are likely to be reimaged repeatedly to enable different groups of workers to reproduce the data [28,31]. Consequently, the promise of 3D digital data has not yet been fully realized. This is not news [2,28,30]. However, most national and international funders have imposed regulations on data access and sharing that are forcing researchers and institutions to finally confront this challenge [32]. These regulations range from funder-mandated full release of all data [32], through declarations that the data are available from authors on request, to no release of supporting data [32]. When data are released, they are deposited in a diversity of online databases (e.g. BIRN, Dataverse, Dryad, EOL, figshare, GigaDB, Github, MorphoBank, MorphoDBase, MorphoMuseuM, MorphoSource, Phenome10 K, Zenodo), institutional and funder repositories, physical museums, and research group websites. At least in part, this diversity of approaches reflects uncertainty about the available repositories for data deposition and the cost of storing the comparatively large files associated with digital imaging-based research. Researchers can also be reluctant to share data that remain part of an active research programme [33], or to share a subset of data that is part of a larger, unpublished package. There is also a lack of consensus and widespread confusion over issues of data ownership and copyright, and conflict that emerges between institutional policies asserting copyright ownership (e.g. public museum or even private collections) and the regulations of funding bodies and publishers with regard to open data. Consequently, sharing or publishing supporting data is often a low priority and has effectively been considered optional when not prescribed by a journal. Partial datasets (e.g. low-resolution visualizations or external surfaces) can be insufficient for reproducibility or even verification. As digital morphology has evolved, most of us in the research community have failed to achieve what might now be considered best practice of open data. The academic world has already taken important steps towards overcoming some of these motivational and practical obstacles. Platforms for both archiving and sharing data online are becoming more commonplace, and can handle large file sizes. The standard in molecular biology is GenBank (https://www.ncbi.nlm.nih.gov/genbank/), where sequence data underpinning studies are accessioned before publication. For other data formats, journals and publishers offer a mixed landscape of policies on data publishing that is in need of standardization [34,35], but many not only mandate data deposition—some are even prepared to bear the associated costs, making data deposition easier and ultimately improving science, both in terms of practice and accessibility. There are also initiatives to integrate data submission with submissions to peer-reviewed journals, requiring (or at least allowing) the submission of data in the article submission process and enabling reviewers to examine supporting data as part of the review process [36]. However, collectively, these initiatives have not been integrated [34], and they have not yet translated into common practice within many subdisciplines in biology, palaeontology and anthropology. If a consensus can be established among authors, repositories, journal editors, peer reviewers and funding agencies, there is the prospect of finally realizing the potential of digital morphology in the open-data era. Here, we make recommendations on the nature and extent of essential and recommended best practice datasets that should be made available to support scientific publications using 3D digital datasets across biological sciences (summarized in tables 1 and 2). We review the requirements of associated metadata, discuss the current range of repositories available for such studies and comment on issues affecting their utility.
A range of methods exist for studying 3D specimens through the creation of two-dimensional (2D) image stacks (i.e. tomography), including X-ray CT (encompassing medical CT, micro-CT and synchrotron tomography), MRI, neutron tomography, optical tomography, histological microtomy and physical tomography [1,3,4,37,38]. All of these techniques generate datasets consisting of up to several thousand parallel sections or slices (tomograms) through a specimen, with each tomogram represented by an image file. Various techniques exist for the construction of 3D digital models from sets of tomograms [1]. Image stacks are the starting point for most tomographic studies. These provide immediate insight into internal and external features, and form the basis for any subsequent construction of 3D models. Image stacks exist in a range of non-proprietary file formats, but the most common include DICOM, TIFF, JPEG, PNG, VOL, RAW and BMP [39]. All such files can be opened and viewed in free software such as ImageJ, Drishti, SPIERS, Horos and 3D slicer [40], and can be converted into different formats, although this can be more difficult with DICOM files, which exist in a multitude of sub-formats, not all of which can be handled by all software. For most purposes, TIFFs (16- or 8-bit) provide the best balance of accessibility, file size and data quality (lossless compression), but any lossless, standard image file-types are sufficient. Most JPEG formats enforce a lossy compression scheme that may degrade over multiple save operations; lossless JPEG formats do exist (JPEG-LS, JPEG 2000), but they are not widely used. These differences underlie the importance of specifying the file standard used [39]. Minimally, image stacks should retain the contrast resolution (bit-depth) and spatial resolution used in the study. In cases where the image stack is derived from K-space filling (e.g. MRI) or a series of angular projections (e.g. X-ray CT), the process of generating the image stack is largely automated and we do not consider it necessary to publish the raw projections. An image stack alone will not contain all the information necessary to make full use of the data. For example, scale is only preserved if the resolution (e.g. voxel size or slice spacing) is encoded in the files, and for some datasets slice spacing is not constant and requires per-slice documentation. In the case of DICOMs, this information is typically retained within the file or can be added to the file with a header tag editor (e.g. ImageJ). Otherwise, a text file detailing the voxel or pixel size and slice spacing is the minimum necessary information that must accompany publication of any image stacks. Additionally, metadata information should include full details of how the images were acquired (including scan settings), and further information on data copyright, repository and accession of specimens scanned and, if appropriate, comments on preparation or specimen storage for biological specimens (table 1). This information is necessary to reproduce studies, as well as to evaluate if better-quality data could be obtained with a different set of parameters [41]. Minimally, these data should be provided in a simple text file (e.g. TXT or VGI) associated with the dataset, regardless of whether the information is provided in any study based on the data. Typically, tomographic studies involve the reconstruction of 3D models from image stacks, in some cases after image segmentation or other preparation (see below). 3D models are normally triangle-mesh geometries generated via isosurfacing (usually known as surface models) [1]. Publication of the 3D models resulting from isosurfacing allows for the interactive examination of specimen morphology in three dimensions. A wide range of free software is available for this task [1,3], although no ideal general-purpose file format exists for complex models (see below). 3D models may have been modified after initial isosurface construction, for example through smoothing, island removal or hole filling. Consequently, the most appropriate model to publish to enable verification is the final model (or models) on which the results of the study are based, or which is used in downstream analyses. The 3D models generated using tomographic data are available in a range of different file formats [1,42]. The choice of file type may be influenced by various factors including file size and whether colour/texture information is required; it is essential that openly accessible, standard formats are used (e.g. STL, PLY or OBJ), but there is no single ‘ideal’ file format. The stereolithography (STL) format is the most widely used standard for publishing 3D triangle meshes derived from tomographic techniques, and it is simple and supported by the vast majority of 3D visualization programs, including freely available software [1]. STL files are also compatible with most modern 3D printers, offering potential for wider applications in specimen conservation, public outreach or teaching [3,43]. However, STL files cannot store data on colour, texture or scale. Where these are an essential part of the study, an alternative format such as PLY, OBJ with MTL or VAXML [1,39,42] will be required. These formats are also recommended for meshes with a high number of triangles, which can result in very large file sizes in the STL format. While some tomographic datasets are reconstructed as 3D models without any modification or markup, this is unusual. Most datasets are subjected at least to segmentation, the semi-automated or manual differentiation of voxels (3D pixels) into distinct regions-of-interest (using, for example, ‘label fields’ in Avizo or ‘masks’ in SPIERS). Some datasets also require semi-automated or manual modification of the data (e.g. through brightness modifications) to better separate specimen from background (we term this ‘editing’). These processes involve a degree of subjective interpretation; this is especially true for palaeontological datasets, which are often very noisy and can require extensive manual intervention to extract maximal information from the original data. Thus, publication of the original tomographic dataset and final 3D model may not be sufficient to enable other researchers to assess the association between the two. Segmenting and/or editing a tomographic dataset can be very time-consuming and therefore difficult to reproduce in practice; without access to prepared datasets, most secondary users would not be able to fully interrogate the data underlying a 3D model. In such instances, prepared datasets should be released. No standard file format exists, but labels and masks can be released in the native formats by the software used to generate them, or as binary image stacks, which can then be readily reconstructed as a 3D model in a variety of software packages [1,42]. Development of back-projection algorithms can improve signal to noise ratio in generated image stacks and, hence, recent open-data mandates at synchrotron facilities require archiving of the radiograph projections, not the resulting slice data [44]. Thus, it may be sensible for authors to archive the raw projection libraries themselves. This is especially important where access to the same specimen may be problematic, or as a precaution in case unique specimens are damaged, lost or destroyed. For physically destructive and optical tomography, tomograms need to be registered (aligned relatively and absolutely in the X, Y and Z planes, either manually or semi-automatically) prior to any reconstruction of 3D models. This adds a potentially subjective step that may have a bearing on downstream analyses, and so we recommend publishing both the original (unregistered) and registered image stacks as best practice. Alternative surface-based methods exist for digitizing only the exterior features of specimens in 3D, most notably laser or structured light scanning [45] and photogrammetry [1,46,47]. For photogrammetry, data begin as 2D photographs, whereas in surface-scanning techniques, the 3D shape is usually directly captured as 3D point clouds, with or without texture capture (colour) for each point. In photogrammetry, a 3D polygonal mesh with texture data is generated and warped onto the 3D surface (typically automatically), giving each triangle a colour value. Scanning methodologies may directly visualize point clouds, or may generate and visualize a 3D triangle mesh, with or without texture mapped onto triangles or vertices. The production of the initial 3D surface from photographs or surface scans is largely automated. The most critical data are the final 3D surface files, which may be fused from the original component meshes (e.g. in STL, PLY or OBJ formats) [39]. In cases where the surface texture (i.e. colour information) is directly relevant to the outcomes of a study, the published 3D models must retain this information (i.e. should be provided in PLY or OBJ formats). Surface models are not normally segmented into multiple geometric objects, so single-file models in PLY or STL format are practical. A text file of metadata should be provided that documents details of the imaging settings and techniques used to generate the 3D model (table 1). Preparation of 3D meshes may involve a range of operations, including trimming irrelevant data, realigning or reorienting components of the mesh, fusion into a single mesh, smoothing, hole filling and/or manual manipulation of the location of individual point coordinates or surfaces. These operations should be detailed in the metadata file. Where such operations are non-trivial and/or involve interpretation, those data (photographs, raw point clouds) are an essential provision, in open and widely accessible formats, where possible. Colour data from the surface can provide useful information to help interpret the specimen (e.g. taphonomic preservation). As best practice, this should be included if available, in PLY or OBJ format. The photographs or data captured by the scanner or the 3D data generated by the photogrammetry software allow verification of the processes used to generate the model and should be included as best practice. For 3D scanning, in some cases it may only be feasible to release the raw data in proprietary formats but, where possible, widely compatible (e.g. STL) surfaces should be exported. For methods that involve the digital alignment of different aspects of a specimen, or significant manual intervention in the model construction, the unfused data should be released as the accuracy of the original alignment may be of variable quality. It is important to consider not only the generation of 3D models, but also the data that may be produced in the course of downstream analyses to which these data are subjected. Common types of analysis include: (i) size and shape analyses through topological and landmark-based techniques such as geometric morphometrics; and (ii) assessment of the functional performance of specimens through computer modelling approaches, such as FEA, multibody dynamics analysis (MDA) or CFD. These studies are often based on 3D models with the data subsequently analysed in specialist software packages [1]. For morphometric approaches, the original landmark coordinates and the rules defining landmark location should be provided as these constitute the raw data for the morphometric analyses. For 2D landmark data, a TPS file or similar format links landmarks to their constituent images. Where 3D landmark data points are collected via a 3D digitizer, it is common practice to tabulate the specimen number of the digitized specimen. Where the analyses are based on 3D surfaces or digital models, it is desirable that the models (surface or volume) used in the analysis should be published in an accessible format (following the guidelines outlined above). Functional analyses typically convert 3D digital datasets into proprietary formats for specific methodologies, such as FEA, CFD and MDA. Free software packages do exist, but typically industry standard commercial packages are employed. These have the advantage of reliability and standardized algorithms underpinning the computational analysis. Specialist software has the disadvantage that it outputs data in proprietary file formats that may not be widely accessible to many potential users. For morphometrics, a text file detailing any corrections or transformations applied to the data and an explanation of the analyses should be published. If the morphometric analysis is conducted in the R environment, an annotated R script is a convenient solution. For 3D functional analyses, the (usually proprietary) files containing the analysis set-up and parameters, either with or without the results files, are required for model verification. This addition enables a user with access to the appropriate software to replicate the analyses. Full metadata should be provided with details of processing techniques used to generate the final model, as well as a description of any parameters specified by the user in the analysis (table 1). Analytical techniques used to investigate the function and biomechanical performance of 3D modelled taxa will produce a range of additional digital data, which should also be made available in order to replicate studies. In the case of FEA, programs use volumetric meshes consisting of a finite number of elements. For MDA and CFD, formats such as the parasolid standard are often essential to perform the analyses. Further parameters and boundary conditions are then defined in specialist software (e.g. Abaqus, Ansys, Strand 7, Adams, Opensim, Gaitsym, COMSOL). Ideally, both the model set-up as well as the result files would be published alongside a study. For commercial packages, viewing software is sometimes available which allows the display of models and results files, but no additional analyses. Some industry software packages have text-editor-readable files that list and detail the location and nature of boundary conditions (e.g. INP files for Abaqus FE software). Researchers have a responsibility to ensure that all of the data necessary to reproduce a published study are made available. As explained above, for 3D digital datasets these data may include original 2D images, prepared/segmented 3D images, 3D geometries and relevant metadata. These datasets can be, in toto, very large by today's standards; over 100 GB per specimen is possible in some scenarios, and there may be some instances where single publications utilize huge numbers of specimens, the storage of which is in itself a project. Publishers and other institutions hosting repositories must manage and facilitate access to the data they host, with these obligations persisting into the future, ideally indefinitely. Museums and other institutions holding original specimens often consider digital data as an intrinsic aspect of the specimen, and request researchers to deposit these data with them. Many have active programs of 2D and 3D digital curation, and normally make data freely available for research purposes. Data access for commercial use is a source of much-needed income, and commercial reuse of data released for research purposes is a genuine concern. However, most museums do not yet have systems, policies or resources in place for the long-term curation and distribution of digital morphological data [30]. This is not surprising given the paradigm shift in the concept of the accessioned specimen brought about by digital morphology, expanding from the physical specimen to a diversity of avatars. Digimorph.org pioneered the curation of digital morphological data for in-house scans generated by the University of Texas High-Resolution CT Facility (UTCT), and there are now a number of general and specialist repositories facilitating the publication and dissemination of supporting data at a variety of scales (electronic supplementary material, table S1). Many journals have agreements with such repositories and will cover charges, even for relatively large datasets. In addition, many funding agencies cover the costs of long-term data storage, and many institutions have developed their own data repositories to manage research data generated by their own researchers. Out-moded promises to make data ‘available on request’ should give way to permanent URL links to 3D image data in biology, anthropology and palaeontology (cf. [35]). A range of repositories are available that cater for 3D digital datasets arising from research in biological sciences (electronic supplementary material, table S1). These can vary greatly in terms of the size and types of data they are willing to accept, as well as the cost of storage. In some cases, the choice of repository may be prescribed by the funding body or journal, but this decision will most often be made by the researcher. Modern facilities for publicly sharing datasets include national data centres (typically supported by a research funding body; e.g. RCUK data centres), multidisciplinary (e.g. Dryad, datadryad.org; figshare, figshare.com; MorphoMuseuM, morphomuseum.com; MorphoSource, morphosource.org; Phenome10 K, phenome10 k.org; Zenodo, zenodo.org) or discipline-specific (e.g. XROMM, xromm.org) repositories, and institutional repositories for data produced in-house (e.g. Bristol University's Research Data Repository, data.bris.ac.uk/data; Natural History Museum London's Data Portal, http://data.nhm.ac.uk). It is not entirely clear that all of these are sustainable in the long term. Traditional repositories of physical specimens can also store and disseminate data, and many are moving towards online access to their digital collections. Digital repositories should have the same qualities as repositories of physical specimens, in that they should ensure the long-term persistence and preservation of datasets in their published form, provide expert curation and stable identifiers for submitted datasets, and facilitate public access to data without unnecessary restrictions. However, by their very nature, they should also ensure that the data are discoverable online, provided with unique, permanent and citable reference codes (e.g. DOIs), associated with relevant metadata (e.g. readme text file), and have links to relevant publications and funding bodies [2,28]. The specific licence used by the repository should be considered. Many facilities currently use the CC-BY-NC licence, which disallows reuse for commercial activities. This may be desirable where there are concerns over activities such as selling 3D prints of museum specimens with no benefit to the institutions charged with maintaining those collections. Some data repositories (e.g. MorphoSource) allow users to specify the most appropriate licence for their data. Authors may prefer to choose the CC-BY licence, which is among the most open creative common licences available and has become the standard for open access publication of journal articles. This licence lets others distribute, edit and build upon the original data, even commercially, as long as they credit the original creator. The CC-0 licence (Dryad default) goes further and allows copyright owners to waive all rights. CC-BY-ND is less attractive, as it allows sharing but does not allow the end user to publish derivatives of the data. 3D digital datasets associated with published studies should be verifiable and fully traceable from production to publication, and later republication. One option is digital watermarking, which provides a means of achieving verification of the authenticity and integrity of data, and is imperceptible to the human eye, but also durable in both digital and printed forms, surviving most image edits, file format conversions, data compression, filtering, partial data removal and smoothing. Another option would be to require users to register with the repository before data can be downloaded and used, a practice already imposed by some repositories (e.g. Dryad, MorphoSource). Registration is usually free and open to everyone, but allows the repository to track data access. When publishing large (e.g. more than 10 GB) 3D digital datasets, it is vital to consider the financial costs, which are typically proportional to the amount of data being stored. Some repositories do not currently charge for accessions (e.g. MorphoSource), but for some, accession charges are not insignificant. The popular online digital repository Dryad (datadryad.org) currently charges $120 per data package of 20 GB plus $50 for each additional 10 GB. Datasets based on synchrotron tomography supporting a single publication can easily run to 100 GB for a relatively small number of scans of individual specimens, and it is possible to envisage future projects, especially synthetic papers and large-scale comparative analyses, generating datasets that are orders of magnitude greater in size. Publishing such datasets can quickly become prohibitively expensive; many journals offer to fully or partially cover the costs of depositing digital datasets, but do not have a clear policy for datasets that are hundreds of GB to TB in size. Applications for research funding are increasingly budgeting for data storage costs, but this does not assist projects making use of pre-existing data, or those where funds for data publication are not available. One way of minimizing costs is by reducing the total size of data published without compromising the quality. Cropping of redundant space around a volume representing the specimen is an obvious first step. Lossless compression of individual image files is an excellent route to reduce data storage for image stacks in certain formats. For example, LZW compression, both lossless and fully reversible, can provide upwards of 40% reduction in file size on eight-bit TIFFs with no evident effect on data quality, but it is not routinely applied. The PNG image format provides a similar level of lossless compression. Most of the JPEG image formats enforce lossy compression that degrades data, and should not be used despite appealingly high compression ratios. Placing files into ZIP archives (e.g. one ZIP file per image stack) also reduces disc space through lossless compression and is more convenient for downloading. However, ZIP and VOL archives are less secure for long-term storage since, if the single file containing a dataset becomes corrupted, the entire dataset will be lost. Corruption of single files within a large dataset is less serious, and at least some repositories have procedures in place to detect and remediate bitrot [31]. We recommend that unarchived copies of the original data are stored and made available where possible. In our enthusiasm for recycling 3D digital data and easing reproducibility of morphological studies based on them, the environmental costs of storage should be considered. Most datasets will be accessed infrequently and so there is no need or justification for their storage on spinning discs. Many repositories make use of automated tape storage which is stable and comparatively low in direct costs for the same reasons that make it environmentally low-cost. An increase in the availability and ease of use of data repositories raises the prospect of making data available from previously published studies where the data were not released at the time of publication. Digital datasets can be uploaded to online data repositories and linked to past publications. At present, there are no policies or mechanisms we are aware of among journals and publishing houses to link archival publications to newly deposited data. However, there is no material technical barrier to salvaging legacy data in this way. Publishers are likely to welcome such an initiative as it would obviously improve data visibility, facilitate reproducibility, and probably rejuvenate old publications in terms of access, citations and, ultimately, their marketability. Obtaining digital characterizations of morphology can be time-consuming and expensive, and researchers rarely exhaust their data with the first publication. Funders and publishers are increasingly removing choice over whether to release supporting data, and so it can seem unfair that the researchers who generated datasets have to subsequently compete to exploit them further. This can be particularly difficult for lone early-career researchers potentially competing with large experienced research groups [33]. One potential solution to this would be the introduction of time-limited embargos, which can already be facilitated by some data repositories. However, such embargos violate the most basic tenet of open data: that of removing barriers to assessing the reproducibility of research [48]. After the point of publication, it is also effectively impossible to police the release of supporting data and, consequently, we see no alternative to the release of data with publication. A possible compromise may be borrowed from the Bermuda [49], Fort Lauderdale [50] and Toronto [51] agreements of the genomics community. These mandate data release at the time they are obtained but, more germane to morphologists, these agreements provide safeguarding for data generators through published, time-limited statements of intent of how they propose to exploit the data [51]. Other researchers are free to exploit the data for other purposes, and for any purpose after the stated period of limitation of the statement of intent [52]. Third-party users with overlapping research interests are expected to proceed respectfully and in dialogue with the data generators to identify a mutually agreeable publication schedule [51]. Invariably, much more is at stake in such projects, and though these informal agreements are rarely violated, they are generally well policed by the peer review process [52], and by the reputational damage suffered by those who choose not to observe these agreements. Practice in the genomics community underscores the point that there is more to gain from open data than the warm glow of altruism [51,53]. Not only has it led to greater and more rapid scientific advance [48,51], it can lead to material personal gain, through proposals for collaborative exploitation of published data, both to achieve stated research objectives, and to achieve new objectives that would not be possible without unforeseen collaborators [51,53]. Citation and access-tracking of published datasets also provide credit to the authors [31]. Attribution of authorship is mandated under CC-BY licences and is in any case integral to the academic culture. Many journals already mandate citation of published datasets, not (or not merely) the publications describing research based upon them; this must become common practice. Further mechanisms for encouraging researchers to share their data should only add to this motivation, such as explicitly evaluating the open sharing of data in hiring, promotion or other reward processes. Nevertheless, data can be associated with ethical sensitivities that may require the withholding, or restriction on public distribution, of data (e.g. anthropology or medical science [54,55]). In such instances, the issues that apply should be clearly defined so that beyond these boundaries researchers and publishers can follow an ethos of open-data publication. Mechanisms already exist to cope with these constraints while still making data available, such as data anonymization and vetted access [51]. While the principle of open data has been mandated by the majority of funders [32], publishers, physical repositories and researchers are all scrambling to meet the resulting challenges. Above all, the competing interests over ownership of digital data need to be resolved between (i) funders who pay for research, (ii) researchers who collect specimens and create the digital datasets, (iii) research facilities where data are collected, (iv) museums that have a duty of care for the physical specimens and (v) research publishers. Funders, researchers and publishers may have converged on an ethos of open data. However, the institutions that are responsible for the physical specimens have not obviously been invited to engage in the development of open-data policy, and yet it is museums that will have to change most in terms of their policies on the nature of what they consider intrinsic aspects of the physical specimens that they hold in their care. One solution for museums might be to comply with research funders’ requirements, and waive copyright over digital representations of their collections, along with its associated income stream. Another solution would be for these institutions, which are those best-placed to inform policy on the curation, storage and distribution of data, to develop digital collections with the stability to match that of their physical inventory. Indeed, with the development of cybertypes [28,29], this may be an inevitable future aspect of the world's leading museums. However, if this readily realizable vision of data repository quality, stability and credibility is to be achieved, it will require the funders who have mandated data deposition to cover the costs of establishing and maintaining such facilities, through block grants, not through piecemeal funding to researchers. If such change is to be achieved, it must happen not only in wealthier countries but worldwide, and thus more amply provisioned funders should provide further means to help other countries improve their data-sharing capacities. Data access is not only important post-publication, to aid reproducibility, but during peer review, so that the results of a study and their interpretations can be verified prior to publication. Providing tomographic or 3D data at the point of journal submission is, in our experience, a comparatively rare phenomenon that the publishing infrastructure is not currently well set up to facilitate. Publishers must develop a more homogeneous policy on open data [34], along with procedures to ensure data sources are acknowledged and linked electronically to the derivative publications [48]. It is also important that systems are developed to ease the submission of such data, and facilitate secure, anonymized distribution of data to reviewers. Dryad offers an integrated submission system where publishers can coordinate submission of a manuscript with submission of data, which can then be accessed securely by referees and editors. For non-integrated journals, an interim solution may be to host data at a temporary, hidden URL that can be forwarded to the reviewers via the journal. Authors may be cautious about sharing such data ahead of an article being accepted for publication, and there should be a clear policy governing the restrictions of use for reviewers. Data sharing is essential in order for the benefits of 3D digital data to be fully realized by the scientific community, as well as for the maximum benefit to be gained from the public and private funding that allows these data to be collected. Not only are the benefits of 3D digital data not currently being fully realized, but failure to publish supporting data is rendering many studies based on 3D digital data at least difficult to reproduce. We have presented a series of proposals for open 3D digital data. These outline the minimal standards of verifiability that studies should meet before they are published. We also present more ambitious standards that we hope can be assumed as normal best practice (table 1). We have all been guilty of failing to meet these standards in the past because of technical and other limitations; however, technology has changed and so must we. There are costs associated with releasing data, both real and in-kind, but these are insignificant in proportion to the real costs of regenerating data, and the reputational costs to individuals, institutions, journals and editors of publishing research predicated upon inaccessible data. The project was conceived by T.G.D., I.A.R., S.L., J.A.C., E.J.R. and P.C.J.D., all of whom drafted the original manuscript, to which all others contributed. We declare we have no competing interests. The authors are funded by BBSRC (P.C.J.D., E.J.R.), The Calleva Foundation and the Human Origins Research Fund (C.S.), European Research Council (A.G., J.R.H., R.B.J.B.), Generalitat Valenciana and MINECO (C.M.-.P.), Leverhulme Trust (A.G., R.B.J.B.), NERC (J.A.C., P.C.J.D., A.G., J.R.H., E.J.R.), NWO (M.R.), National Science Foundation (A.G., A.P.S., S.Y.S.), 1851 Royal Commission (I.A.R.), Royal Society Wolfson Merit Award (P.C.J.D.) and the Swedish Research Council (S.B.). We thank Zosia Beckles (data.bris), Else-Marie Friis (NRM, Stockholm), Mark Hahnel (figshare), Iain Hrynaszkiewicz (Springer Nature), Elizabeth Hull (Dryad), Phil Hurst (Royal Society Publishing), Rhiannon Meaden (Royal Society Publishing), Sowmya Swaminathan (Springer Nature), Stuart Taylor (Royal Society Publishing) and Sally Thomas (Palaeontological Association) for discussion. FootnotesElectronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3740174.v1. © 2017 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited. References
Page 13It has long been recognized that the terrestrial tetrapod fauna of the early Permian differs substantially from that of the middle and late Permian. The early Permian fauna is dominated by pelycosaurian-grade synapsids, with the carnivorous Sphenacodontidae and herbivorous Edaphosauridae being the most abundant and diverse large amniotes [1–4]. The eureptile clade Captorhinidae is also prevalent at this time, and there is a diverse array of amphibians [5]. During the middle and late Permian, this fauna is replaced by a therapsid-dominated fauna with therapsid synapsids representing the most common large carnivores (dinocephalians, therocephalians, gorgonopsians) and herbivores (dinocephalians, dicynodonts) [6–9]. Parareptile diversity increases [10] while the amphibian fauna is noticeably reduced [5,11]. These faunal changes were accompanied by a shift towards more complex ecosystems with more trophic levels [2]. A detailed understanding of this transition has been hindered so far by the geographical patchiness of the fossil record. The lower Permian terrestrial record is dominated almost entirely by palaeoequatorial records from North America and Western Europe, whereas the middle and upper Permian record is known primarily from the palaeotemperate localities of Russia and South Africa [12,13]. There is very little overlap between the North American and Russian faunas [14], making it difficult to establish over what time scale the transition took place. It might even be argued that the apparent change in faunas in reality represents the geographical shift in the record; thus, the lower Permian fauna would represent the equatorial fauna while the middle Permian fauna would represent the temperate fauna. Alternatively, it has also been suggested in the past that the faunal transition was accompanied and possibly driven by a mass extinction event dubbed Olson's Extinction [3,15]. Olson's Extinction has been demonstrated in global diversity estimates using sampling correction [4,5,10,16], but even so there is still disagreement. Benson & Upchurch [11] suggested that the apparent mass extinction could be an artefact of the geographical shift in sampling localities, arguing that the sampling of more species-rich low-latitude localities in the Kungurian, and more species-poor higher-latitude localities in the Roadian, could have produced the apparent drop in diversity across the boundary. Fortunately, these two conflicting hypotheses make explicit predictions that can be tested against the fossil record. For the (I) artefact hypothesis (due to sampling of different latitudes) one would expect (i) substantially different faunas in different latitudes even in the early Permian, with faunal composition in high latitude localities more similar to middle Permian localities; and (ii) a latitudinal biodiversity gradient with higher species richness in low latitude localities. For the (II) transition and mass extinction hypothesis, one would expect (i) the lower Permian faunal composition of high latitudes to be more similar to contemporary low latitude faunas than to high latitude middle Permian faunas; and (ii) the decrease in species richness across the early–middle Permian boundary should be apparent even in a single latitudinal bin. Here, using sampling correction techniques [17,18], recent discoveries and up-to-date knowledge of the taxa and biostratigraphy, we can address the following questions. (i) Were there differences between the faunas at different latitudes during the early Permian? (ii) When did the transition between the lower and middle Permian faunas occur, and did it occur simultaneously at different latitudes? (iii) Were there substantial differences in the patterns and changes in species richness at different latitudes? Fortunately, there is temporal overlap between the palaeoequatorial and palaeotemperate formations across the Kungurian–Roadian boundaries, allowing a comparison of the faunas, at least in Laurasia. The Kungurian-aged Inta Assemblage from the Timano-Pechora region of Russia (palaeotemperate Laurasia) can be compared with the considerably better-sampled palaeoequatorial Laurasian fauna of the USA. The Inta fauna has not been extensively discussed in connection with the faunal turnover, despite being the only example available of a lower Permian terrestrial fauna from Russia, as it is an extremely poorly sampled fauna with only 14 specimens representing eight species, and some of these have been lost [19]. Recent examinations of the stratigraphy have supported a correlation between Inta and the faunas of the Vale and Choza formations of Texas and the Hennessey Formation of Oklahoma [14], representing the Laurasian equatorial faunas. The Roadian data also come from Russia and the USA. Although there has been substantial debate surrounding the biostratigraphy of the terrestrial faunas across the lower–middle Permian transition [1,12,14,19–23], it is now widely accepted that the San Angelo Formation of Texas, the Chickasha Formation of Oklahoma and the faunas of the Kazanian series from Russia (the Golyusherma Subassemblage) are contemporaneous and of Roadian age (but see discussion in the electronic supplementary material text). The Cala del Vino Formation of Sardinia is also included as an equatorial Laurasian locality (see electronic supplementary material text), despite only having produced a single tetrapod species: the caseid Alierasaurus ronchii [24,25]. These faunas were used in a quantitative comparison between faunas across the lower–middle Permian boundary in both palaeoequatorial and palaeotemperate latitudes, at least in Laurasia. Unfortunately, there is no information from Gondwana during the Roadian. Recent discoveries in the Pedra de Fogo Formation of northeastern Brazil [26] provide data on a lower Permian palaeoequatorial fauna of Gondwana, but there is no post-transition data until the late Permian (see the electronic supplementary material text for a more detailed discussion of the Pedra de Fogo Fauna). A recent modification of the Forbes similarity index [27], hereafter referred to as Forbes* index, was used to compare the faunas in the four latitudinal bins: palaeotemperate and equatorial in both the Kungurian and Roadian. The Inta fauna (Kungurian palaeotemperate) is found to show greater similarity to the Kungurian equatorial formation than to the Roadian palaeotemperate fauna from Golyurshma (table 1). The difference is not found to be significant, probably due to the small sample size from Inta. Meanwhile Golyusherma group shows greatest similarity with the palaeoequatorial Roadian fauna.
These results suggest that Inta is best interpreted as a pre-transition fauna rather than there being a distinct temperate fauna prior to the Kungurian–Roadian boundary, a result largely driven by the amphibians present. Intasuchus silvicola and Syndyodosuchus tetricus are both archegosauroid temnospondyls [28], and although this clade contains taxa from the USA, Russia and South America throughout the Permian, Ruta et al.'s [28] supertree found both species most closely related to lower Permian taxa from the USA. The other temnospondyl amphibian from Inta, Clamorosaurus nocturnes, is an eryopid, a clade abundant in contemporary formations of the USA (see below) and completely unknown from the middle Permian. Amniotes from Inta are mostly uninformative, the most abundant of which are captorhinids, previously referred to Riabininus uralensis but now recognized as undifferentiated from better-known American forms [29]. Captorhinids are found globally throughout the Permian in both temperate and equatorial regions, so this taxon does not provide a useful point of comparison. Riabininus is a single tooth-rowed captorhinid [29], a group most commonly found in the early Permian of the USA, but not unknown from the upper Permian palaeotemperate latitudes [30]. The only other amniote worth mentioning is an indeterminate form previously assigned to Gnorhimosuchus satpaevi, which has been considered variously a bolosaurid and captorhinid [31,32]. Like captorhinids, bolosaurids are known from both the early and middle Permian of the USA and Russia [33–35], so are also not helpful in this comparison. The amphibians of Inta indicate a fauna similar to that of the equatorial early Permian (i.e. a pre-transition fauna). Unfortunately, there are no fossils of large amniote taxa known (whether due to poor sampling or taphonomy is unclear), but the ichnotaxa Dromopus and Dimetropus are known from the area [36], with the latter traditionally associated with sphenacodontid synapsids [37–39]. Although recent examination shows Dimetropus to be more diverse than previously thought [40], it is still considered representative of large pelycosaurian-grade synapsids, and both Dimetropus and Dromopus are both highly characteristic of the early Permian. While the Kungurian of Russia appears to have contained a pre-transition fauna, the Kazanian (Roadian) aged palaeotemperate fauna of Russia is decidedly middle Permian in character. The most abundant taxa are dinocephalian therapsids (figure 1c), which continued to dominate the terrestrial realm both in South Africa and Russia until the end of the middle Permian [9]. The most abundant taxon in the Golyusherma subassemblage is the herbivorous dinocephalian Parabradysaurus sileantjevi, but also present is the anteosaur Microsyodon orlovi [41] and the enigmatic Kamagorgon ulanovi (of uncertain affinity but considered a dinocephalian by Ivakhnenko [42]). There is a large amount of Kazanian material assigned either to Phthinosaurus borissiaki or to the clade to which it supposedly belongs, the Phthinosuchia [41,43]. This clade is of untested affinity and monophyly, but the most recent opinion on Phthinosuchia assigns it to Dinocephalia [42]. Also abundant in the Kazanian of Russia are the bolosaurid Belebey and the melosaurine archegosaurids Koinia silantjevi and Melosaurus kamaensis (figure 1c). As mentioned above, bolosaurids are known throughout the early and middle Permian in palaeoequatorial and palaeotemperate regions, but Melosaurinae is known only from the middle Permian of Russia [32]. Figure 1. Rank-abundance distributions of the tetrapods from (a) late Kungurian equatorial Laurasian bin, (b) early Roadian equatorial Laurasian bin and (c) early Roadian temperate Laurasian bin. Legend in (a) also applies to (b).
Interestingly, the Roadian marks the separation between the equatorial fauna of the USA and the temperate fauna of Russia. The Roadian palaeotemperate bin might exhibit greater similarity with the contemporary equatorial bin than the earlier temperate bin, but the Roadian palaeoequatorial bin shows the greatest similarity to the Kungurian equatorial bin (table 1). The faunas of the Hennessey, Vale and Choza formations (Oklahoma and Texas) represent the best-sampled assemblages of the late Kungurian. Their fauna is largely similar to that found in the earlier Cisuralian, except for a replacement in the large herbivore guild of Edaphosauridae by another pelycosaur family, Caseidae [1,8,44–46], the most abundant of which is Cotylorhynchus romeri (figure 1a). Edaphosaurus was still present, but is represented only by some neural spine fragments [47] and is clearly not a substantial part of the fauna. As in other lower Permian formations in the USA, the largest terrestrial carnivores are pelycosaurian-grade synapsids from the family Sphenacodontidae (figure 1a), most of which have been assigned to Dimetrodon giganhomogenes. Other abundant, though smaller, terrestrial taxa include the dissorophoid temnospondyl Tersomius and several taxa from the eureptile family Captorhinidae (figure 1a). By far the most common tetrapod taxa in these formations are aquatic amphibians (figure 1a). The most abundant of these are small bodied—the lysorophid Lysorophus tricarinatus, the dvinosaur Trimerorhachis and the nectridean Diplocaulus magnicornis—but a large aquatic carnivore also occurs with high frequency: the eryopid temnospondyl Eryops. The Roadian-aged San Angelo and Chickasha formations represent similar environments to the Kungurian formations from Oklahoma and Texas: fluvial/deltaic [48,49]. Nevertheless, there are differences between the Kungurian and Roadian faunas in the USA. One example is the substantially reduced presence of amphibians (figure 1); a decrease is observed in both abundance and diversity (a fall from 14 species to 5). The few amphibians that are present are from clades characteristic of the early Permian, including a microsaur (Cymatorhiza kittsi), a nectridean (Diplocaulus parvus) and two dissorophoids (Fayalla chicashaensis and Slaugenhopia texensis). It has previously been suggested that the Roadian ecosystems contained a novel food chain, in which the dominant carnivores were varanopid synapsids while the previous top predators, the sphenacodontids, declined to extinction [4,14]. While it is true that varanopids increase in diversity (species richness and morphological diversity) across the Kungurian–Roadian boundary [4], the data presented here show that sphenacodontids remain the most abundant large carnivores in the Roadian of the USA. There are two varanopid species present in the USA at this time, Varanodon agilis and Watongia meieri, with lengths of 1.5 m or more [48,50], but each is represented only by a single specimen. On the other hand, 18 specimens are assigned to Sphenacodontidae (figure 1b), although only one species has been named (Dimetrodon angelensis). The primary consumers in this ecosystem are caseids and captorhinids (figure 1b), again similar to the lower Permian fauna. The only taxon known from Sardinia is a caseid. In summary, the faunal transition observed across the Kungurian–Roadian boundary does not appear to be an artefact of the latitudinal shift in sampling. Quantitative and qualitative examination of the tetrapod fauna of Inta indicates it is more similar to the contemporary equatorial fauna than the later fauna of the same latitude. In Russia, a transition to the therapsid fauna already appears to have occurred by the Roadian (see electronic supplementary material for discussion of Mezen). However, aside from the much-reduced diversity of amphibians, the Roadian of the USA contained much the same set of clades as was present in the Kungurian; it is a pre-transition fauna. Hence, it seems that the transition between the faunas was not consistent in the different latitudinal bins. It is unfortunate that there are no middle Permian equatorial formations aged later than the Roadian, preventing further study of the progress of the transition. It is not until the late Permian that equatorial data are again available, by which time the amniotes present were very much characteristic of the post-transition fauna. The equatorial amphibian fauna, however, retained many characteristics of the Cisuralian fauna as late as the Changhsingian (see electronic supplementary material for a more detailed discussion of the late Permian faunas in different latitudinal bins). The similarity of the amphibians of Moradi and Ikakern to the taxa present in the early Permian has been noted previously, as has their unusual morphology [51], both of which suggest the long isolation of these localities leading to a highly autapomorphic set of species. Sidor [51] suggested that the extreme aridity of the equatorial environments prevented the migration of amphibian taxa into or out of Moradi and Ikakern. The amniotes, not limited by a tie to water, were better able to disperse, leading to more homogeneous amniote faunas. A ubiquitous pattern governing the distribution of modern species is the decreased diversity in higher latitudes relative to lower latitudes, a pattern termed the latitudinal biodiversity gradient (LBG) [52,53]. However, this pattern appears to have only been pervasive since the Neogene, and to have varied greatly in deep time [54–57]. Mannion et al. [58] recently suggested that the modern LBG was only found in periods of glaciation. During hothouse periods when there were no polar icecaps, peak biodiversity was concentrated in temperate latitudes. The early–middle Permian is an important data point when examining the LBG. During this period, the last period of glaciation prior to the Neogene was ending, and the extensive southern icecap present during the Carboniferous and Cisuralian disappeared [59]. The southern tundra appears to have been present until the middle Permian, at which point it was replaced by a cold temperate climate from 30–60°S [60]. The transition between the lower and middle Permian tetrapod faunas coincided with a trend towards warmer and dryer climate [60], providing a perfect opportunity to examine the impact of changing climate on the LBG. A modern LBG was posited by Benson & Upchurch [11] as an explanation for the apparent decrease in global tetrapod biodiversity across the Kungurian–Roadian boundary, dubbed Olson's Extinction (artefact hypothesis, above). They suggested that the observed drop in species richness was an artefact of the increased sampling of the temperate latitudes and reduced sampling of equatorial latitudes; if a modern LBG was present in the Permian, it would be a lower-diversity fauna sampled in the Roadian. However, in all time bins tested, both subsampled diversity (figure 2a) and the diversity calculated using TRiPS (electronic supplementary material, data S16) show a reduced diversity in equatorial latitudes relative to temperate latitudes. The gradient is weak during the Kungurian but strengthens during the middle and late Permian. In Laurasia during the Roadian and Wuchiapingian, the diversity in the temperate latitudes was more than double that of the equatorial latitudes. During the Changhsingian the gradient seems to strengthen still further. Given the continual trend towards increased temperature throughout the Permian [61], the data presented here support the theory of Mannion et al. [58]: the trend towards a hothouse climate shifts the diversity peak towards temperate latitudes, leading to the establishment of an inverted LBG already in the early Permian. This aspect of the artefact hypothesis proposed by Benson & Upchurch [11]—that Olson's Extinction was an artefact of different sampling of the different latitudes—requires a modern LBG, whereas it is shown here that an inverse gradient was prevalent throughout the Permian. The shift in sampling from palaeoequatorial to palaeotemperate localities was a shift to more species-rich localities, and therefore cannot be used to explain the observed drop in diversity. Figure 2. Subsampled diversity curves. (a) Comparing latitudinal bins and (b) comparison through time.
Nevertheless, our detailed examination of individual faunas shows that Olson's Extinction was a more complicated event than has been previously suggested. The strengthening of the LBG across the Kungurian–Roadian boundary appears to have been caused by a sharp drop in species richness in the palaeoequatorial latitudes; that of the temperate faunas shows no decrease (figure 2; electronic supplementary material, data S16). The equatorial and temperate faunas appear to have responded differently to the climatic changes occurring across the Kungurian–Roadian boundary. The low-latitude extinctions may be due to the disappearance of the everwet biome prevalent since the Carboniferous [13]. The formations in the temperate latitude exhibit a complete faunal turnover across this boundary, but the species richness of the pre- and post-transition faunas was largely unchanged. On the other hand, the equatorial fauna remained unchanged in composition until later in the Permian, but the species richness dropped noticeably (figure 2). While the overall decline in global tetrapod biodiversity across the Kungurian–Roadian boundary has been borne out by sampling correction [4,5,10,11,16], in temperate latitudes the turnover appears to have been rapid, and a new equally diverse fauna replaced the old by the early Roadian. In the equatorial latitudes, replacement was slow; during the early Roadian the tetrapod fauna was merely a depauperate version of the Kungurian fauna. Equatorial species richness does not appear to have recovered even by the late Permian. With an up-to-date biostratigraphy and modern quantitative methods, it is possible to empirically test the two hypotheses explaining the apparent faunal shift during the middle Permian: (i) the artefact hypothesis and (ii) the transition and mass extinction hypothesis. As discussed in the introduction, each hypothesis makes two predictions, one regarding the relative similarity of the faunas in different latitudinal bins (tested using the Forbes* similarity metric) and one regarding the relative species richness in the different latitudinal bins (tested using sampling corrected diversity estimates). In the case of the artefact hypothesis, neither prediction is borne out: (i) the lower Permian faunas at different latitudes are not substantially different, and in fact the Kungurian temperate faunas are more similar to contemporary equatorial assemblages than the Roadian temperate faunas; and (ii) the diversity estimates do not show higher species richness at low palaeolatitudes, but instead in the palaeotemperate latitudes. While it therefore appears that the faunal transition is not an artefact of sampling biases, the event does appear to have been more complex than previously considered. The faunal transition was rapid at higher latitudes, with the middle Permian fauna established by the Roadian, while conversely at lower latitudes the transition appears to have been slower, with equatorial Roadian faunas still containing a similar assemblage of taxa to that present in the Kungurian, albeit less diverse. A post-transition equatorial fauna is not seen until the late Permian (although the lack of equatorial samples from the Wordian and Capitanian should be noted). Olson's Extinction has been a controversial event since it was first described, but detailed examination of the faunas shows it was an important episode in the early evolution of terrestrial ecosystems. Although hindered by a geographically patchy record, the data are available to reveal the pattern of the mass extinction, including taxonomic selectivity, geographical variations and the establishment of the post-extinction fauna. This study focuses on the four time bins indicated in electronic supplementary material, table S1, selected as the only time periods when terrestrial specimens are known from localities from both temporal and equatorial palaeolatitudes. The localities in each time bin were assigned to one of four latitudinal bins: equatorial Laurasia, equatorial Gondwana, temperate Laurasia and temperate Gondwana. The boundary between equatorial and temperate is set at palaeolatitude 25° north and south. The formations present in each time and latitudinal bin are shown in electronic supplementary material, table S1. A database of the specimens found in each latitudinal and time bin and their taxonomic assignments was collated from a variety of sources, including the published literature, the Paleobiology Database [62] downloaded via the Fossilworks website (http://fossilworks.org/), museum catalogues and direct observation. The data itself and full details on the data's provenance and treatment are provided in the electronic supplementary material text and data. These data were used to calculate a diversity (species richness) estimate for each latitudinal and time bin. Due to the high heterogeneity of sampling during the Permian [4,6,7,11,63], sampling correction was implemented using the Shareholder Quorum Subsampling method (SQS) [17] and TRiPS [18]. The former method was implemented in R version 3.1.2 [64] using version 3.3 of the script available on the website of John Alroy (http://bio.mq.edu.au/~jalroy/SQS.html). In each bin, diversity was calculated at six quorum levels, from 0.4 to 0.9 at intervals of 0.1 (note that diversity cannot be calculated at all quorum levels in all samples; the maximum quorum possible is the Good's U-value of the sample). 10 000 subsampling trials were carried out at each quorum level. Following recommendations of Alroy [17], the most abundant species in each sample was not included when calculating the proportion of shares contributing to the quorum, although it was included in diversity counts if drawn. TRiPS is an approach that estimates richness and confidence intervals by modelling sampling within each time interval as a homogeneous Poisson process [18]. This analysis was performed in R using the script provided by Starrfeltt & Liow [18]. The results shown in the paper are using SQS, while those of TRiPS are presented in the electronic supplementary material, data S16. Pairwise analyses of faunal similarity were carried out on the Laurasian latitudinal bins across the Kungurian–Roadian boundary using the recently published Forbes* similarity measure [27]. This is based on an older metric [65] that ranges from 0 to 1, with 0 meaning no taxa are shared between two assemblages and 1 meaning the assemblages are identical or one is a subset of the other. Alroy [27] modified this metric to make it suitable for palaeontological data, where sampling is incomplete and the total number of taxa present was unknown. Since, with very few exceptions, the species present are endemic to each formation, the similarity of the faunas must be assessed based on sharing taxa with close phylogenetic relationships. Unfortunately, the available distance/similarity metrics do not yet incorporate phylogenetic information, and many of the Russian taxa have yet to be incorporated into a phylogenetic analysis. Until these issues are resolved, we follow the recent practice of Sidor [51] of performing the similarity analysis at the family level (or, where family-level taxonomy is not so well established, e.g. within Dinocephalia, at higher taxonomic levels). The Forbes* metric was calculated using a custom R script (see electronic supplementary material) written from the equations presented by Alroy [27]. To provide a measure of the statistical significance of the similarity/distance between the assemblages, 1000 simulated datasets were constructed, where assemblages with the same sample size as the observed data had their taxa randomly assigned to families. These simulated datasets were compared pairwise with the Forbes* metric and the results were compared with the observed data. Fossil occurrence data are uploaded as online electronic supplementary material. N.B., M.O.D. and B.S.R collected data. N.B. analysed data. N.B., M.O.D., B.S.R. and J.F. wrote the paper. The authors declare no competing interests. This study was funded by a Deutsche Forschungsgemeinschaft grant (number FR 2457/5-1), awarded to J.F., the DST/NRF Centre of Excellence in Palaeosciences, the NRF African Origins Platform, and the Scatterlings of Africa programmes of the Palaeontological Scientific Trust (PAST). We would like to thank Roger Benson and Roger Close, as well as the Fröbisch working group, for helpful comments and discussion. Alexander Dunhill and an anonymous reviewer provided many helpful suggestions on how to improve an early draft of the manuscript. This is Paleobiology Database Official Publication number 278. We are grateful to John Alroy and Johannes Müller for uploading the relevant data to the database via the Fossilworks platform. FootnotesElectronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3741413. References
Page 14You have access Proc. R. Soc. B280, 20122495 (28 November 2012; Published online 22 January 2013) (doi:10.1098/rspb.2012.2495) There were two minuses missing in the Results section of the original article. The correlation between circulating testosterone levels and body adiposity is negative (r69 = −0.54, p < 0.001), not positive as previously indicated. Similarly, the correlation between circulating testosterone levels and facial adiposity is negative (r69 = −0.52, p < 0.001), not positive as previously indicated. These results indicate that men with higher body and facial adiposity have lower testosterone levels, not higher testosterone levels as originally indicated. Page 15You have accessCorrection Proc. R. Soc. B284, 20170125. (8 March 2017; Published online 15 March 2017) (doi:10.1098/rspb.2017.0125) We recently discovered an error in the section summarizing author contributions. Below, we provide a corrected version of author contributions. I.K. and J.R. conceived of and designed experiments; I.K. carried out experiments; I.K. analysed data with inputs from J.R. and D.A.; I.K., J.R. and D.A. wrote the manuscript. All authors gave final approval for publication. FootnotesPage 16
You have to enable JavaScript in your browser's settings in order to use the eReader. Or try downloading the content offline DOWNLOAD |
Why is sympatric speciation less likely to occur than allopatric speciation?
Postagens relacionadas
direito autoral © 2024 cemle Inc.
