Beard, K. C. Gliding Behavior and palaeoecology of the alleged primate family Paromomyidae (Mammalia, Dermoptera). Nature 345, 340-341 (1990). Show Beard, K. C. The Hunt for the Dawn Monkey: Unearthing the Origins of Monkeys, Apes, and Humans. Los Angeles: University of California Press, (2004). Bloch, J. I. & Boyer, D. M. Grasping Primate Origins. Science 298, 1606-1610 (2002). Bloch, J. I., Silcox M. T., et al. New Paleocene skeletons and the relationship of plesiadapiforms to crown-clade primates. Proceedings of the National Academy of Science 104, 1159-1164 (2007). Carlsson, A. Uber die Tupaiidae und ihre Beziehungen zu den Insectivora und den Prosimiae. Acta Zoologica 3, 227-270 (1922). Cartmill, M. Rethinking Primate Origins. Science 184, 436-443 (1974). Clemens, W. A. Purgatorius (Plesiadapiformes, Primates?, Mammalia), a Paleocene immigrant into Northeastern Montana: stratigraphic occurrences and incisor proportions. In Fanfare for an Uncommon Paleontologist: Papers in Honor of Malcolm C. McKenna. eds. Dawson, M. R. & Lillegraven J. A. (Pittsburgh: Bulletin of the Carnegie Museum of Natural History, 36, 2004) 3-13. Eriksson, O., Friis, E. M. et al. Seed size, fruit size, and dispersal systems in angiosperms from the Early Cretaceous to the Late Tertiary. American Naturalist 156, 47-58 (2000). Fu, J.-F., Wang, J.-W. et al. The new discovery of the Plesiadapiformes from the early Eocene of Wutu Basin, Shandong Province. Vertebrata PalAsiatica 40, 219-227 (2002). Ho, S. The molecular clock and estimating species divergence. Nature Education 1 (2008) http://www.nature.com/scitable/topicpage/the-molecular-clock-and-estimating-species-divergence-41971 Hoffstetter, R. Phylogénie des Primates. Bulletins et Mémoires de la Société d'Anthropologie de Paris t4, XIII, 327-346 (1977). Janečka, J. E., Miller, W., et al. Molecular and genomic data identify the closest living relative of Primates. Science 318, 792-794 (2007). Kay R. F. & Cartmill M. Cranial morphology and adaptations of Palaechthon nacimienti and other Paromomyidae (Plesiadapoidea, ?Primates), with a description of a new genus and species. Journal of Human Evolution 6, 19-53 (1977). Kay, R. F., Thorington, R. W. et al. Eocene plesiadapiform shows affinities with flying lemurs not Primates. Nature 345, 342-344 (1990). Luckett, W. P. Comparative Biology and Evolutionary Relationships of Tree Shrews. New York: Plenum Press (1980). MacPhee, R. D. E., Cartmill, M. et al. Craniodental morphology and relationships of the supposed Eocene dermopteran Plagiomene (Mammalia). Journal of Vertebrate Paleontology 9, 329-49 (1989). Martin, R. D. Towards a new definition of Primates. Man 3, 377-401 (1968). Novacek, M. J. Mammalian phylogeny: shaking the tree. Nature 356, 121-125 (1992). Rose, K.D. , Chew, A.E., et al. Earliest Eocene mammalian fauna from the Paleocene-Eocene thermal maximum at Sand Creek Divide, Southern Bighorn Basin, Wyoming. University of Michigan Papers on Paleontology 36, 1-122. (2012) Sargis, E.J. New views on tree shrews: the role of tupaiids in primate supraordinal relationships. Evolutionary Anthropology 13, 56-66 (2004). Silcox, M. T. New discoveries on the middle ear anatomy of Ignacius graybullianus (Paromomyidae, Primates) from ultra high resolution X-ray computed tomography. Journal of Human Evolution 44, 73-86 (2003). Silcox, M. T. The biogeographic origins of Primates and Euprimates: East, West, North, or South of Eden? In Mammalian Evolutionary Morphology: a Tribute to Frederick S. Szalay. eds. Dagosto M. & Sargis, E. J. (New York: Springer-Verlag, 2008) 199-231. Silcox, M. T. & Gunnell, G. F. Plesiadapiformes. In Evolution of Tertiary Mammals of North America Vol. 2: Marine Mammals and Smaller Terrestrial Mammals. Eds. Janis C. M., Gunnell, G. F. & Uhen, M. D. (Cambridge: Cambridge University Press, 2008) 207-238. Silcox M. T., Bloch J. I., et al. Euarchonta. In The Rise of Placental Mammals: Origins and Relationships of the Major Extant Clades. Eds. Rose K.D. & Archibald J. D. (Baltimore: Johns Hopkins University Press 2005) 127-144. Silcox M. T., Bloch J. I., et al. Cranial anatomy of Paleocene and Eocene Labidolemur kayi (Mammalia: Apatotheria) and the relationships of the Apatemyidae to other mammals. Zoological Journal of the Linnean Society 160, 773-825 (2010). Soligo, C. & Martin, R.D. The first primates: a reply to Silcox et al. Journal of Human Evolution 53, 325-328 (2007). Springer M. S., Murphy W. J., et al. Placental mammal diversification and the Cretaceous-Tertiary boundary. Proceedings of the National Academy of Science 100, 1056-1061 (2003). Springer M. S., Stanhope M. J. et al. Molecules consolidate the placental mammal tree. Trends in Ecology and Evolution 19, 430-438 (2004). Stauffer, R.L., Walker, A. et al. Human and ape molecular clocks and constraints on paleontological hypotheses. The Journal of Heredity 92,469-474 (2001). Sussman R.W. Primate origins and the evolution of angiosperms. American Journal of Primatology 23, 209-223 (1991). Sussman R. W. & Raven P. H. Pollination of flowering plants by lemurs and marsupials: a surviving archaic coevolutionary system. Science 200, 731-736 (1978). Szalay, F. S. The beginnings of primates. Evolution 22, 19-36 (1968). Szalay F. S. & Drawhorn G. Evolution and diversification of the Archonta in an arboreal milieu. In Comparative Biology and Evolutionary Relationships of Tree Shrews. ed. Luckett , W. P. (New York: Plenum Press ,1980) 133-169. Van Valen L. M. & Sloan R. E. The earliest Primates. Science 150, 743-745 (1965).Tropical angiosperms have been a major source of nutrition for primates throughout their evolutionary history [1]. Almost all of the more than 500 extant primate species [2] are highly dependent upon the fruits, leaves, seeds, flowers, nectar, bark, and other parts of flowering plants to meet their energetic and nutritional demands for survival, growth, development and reproduction [3–5]. Consequently, primates may have co-adapted with angiosperms, particularly with those plant families they have interacted with most intensely over time [6]. This dependence and apparent competition for angiosperm food resources has likely also shaped phylogenetic patterns in diet across the primate family, where even closely related lineages may vary greatly in dietary specialization (e.g. obligately frugivorous chimpanzees versus opportunistically frugivorous gorillas). Thus, the evolution of primates is intrinsically interconnected with the evolution of angiosperms over the last 55 million years [7–9], and the mutualistic interactions (e.g. seed dispersal) between primates and plants may have even shaped the diversification of different primate lineages through time [10]. The nutritional composition and mechanical properties of different plant parts represents a crucial aspect of the consumption of angiosperms by primates, profoundly shaping primate ecology, morphology, physiology and behaviour. Fruits are important sources of sugar, whereas leaves are major sources of protein [1,3]. As a consequence, more frugivorous primates tend to target daily protein thresholds, while more folivorous species tend to target daily caloric thresholds [11,12]. As fruit resources may be spatially and temporally patchy, more frugivorous species tend to have larger home range sizes compared to more folivorous species [1]. Furthermore, primates with more folivorous diets tend to have a larger body size to cope with the high fibre content of their diets, given their large guts that allow them to more effectively digest their cellulose-rich diets and more efficient metabolism to subsist on lower quality foods [13,14]. While much is known about the role of plant parts and nutritional content from local field studies of primate dietary ecology, no large-scale synthesis of the taxonomic composition of primate diets is available. Studying the relative proportion of different plant families in primate diets may reveal important insights into how variation in plant morphology and chemistry, including nutrients, toxins and other secondary compounds, across plant species, genera and families relates to variation in primate physiological, morphological and behavioural traits. For example, the leaves and seeds of Fabaceae tend to be higher in protein content compared to other plant families due to the ability of most species in the family to fix nitrogen [15]. Such phylogenetic conservatism in traits at the plant family level [16] similarly extends to toxic secondary compounds (e.g. tannins, alkaloids), which may require special adaptations that allow primates to circumvent these defences. Primate traits may also be shaped by the temporal and spatial availability of resources provided by different plant families (e.g. phenology of fruiting, size of fruit crop). Here, we compiled a novel global dataset of wild primate diets to synthesize our knowledge of primate's use of plant families and their parts (e.g. fruit, leaf). For widely consumed plant families, we tested for phylogenetic patterns in consumption across primate taxa and to what extent the consumption of these families was associated with primate dietary specialization, home range size and body mass. We hypothesized that (i) only a few plant families would form a core component of primate diets, while most others would be consumed to a small degree, with different plant families used for different parts, given variation in availability and phytochemistry; (ii) a phylogenetic signal would exist in plant family consumption given dietary variation across primate lineages (e.g. folivorous colobines will consume high amounts of legumes given their high protein content of leaves) and (iii) the most consumed plant families would influence primate traits such as body mass and home range size given their importance in terms of availability and phytochemistry. Although other primate dietary databases exist (e.g. [17]), we are not aware of any open-access dataset that covers the full breadth of the primate phylogeny with specific descriptions of food items for each plant family and their parts. To evaluate the relative contribution of different plant families across wild primate diets, we compiled quantitative studies of diet composition across primate species using a structured literature search. We obtained a preliminary set of articles using Web of Science in May 2019 using the following search terms: (Food* OR Feed* OR Diet*). We started with articles that were published in four well-known primate journals (i.e. American Journal of Primatology, International Journal of Primatology, Folia Primatologica, and Primates), resulting in 2091 candidate articles. This initial search resulted in dietary studies for 71 primate species, which we further supplemented with references from a recently compiled primate trait database paper [17], as well as an expanded search to increase the taxonomic coverage of our database. The title and abstracts of articles were then screened using the R package ‘metagear’ [18]. Only articles that contained dietary studies with a quantitative and taxonomic breakdown of plant food items were considered. We defined a study as any quantitative summary of the diet based on direct observations of a population or group of a primate species over a discrete time period. Our criterion thus excluded studies that relied on indirect measures such as faecal sampling or metabarcoding approaches. When data from multiple populations or groups of a primate species or multiple species were reported in a single reference, we treated them as separate dietary studies. All dietary records for plant and non-plant food items reported in each study were compiled. For dietary records of plant food items, we additionally recorded the taxonomic identity (plant genus and species names) and plant part (e.g. leaves, fruits) if provided (e.g. ‘Ficus natalensis young leaves'). Each dietary record was defined as the amount of time spent feeding on a given food item as a percentage of the overall amount of time observed feeding, the number of times a population was recorded feeding on a given dietary item as a percentage of the total number of feeding events observed, or an estimate of the relative weight of consumption compared to all items consumed, depending on the methodology of each study (e.g. [4]; electronic supplementary material, table S1). Diets were typically quantified either using scan or continuous sampling techniques (e.g. [19]). Our dataset consisted of dietary information on 119 primate species, derived from 9270 dietary records across 232 studies in 149 references (electronic supplementary material, table S1). However, we omitted studies for which (i) the sampling period was shorter than six months in duration (to capture seasonal variation), (ii) the majority of the diet (greater than 50%) was not known or reported or (iii) the plant family for the majority of the plant portion of the diet (greater than 50%) was either not identified by the author(s) of the study or could not be taxonomically verified by us during collation. The filtered dataset comprised dietary records of 112 primate species across 220 studies from 141 references. Studies of primate species were geographically extensive and were distributed across Central and South America, the African continent, Madagascar and Asia (electronic supplementary material, figure S1). The median duration of studies in the filtered dataset was 12 months (s.d. = 8.7, range = 6–66). For the majority of primate species, dietary information was derived from only one study (s.d. = 1.3; 60 out of 112 species). For the remaining, dietary information of primate species was derived from between two to eight separate studies. Of the 8981 dietary records in the filtered dataset, 8744 are associated with plants. Of these, 84% contained information on the plant part consumed. Of those dietary records (n = 7327), a total of 45 unique plant part terms were reported, sometimes in combination with many records (e.g. fruits and leaves). The large majority of plant parts identified were associated with leaves (n = 3198) and fruits (n = 3295), although a large number of records (n = 2913) was associated with other plant parts, including seeds and flowers (electronic supplementary material, table S2). Taxonomic identity for each plant dietary item was standardized against Plants of the World Online [20]. To quantify the relative importance of each of the plant families, we calculated the proportion of primate species that consumed each plant family. For all families consumed by more than 50% of primate species, we additionally quantified the percentages of each family (regardless of plant part) in the diet of each primate species (averaged across multiple studies of the same species) and the relative consumption of different plant parts (leaf, fruit and other) for each plant family (for more methodological details, see supplementary material and electronic supplementary material, figure S2). To evaluate if closely related primate species would feed on a given plant family to a similar degree, we tested if the consumption of the widely consumed plant families (by more than 50% of primate species) showed a phylogenetic signal across the primate phylogeny. We analysed the average dietary percentages (regardless of plant part) for these families by calculating Pagel's λ [21] and Blomberg's K [22], two of the most widely used measures of the phylogenetic signal [23,24]. For both Pagel's λ and Blomberg's K, values of 0 indicate that the degree of consumption in a given plant family was random with respect to phylogeny, whereas values at or above 1 indicate that consumption was at or greater than expected under a Brownian motion model. We additionally analysed the occurrence of each of the widely consumed families in the diet (yes or no) using a binary measure of phylogenetic signal, D, which is based on a Brownian threshold model [25]. This allowed us to test if consumption of any given family was more likely within closely related primates, regardless of the degree of consumption actually observed. To perform our tests of phylogenetic conservatism, we generated the maximum clade credibility tree using the posterior distribution of phylogenies provided by [26], using TreeAnnotator v. 1.10.4 [27]. Species names in our database were first reconciled with tip labels of the tree using the Integrated Taxonomic Information System database (www.itis.gov) and the IUCN (IUCN v. 3, 2019). In total, 108 out of the 112 primate species in the filtered dataset were represented in the phylogenetic tree. Pagel's λ and Blomberg's K, including significance tests, were performed using the ‘phylosig’ function of the ‘phytools’ R package v. 0.6-99 [28]. D and significance tests were performed using the ‘phylo.d’ function of the ‘caper’ R package v. 1.0.1 [29]. For each test, correction for multiple testing [30] was performed using the ‘p.adjust’ function in R. We estimated the degree of frugivory and folivory for each primate species using the same approach used for calculating the relative proportion of different plant part categories for each plant family: we calculated the average dietary percentages for each plant part category for each primate species (n = 97; see electronic supplementary material, figure S3). Supporting the reliability of our dataset and approach, our estimates of frugivory and folivory were highly congruent with two widely used trait databases that semi-quantitatively score the degree of frugivory and folivory from the natural history literature, EltonTraits [31] and MammalDiet [32] (see electronic supplementary material, figures S4 and S5). We used structural equation models (SEM) [33] to evaluate potential causal relationships between the degree of frugivory or folivory, home range size and body mass. This was performed for 92 primate species for which estimates of both home range size and body mass were available [17]. Body mass estimates were an average of both sexes. Home range size and body mass were log-transformed. Where multiple estimates of either home range size and body mass were available, the average was used. Because we were interested in the role of the most important plant families in mediating these relationships, we additionally quantified the taxonomically identified proportion of either the leaf or fruit diet that was made up by the most consumed plant family for these two dietary categories (family percentage). We created two SEMs, one for frugivory and one for folivory, and both included three main components: (i) home range size as a response variable, with body mass, degree of frugivory or folivory, and family percentage as dependent variables, (ii) body mass as a response variable and degree of frugivory or folivory as well as a family percentage as dependent variables and (iii) a hypothesized path where family percentage was a function of the overall degree of folivory or frugivory. Because relationships may be driven by phylogenetic conservatism in traits (e.g. body mass, degree of frugivory/folivory due to dietary adaptations), the three SEM components (corresponding to linear regression models) were alternatively fitted using phylogenetic generalized least-squares (PGLS) models using the ‘gls’ function from the ‘nlme’ R package (v. 3.1-140, [34]) with the correlation argument set using the ‘corBrownian’ function from the ‘ape’ R package v. 5.3, with the ‘gamma’ parameter set to 1. This approach assumes that the correlation between points is consistent with a Brownian model on the phylogenetic tree. We omitted two species in our phylogenetic SEM as they were not in the phylogenetic tree. Non-phylogenetic SEMs were implemented using the ‘lavaan’ R package v. 0.6-5 [35]. Phylogenetic SEMs were constructed from component PGLS models using the piecewiseSEM R package v. 2.1.0 [36]. Across the 112 primate species from the filtered database covering 141 publications on primate diets, 205 angiosperm families were recorded (figure 1a), out of the approximately 416 angiosperm families worldwide [37]. Plant parts consumed included fruits, leaves, seeds, flowers and other items (e.g. exudates, tubers, bark) (electronic supplementary material, table S2). However, most of these 205 plant families were consumed by a limited number of primate species, with only 10 families being consumed by more than half of all species (figure 1b). Most families only had one or two genera (and four species or fewer) recorded. In general, only a few genera made up the most species recorded in each family. The three genera with the greatest number of species consumed were Ficus (Moraceae, 114 species), Pouteria (Sapotaceae, 40 species) and Inga (Fabaceae, 38 species) (electronic supplementary material, figure S8). Moraceae and Fabaceae were the two most consumed plant families in terms of both the number of primate species documented to consume at least one food item from that family and median percentage of diet across all diets. Out of the 112 primate species in our final database, 91 were observed feeding on Moraceae (81.3%) and 88 on Fabaceae (78.6%), while the median percentage of Moraceae and Fabaceae across primate diets was 7% and 6%, respectively (electronic supplementary material, table S3). Consumption of different plant families was also associated with different plant parts. Fruits of Rubiaceae, Moraceae, Sapotaceae, Myrtaceae and Anacardiaceae were consumed to a greater degree than their leaves (ratio of fruit : leaf consumed greater than 1.5), whereas only leaves of Fabaceae were consumed more than fruits (ratio of leaf : fruit consumed greater than 1.5) (figure 1c).  Figure 1. Relative consumption of plant families across primate diets. A total of 205 plant families were recorded across the diets of 112 primate species. (a) Proportion of primate species consuming each plant family with the 10 families that were consumed by more than 50% of primate species highlighted in red. (b) Proportion of primate diets that came from the 10 families that were consumed by more than 50% of primate species regardless of plant part (displayed in descending order). Bold black lines represent the median percentage across all primate species, boxes represent the interquartile range and points represent species that were more than 1.5 times the interquartile range. (c) Relative proportion of different plant parts for the 10 families that were consumed by more than 50% of primate species (displayed in descending order), averaged across the diets of primate species recorded to consume them. ‘Others’ represents dietary records from the ‘seeds’, ‘flowers’, and ‘others’ categories (electronic supplementary material, table S2). (Online version in colour.)
The degree of consumption was random with respect to the phylogeny of primates for plant families consumed by at least half of the primate species (electronic supplementary material, table S3). A phylogenetic signal was also not detected when we considered the occurrence of consumption (binary variable; ‘yes’ or ‘no’) for most plant families. However, there was one exception. The occurrence of Fabaceae (D = 0.57, p = 0.03) consumption showed some non-random phylogenetic signal (i.e. D < 1), although such a signal was fairly weak and not as strong as that expected under a Brownian threshold process (i.e. D = 0). Across primates with data on plant parts consumed (n = 97), fruits were the most consumed item (median = 36%, s.d. = 25%, range = 0–92%), followed by leaves (median = 17%, s.d. = 27%, range = 0–92%), others (e.g. seeds, flowers) (median = 14%, s.d. = 19%, range = 0–88%) and animals (median = 0%, s.d. = 10%, range = 0–59%). A median of 5% of primate species diet was unidentified or not reported (s.d. = 12%, range = 0–50%; electronic supplementary material, figure S3). For our analysis of the potential influence of plant families on primate traits, we focused on Moraceae fruits and Fabaceae leaves (figure 2b,c), as they were identified as the two most widely consumed plant families that formed a disproportionately large proportion of the fruit and leaf diets of primates, respectively. In our SEM, we found that frugivory primarily influenced home range size indirectly through a negative effect on body mass (path coefficient = −0.344, figure 3a,b). The more frugivorous a primate species, the smaller their body mass and thus the smaller their home ranges. The proportion of Moraceae in the fruit diet had a significant direct negative effect on home range size (path coefficient = −0.253), suggesting that with the degree of frugivory being equal, primate species that consume a greater proportion of Moraceae fruits are associated with smaller home ranges, independent of body mass. After correcting for phylogeny, body mass was still correlated with home range size, but the effect of frugivory on body mass (path coefficient = −0.109) and the effect of Moraceae on range size (path coefficient = −0.125) were no longer statistically significant (figure 3b). However, under this model, frugivory had a significant positive correlation with Moraceae fruit consumption (path coefficient = 0.384). Figure 2. Relative consumption of Moraceae and Fabaceae across primate diets. (a) Average proportion of Moraceae (inner ring) and Fabaceae (outer ring) in the diets of 108 primate species depicted as a heatmap with darker colours representing a higher proportion. (b) Common chimpanzee (Pan troglodytes) feeding on fig (Ficus) fruits (Moraceae). (c) Red colobus monkey (Piliocolobus tephrosceles) feeding on legume (Albizia) leaves (Fabaceae). Primate silhouettes obtained from PhyloPic (phylopic.org) and are in the public domain, except for the image of Lemur by Roberto Díaz Sibaja under a CC BY 3.0 Unported license (https://creativecommons.org/licenses/by/3.0/) and image of Propithecus by Terpsichores under a CC BY-SA 3.0 Unported License (https://creativecommons.org/licenses/by-sa/3.0/). Photo credit (b,c) to Julie Kearney Wasserman. (Online version in colour.)
Figure 3. Structural equation models depicting relationships between primate diets and traits. (a) Phylogenetic distribution of primate trait variation and diet depicted by a heat map. Colours represent standardized values (z-values), defined as the number of standard deviations that a value is away from the mean value of each trait (dark red indicating higher than average values whereas dark blue indicated lower than average values). (b) Structural equation model illustrating the relationships between primate body mass, home range size, degree of frugivory and % Moraceae fruits in the diet. (c) Structural equation model illustrating the relationships between primate body mass, home range size, degree of folivory and % Fabaceae leaves in the diet. Arrows represent direct effects with standardized path coefficients from non-phylogenetic SEMs above and PGLS models in parentheses below. Statistically significant relationships are marked with an asterisk (*) with red representing a negative relationship and green a positive relationship. For each response variable, the amount of explained variance (R2) is given. (Online version in colour.)
Folivory was positively associated with larger body mass (path coefficient = 0.409) but negatively associated with home range size (path coefficient = −0.446), while body mass was associated with larger ranges (path coefficient = 0.668, figure 3a,c). This meant that body size kept equal, the more folivorous a primate, the smaller its home range size. However, the effect of folivory on body mass was not significant in the phylogenetic analysis (path coefficient = 0.076). The proportion of Fabaceae leaves in the leaf diet had no effect on either body mass or home range size in either model (figure 3c). Using an unprecedented dataset covering the relative contribution of all plant families in the diets of more than 100 species of primate worldwide, we show that primates consume a wide range of dietary items, including fruits, leaves, seeds, exudates, tubers, bark and flowers, from at least 205 plant families. However, only 10 plant families were consumed by more than half of the primate species. The largest component of primate diets was fruit (about 38% of diet), followed by leaves (about 26%), with fruit being the most consumed plant part within seven of the 10 widely consumed families: Moraceae, Sapotaceae, Rubiaceae, Anacardiaceae, Apocynaceae, Myrtaceae and Sapindaceae. By contrast, the only widely consumed family to be mainly consumed as leaves was Fabaceae. All consumed plant families showed large variation in their importance to specific primates. Most families were consumed at low levels across the primate order and only highly consumed by a few primate species. For example, the third most widely consumed family, Sapotaceae, was eaten by 66% of species, but only made up greater than 10% of diet for 17 species. For one specific primate, the highly folivorous Preuss' red colobus (Piliocolobus preussi), Sapotaceae plant parts composed almost 30% of the diet (electronic supplementary material, table S4). This extreme variation in consumption patterns for these plant families across primate species is probably due to a combination of the relative availability of different plant families in different geographic areas at both local and regional scales, variation in phytochemistry and morphology of various parts of particular plants in those areas, and the dietary adaptations of the primate species. One limitation of our dataset is that interannual and intergroup variation may not be fully captured for some species. Most species had only one study that was conducted over a period of greater than six months, highlighting the need for future dietary studies to track multiple primate populations over multiple years. Nevertheless, the consumption of Moraceae and Fabaceae emerged as highly prevalent components in primate diets (in terms of the number of primate species that consumed them, as well as their proportion across primate diets). Species of both families were consumed by multiple primate species from all major geographic areas, clades and dietary niches, unlike other plant families. Moraceae is among one of the most diverse angiosperm families (ca 1100 species [38]) and is found across a wide variety of habitats. Its fruits and leaves are consumed by a wide variety of primate species, including spider monkeys, chimpanzees, gibbons, howler monkeys, colobus monkeys and gorillas [39–44]. Ficus, the largest genus in Moraceae, is also consistently among the most diverse genera in lowland tropical forests across all three tropical biogeographic realms [45] and may thus explain the large contribution of Moraceae across primate diets (electronic supplementary material, figure S8). Furthermore, young leaves and ripe fruits of Ficus are available throughout the year since individuals tend to be asynchronous in phenology, thus offering constant food supply for primates to meet their basic nutrient needs, including high levels of calcium [46–50], and making them a keystone resource for many frugivores [43,48,51–53]. With more than 700 genera and 20 000 species, Fabaceae is the third largest angiosperm family [54]. It is the most commonly found plant family in tropical rainforests of Latin America and Africa [55,56]. A striking adaptation found among some members of this family is their symbiosis with nitrogen-fixing microbes, which results in high levels of protein [54]. Since high protein food items are often prioritized by folivorous primates, legumes may be particularly attractive food items [11,57–60]. Overall, the reliability and quality of Fabaceae and Moraceae species as food sources likely explains their prevalence and importance and highlights their role as keystone families in primate diets globally. The degree of consumption of Fabaceae and Moraceae was random with respect to the primate phylogeny, although we detected a weak association between the occurrence of consumption of Fabaceae and the primate phylogeny. The general lack of a phylogenetic signal in family-specific consumption patterns contrasts with the fact that the degree of frugivory or folivory in primates is associated with a broad suite of traits, many of which are strongly constrained by evolutionary history [24]. For example, more frugivorous primates tend to have larger brains compared to more folivorous species [61,62]. Large brain size may be particularly advantageous when exploiting asynchronously fruiting plants such as Ficus [47,48,63,64], which must be tracked at the individual tree level rather than across all individuals of the same species, thus requiring increased demands on spatial and temporal memory [62,65]. By contrast, different primate lineages may possess clade-specific dietary adaptations to deal with the potentially high fibre or toxin content of a leaf-based diet, such as the complex foregut of colobines, enlarged colon for hindgut fermentation in howlers, and specialized teeth of pithecids [66–68]. However, family-level patterns of consumption may not necessarily mirror patterns of folivory/frugivory among primates. Plant families are rarely consumed for a single plant part (e.g. fruit versus leaves), and the degree of consumption of different plant families may be more strongly driven by local or regional differences in the abundance of those plant families than variation in phytochemistry or morphology across plant families. Our SEM results support the hypothesis that space use by animals increases with larger body mass [69] and that more folivorous primates tend to occupy smaller home ranges than more frugivorous species [13]. This may be explained by the greater absolute nutritional requirements of larger bodied primates and the less patchy distribution of leaves compared to fruits [1,70]. Furthermore, folivores with larger body mass can more easily subsist on a low-quality diet due to lower energy demand per gram of body and increased gut capacity for fermenting fibre [13]. We did not, however, find any effect of the proportion of fabaceous leaves in the leaf diet on either body mass or home range size. This was surprising since fabaceous leaves are generally regarded as high-quality foliage relative to leaves from other families, as well as being widely available in tropical ecosystems [54,55]. Thus, primates that consume fabaceous leaves to a large degree may be expected to have smaller home range requirements or may be able to more efficiently digest them without having to be larger in size. This discrepancy could reflect (i) high variability in the nutritional quality of leaves among Fabaceae species, (ii) that secondary compounds may be counteracting the nutritional benefits or (iii) that high-quality leaves of legumes may be more patchily distributed than expected due to greater competition thus negating any effect of individual plant availability on home range size [71,72]. When analysing patterns of frugivory, we found that primates with a larger percentage of Moraceae in their fruit diet had smaller home range sizes. This may be explained by the greater availability of Moraceae fruits year-round compared to other fruiting tree species [47,48], suggesting this family to be uniquely important as a keystone resource to many primates [43,50,53]. When controlling for phylogeny, however, the positive relationship between folivory and body mass and negative relationship between Moraceae fruit consumption and home range size were no longer statistically significant. This suggests that these trends are primarily driven by certain primate clades (electronic supplementary material, figures S6 and S7). The negative effect of Moraceae consumption on home range size appears to be primarily driven by the cercopithecids. We acknowledge, however, that any global comparative approach of primate diets is limited by the available information on feeding records and does not usually reflect the nutritional content of the food items consumed. In addition, the quality of the same plant parts will vary within their respective families. A deeper understanding of how morphological and phytochemical variation, as well as the relative abundance, of different plant resources across plant species, genera and families shape the traits of specific primate species will add additional insights into plant–primate evolutionary relationships. Plants from the Moraceae and Fabaceae families form a core component of primate diets worldwide and may potentially be keystone resources for primates. In particular, the consumption of Moraceae may be associated with certain primate traits, such as home range size, which we argue may be due to the year-long reliability of moraceous fruits. We thus show that primates may be intricately and subtly shaped by the plants that they consume. The high proportion of fruit in the diet of primates also highlights their importance as crucial seed dispersers for many angiosperm species across ecosystems [2]. Our study also shows that comparative approaches may provide further insights into the role of diet in primate ecology and evolution. We encourage primatologists to continue with field studies on primate diets that report valuable and detailed information on plant dietary items clearly organized by plant species and part to allow for future comparisons across studies (see [4]). Full dataset and metadata, phylogeny, and R scripts used to perform data cleaning, analyses and figure generation for this study are available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.sbcc2fr40 [73]. The full list of references compiled for our dataset is also available in the electronic supplementary material [74]. J.Y.L.: data curation, formal analysis, funding acquisition, investigation, methodology, resources, software, validation, visualization, writing-original draft, writing-review and editing; M.D.W.: conceptualization, funding acquisition, investigation, project administration, resources, supervision, writing-original draft, writing-review and editing; J.V.: data curation, investigation and writing-original draft; M.D.: data curation, investigation and writing-original draft; W.D.K.: funding acquisition, investigation, methodology, project administration, resources, supervision, writing-original draft, writing-review and editing. All authors gave final approval for publication and agreed to be held accountable for the work performed therein. We declare we have no competing interests. J.Y.L. thanks the NTU Presidential Postdoctoral Fellowship. M.D.W. thanks the Tomlinson Postdoctoral Fellowship at McGill University where the ideas for this paper first materialized, the AnthroTree Workshop supported by the National Science Foundation (NSF, BCS-0923791) and the National Evolutionary Synthesis Center (NSF, EF-0905606) for training in phylogenetic methods, and Indiana University for funds that supported the development of this project to completion. W.D.K. acknowledges support from the Netherlands Organization for Scientific Research (grant no. 824.15.007) and from the University of Amsterdam via the Faculty Research Cluster ‘Global Ecology’. We would like to thank the countless dedicated primatologists whose work have made this study possible. We would also like to thank Katharine Milton, Colin Chapman, Lauren Chapman, Charles Nunn, Karline Janmaat, Jessica Rothman, Ria Ghai and Richard Mutegeki for discussions and insight that have helped us develop these ideas over the years. FootnotesElectronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5448694. References
Page 2The offshore banks and seamounts around New Caledonia are recognized for their rich diverse marine fauna, occurring both at sublittoral (0–200 m) and upper bathyal depths (200–1000 m) [1]. This fauna includes several phylogenetic relicts, which appear to be survivors of ancient faunal groups that have otherwise gone extinct. Examples include: one of the few surviving species of nautiloid cephalopods (now the marine emblem of New Caledonia); three of the eight extant species of cyrtocrinids, an order of crinoids that was highly diverse during the Mesozoic [2]; and one of the two modern species of glypheoid lobsters that otherwise went extinct in the Eocene [3]. This fauna has been discovered though the exceptional tropical deep-sea benthos (TDSB) programme led by the Muséum national d'Histoire naturelle (MNHN) and Institut de Recherche pour le Développement (IRD) who have conducted over 40 scientific surveys of the deep-sea environment around New Caledonia (http://expeditions.mnhn.fr/program/tropicaldeep-seabenthos). Here, we describe a new species, genus and family of Ophiuroidea (brittle-stars) collected by the EXBODI expedition from the Durand Bank on the Loyalty Ridge, located in the south-west of the New Caledonian Exclusive Economic Zone. A recent phylogeny of living Ophiuroidea based on 265+ kbp of DNA sequence [4,5] allows us to place this new lineage within the tree of life with precision. The discovery that the external micro-structure of lateral arm plates (LAPs) is phylogenetically conserved on ophiuroids [6–9] facilitates a comparison with fossil taxa. The holotype (and only known specimen) of the new species was collected on 13 September 2011 by 4 m beam trawl from the RV Alis on the EXBODI expedition (stn CP3849) on Banc Durand, New Caledonia (22°2.7′ S, 168°40.8′ E, 360–560 m). This bank is on the southern Loyalty Ridge and is likely of volcanic origin [10]. The bank is now partially covered in sediment on the flanks and coral reef on the flat summit, which lies just below the sea surface at approximately 10 m depth [11]. The specimen was fixed and preserved in ethanol and is stored at the MNHN, Paris, registration no. MNHN IE.2007.6821. The microfossils were hand-picked from washing residues of sediments from the Lower Jurassic (Toarcian, Serpentinum and Bifrons Ammonite zones, 180 Myr) of Feuguerolles, Normandy, France (see [12] for site details). The fossil specimens are stored at the Natural History Museum of Le Mans (France) (collection acronym MHNLM). We sequenced genomic DNA from the new specimen using the exon-capture procedure, as described in Hugall et al. [13] and http://dx.doi.org/10.5061/dryad.db339, with a minimum coverage limit of five reads. Following the approach taken in our previous phylogenomic analyses [4,5,14], we included this new data in a selection of 195 species representing all major lineages and family crown clades (electronic supplementary material, table S1). After standard data filtering [13] the final dataset comprised 265 158 sites in 1496 exons from 416 genes, with the 195 taxa and is 94% (median) complete. For computationally intensive Bayesian dating analyses we used a subset of 79 taxa (phylogenetically spread but focused on the Ophiacanthida) with a subset of 219 exons, one per gene, amounting to 45 981 sites (see electronic supplementary material, table S2 and references [5] and [15] for further details). We generated phylogenetic trees using maximum-likelihood (RAxML v. 8.1.20) [16], Bayesian (BEAST v. 2.4.7) [17] and species tree approaches (ASTRAL II v. 5.5.10) [18], the latter used to assess artefacts of sequence concatenation [18]. All concatenated data RAxML and BEAST analyses used a three-part codon position GTR + G + F site rate heterogeneity model with base frequency estimation, as selected via ModelFinder in IQ-TREE [19]. Phylogenetic trees for both 195 taxa and 79 taxa are generated in RAxML by running 200 non-parametric bootstraps (the –f i command) that retain branch lengths per bootstrap. These analyses were run on the CIPRES supercomputing facility. As with our previous phylogenomic analyses, RAxML trees were post hoc rooted by splitting the branch between the two superorders Euryophiurida and Ophintegrida in proportion to the split estimated by the 425-gene echinoderm transcriptome phylogenomic study of O'Hara et al. [6]. These trees were then used in dating analyses described below. Following O'Hara et al. [4] we also inferred a species tree using ASTRAL II, with local posterior support values [20], from 353 (the most data rich, electronic supplementary material, table S2) separate gene trees, drawn from the full 265 kb 195 taxa dataset. We used genes (rather than exons) as independent loci [21]. Unrooted gene trees were generated by RAxML (−f d command) using a single partition GTR + G + F model (due to reduced information per gene as compared to the full concatenated data). These gene datasets comprised (median) 189 taxa and 513 sites, with each of the 195 taxa having (median) 343 gene trees (330 for specimen IE.2007.6821). The ASTRAL species tree figure is presented as rooted for ease of comparison. To investigate the accuracy and precision of dating the new lineage, we used both penalized likelihood rate smoothing (PLRS) [22] and Bayesian relaxed-clock methods (BEAST). Both the PLRS and BEAST analyses used slightly modified versions of previously published [4,5] calibrations (electronic supplementary material, table S3): 11 fossil-based node minimum constraints spread across most major lineages (invoked as either uniform, γ or exponential priors), plus a root secondary calibration normal prior (mean = 260.0, σ = 7.0 Ma) based on the results of O'Hara et al. [6]. PLRS is a widely used maximum-likelihood method for producing time-calibrated chronograms from phylogenetic trees but is not well-suited to providing uncertainty limits compared to Bayesian methods [23,24]. In order to provide some measure of uncertainty, we included three sources of variation: (1) 200 non-parametric bootstrap trees to account for branch length and topology uncertainty, (2) calibrated (in addition to minimum constraints) with a root age drawn from a normal prior distribution (above, as also used in BEAST analyses) and (3) smoothing factor drawn from a γ prior distribution (α = 1.5, β = 4.0; median = 5.8), based on smoothing factor cross-validation distribution fit (reflecting uncertainty in the relaxed-clock model, [22]). PLRS was run in r8s v. 1.7 using the ADD penalty function and TN optimization. This approach was taken for both the 195 and 79 taxa trees, and trees and age distributions compared, focusing on the stem age distribution of the target sample (IE.2007.6821). These sets of 200 PLRS chronograms were then summarized as a maximum clade credibility (MCC) tree using FigTree v. 1.3.1 (http://tree.bio.ed.ac.uk/software/figtree). Bayesian dating estimates were obtained using BEAST2 v.2.4.7 using sequence evolution data, relaxed-clock and calibration prior models similar to those used before [4,5,14]: a codon position GTR + G + F sequence evolution model, with lognormal relaxed-clock model and Yule tree prior (both using 1/X parameterization). Analyses used a PLRS tree as a starting tree and were a priori rooted by enforcing the reciprocal monophyly of the two superorders. Two independent runs of 100 M generations (1/10 000 sampling) were executed with 20% burnin, and assessed for convergence; each run returning posterior likelihood ESS greater than 100 [25]. The two runs used the alternate standard ‘uniform' or ‘exponential' versions of the fossil-based calibration constraints (electronic supplementary material, table S3). MCC trees and age distributions were evaluated separately and combined, focusing on the stem age distribution of the new species, and compared to the PLRS results. The full PLRS and BEAST MCC trees were presented in electronic supplementary material, figures S1 and S2. Light images of the holotype were taken at Museums Victoria with a Visionary Digital Integrated System, using a Canon 5D Mark II camera with EF100 mm and MP-E65 mm macro-lenses, and montaged using Zerene Stacker v. 1.04 software. Ossicles of the arm skeleton were extracted from a proximal arm portion of the holotype, macerated in household bleach and rinsed in tap water, then mounted on aluminium stubs and gold-coated [8]. Scanning electron microscope (SEM) images were taken with a JEOL Neoscope JCM-5000 at the Natural History Museum Luxembourg. Micro-CT scanning was performed with a Phoenix Nanotom m operated using xs control and Phoenix datos|x acquisition software (Waygate Technologies). The specimen was mounted by wrapping in ethanol soak gauze, secured by bubble wrap, and placed within a plastic specimen jar. An X-ray energy of 50 kV and 300 mA was used for two separate scans, one at a coarser voxel resolution of 25 µm to capture the full body of the specimen, and a higher resolution scan at 10 µm focusing on a region of interest encompassing the main body. Scans were run for 10 min in a fast scan mode collecting 1200 projections through a 360° rotation of the specimen. Volume reconstruction of the micro-CT data was performed using Phoenix datos|x reconstruction software (Waygate Technologies) applying an inline median filter and ROI filter during reconstruction. The data were exported as 16-bit volume files for imaging and analysis in Avizo (Thermo Fisher Scientific). This published work and the nomenclatural acts it contains have been registered in ZooBank, the proposed online registration system for the International Code of Zoological Nomenclature (ICZN). The ZooBank LSID (Life Science Identifiers) for this publication is (urn:lsid:zoobank.org:pub:7DC3B82F-C07A-483B-A1C5-8FDEEB718938). LSIDs for the new taxa are given below. From the only known specimen we obtained 89% of our exon-capture target of 265 158 bp. RAxML bootstrap support for the position of the new species as sister to the rest of the suborder Ophiacanthina was 70%/79% in the 195/79 taxa analyses (figure 1; electronic supplementary material, figure S1a,S2a). Both the ASTRAL species tree (electronic supplementary material, figure S1b) and BEAST relaxed-clock tree (electronic supplementary material, figure S2b) returned this position, with posterior supports of 1.00. The remainder (minority) of the RAxML bootstraps placed the new species as sister to Ophiodermatina (rather than Ophiacanthina), hence complete support for it lying within the Ophiacanthida order crown as defined by O'Hara et al. [4,14]. As support for nodes above and below this rank is in all cases unanimous, we can exclude the new species from lying within or on the stem of any single described family (in particular the Ophiobyrsidae). Figure 1. Phylogeny of the Ophiuroidea showing the position and estimated node age range (bars on nodes) of Ophiojura for the RAxML/PLRS 195-taxon analysis. Width of the collapsed clades not scaled to richness. All nodes are fully supported on the 195- and 79-taxon RAxML (bootstrap), BEAST (posterior) and ASTRAL II (local posterior), except the Ophiojura node which has values of 70, 79, 1.00 and 1.00, respectively. (Online version in colour.)
Our previous fast bootstrap GTRCAT model analysis [5] returned approximately 90% support for the position of the new species, broadly in proportion to the non-parametric bootstrap as is theoretically expected of such methodology (70–95% [16,26]). The difference between thorough bootstrap analyses likely reflects taxon sampling [27]. The significant species tree local posterior support (1.00) compared to the slightly equivocal bootstrap may be expected of a short internode due to a short time interval relative to the mutation/fixation process [21]. Hence the new species can confidently be placed as the sister to all currently known Ophiacanthina and therefore must belong to a distinct family. Results of the dating are summarized in figure 2. The dating of the new species node was very similar between different methods, taxon sets and priors, viz PLRS (n = 195 taxa): median = 180 (95% CI = 163–215), PLRS(79): 178 (154–205), BEAST-Uniform(79): 187 (149–211), BEAST-Exponential(79): 181 (143–215) Myr. This consistency suggests that methodological issues are not critical but the results have limited precision. While we have considerable data to estimate topology and branch length, our calibrations are broad, and perhaps more importantly, there is considerable apparent rate variation among major lineages—(see non-clock trees in [4–6,13]). Such variation inevitably creates complex solutions for any relaxed-clock model, even when there are precise calibrations [23,24]. Some of this can be seen in the age distributions (figure 2) with a bimodal pattern in PLRS and broad distributions and relatively poor ESS per MCMC chain length in the BEAST results. Nevertheless, in all cases the 95% CI fall mainly within the Jurassic Period and it is the oldest single species (longest branch, electronic supplementary material, figure S3) in our global sample of nearly 1000 ophiuroid species [5]. Figure 2. Ophiojura lineage age dating estimate distributions. PLRS transformed RAxML trees from the full 265 kb exon data with either 195 or 79 taxa, BEAST analyses used 79 taxa with the reduced 46 kb exon data, using either the standard ‘uniform' prior or ‘exponential’ prior age calibration sets. Drawn from 200 PLRS trees and 1000 BEAST posterior tree samples in 10 Myr bins. (Online version in colour.)
Phylum: Echinodermata Class: Ophiuroidea Order: Ophiacanthida Suborder: Ophiacanthina Family: Ophiojuridae fam. nov. urn:lsid:zoobank.org:act: 5854ACA1-6994-4379-88DF-FA2F3F53E891 Genus: Ophiojura gen.nov. urn:lsid:zoobank.org:act: EE0F43F8-0773-4641-A47A-D9BDF0AD4386 Species: Ophiojura exbodi sp. nov. urn:lsid:zoobank.org:act: D83078F6-A23D-41EF-B2BA-3094159B2F5A. Registration no.: MNHN IE.2007.6821. EXBODI stn CP3849, Banc Durand, New Caledonia, 22°2.7′ S, 168°40.8′ E, 360–560 m, 13 September 2011. The genus and family are named after the Jura Mountains (noun), type locality for the Jurassic period, with the prefix ‘Ophio' derived from the Ancient Greek word for serpent. The species is named after the EXBODI expedition (noun) that collected the holotype, organized by Dr Sarah Samadi of MNHN, Paris. LAPs higher than long, with a long proximal extension that curves around the large open tentacle pore (figure 3f,h). Arm spine articulations a distinctive ‘pig-snout'-shape (figure 4a,c), with a large muscle opening bordered ventro-proximally by a small, thin ventral lobe, a much larger, swollen lobe dorsally bordering both the muscle opening and the slightly smaller and widely separated nerve opening. Distinctive hook-shaped arm spines on distal arm segments, with nerve and muscle attachments, two rows of thorns and a large, recurved terminal claw. Contiguous ventral arm plates (VAPs), longer than wide, constricted laterally around the tentacle pores (figure 3f). Dorsal arm epidermis contains perforated thin ossicles over the vertebrae but no solid plates (figure 3e). Vertebrae with zygospondylous (yoke-shaped) articulation (figure 4h,i). Oral and tooth papillae with a spine-like core surrounded in soft tissue (figure 3d). Three peristomial plates that sit dorsally over the water vascular ring/neural groove (figure 3g). Figure 3. Ophiojura exbodi (MNHN IE.2007.6821) holotype. (a) Ventral view of disc, (b) ventral view of animal including remaining arm fragments, (c) dorsal view of arms, (d) teeth papillae, (e,f) partially macerated arm segments from mid arm, (e) in dorsal view, with inset showing small perforated epidermal plates, (f) ventral view. Rendered CT scan of the (g) dorsal and (h) ventral surfaces. AdSh, adoral shield; DOP, distal oral plate element; LAP, lateral arm plate; LOPa, lateral oral papillae, OSh, oral shield; PEP, peristomial plate; POP, proximal oral plate element; TP, tentacle pore; VAP, ventral arm plate. (Online version in colour.)
Figure 4. Ophiojura exbodi (MNHN IE.2007.6821) holotype. (a–k) SEMS of internal ossicles from the proximal to middle of the arm. (a–d) Lateral arm plates with dorsal side upwards and distal side with arm spine articulations, (a) external, (b) internal views of proximal plate, (c) detail of lower spine articulations, (d) external view from the middle of the arm. (e,f) Arm spines, (e) dorsal-most and (f) ventral spines. (g) Ventral arm plate in external view, with distal side upwards. (h–k) Vertebrae, (h) distal with dorsal side up, (i) proximal with dorsal side up, (j) dorsal with proximal side up and (k) ventral views (with proximal side up). (l,m) Distal arm segments with dorsal side up and distal side to the right, showing hook-shaped arm spines.
Disc diameter estimated (from remnant pieces) to be at least 30 mm, longest arm piece greater than 80 mm (from mouth centre, 45 segments). Eight arms (figure 3a,b), 3.4 mm across at base. Colour pale tan/white, no indication of any pigmented markings. Dorsal disc mostly absent, small remaining pieces attached to arms are composed of thick dermis typically without spines or plates, although one piece has a tiny oval plate that may be the radial shield. No sign of any remaining genital plate. Jaws elongate, 3.7 mm long, 3 mm wide at base. Adoral shields long and narrow, extending from half the length of the jaw past the oral shield, meeting proximally, widely separated distally, with a truncate to concave distal margin (figure 3h). Ventral jaw surface proximal to adorals and between lateral oral plates is covered in skin. Oral shields 2× wider than long, convex proximally, truncate distally with rounded lateral angles, depressed in centre, four are larger (presumably madreporites), almost as long as wide, rounded, orientated obliquely to jaw (figure 3h). Cluster of four to six tooth papillae at jaw apex, each with an inner-spine-like core that is covered in a conical skin sheath, the inner-spine terminating in one or more thorns that emerge from the skin layer (figure 3d). A row of five teeth in vertical series widened at base but tapering to a point, followed by a cluster of three spine-like dorsal teeth (electronic supplementary material, figure S5a,b). Exposed proximal end of the oral plates long and rounded on the ventral surface, forming elongated jaws and narrowly articulating radially before their proximal extent leaving a cavity between the oral and dental plates. The distal section of the oral plates have large dorsal extensions for inter-vertebral muscle attachments (higher than succeeding vertebrae) and a single dorsal foramen for the buccal water vascular canal on the abradial side that bifurcates within the oral plate to serve both the first (dorsal) and second (ventral) oral tentacles (electronic supplementary material, figures S4 and S5). The seven to eight oral papillae occur deep in the jaw slit along the lateral ridge of the oral plate, interleaving with series on opposing jaw. These papillae also have a spine-like core with a triangular skin sheath that is webbed to their neighbour at their bases. Distal third of jaw slit without oral papillae. One remaining spiniform adoral shield spine still present adjacent to 2nd oral tentacle pore, well outside of the jaw silt. The water vascular ring/neural groove protected dorsally by five overlapping plates (electronic supplementary material, figures S4 and S5), two falcate plates [28] on the proximal side of the groove that also protect the dorsal side of the first tentacle, two supplementary peristomial plates that sit over the groove, and a central droplet-shaped peristomial plate that sits over the external interradial muscle [29] and possibly the stone canal. No dorsal arm plates (DAPs, figure 3c), although some microscopic perforated ossicles appear to be embedded in the epidermis above and near the vertebrae (figure 3e). LAPs roughly ovoid, higher than long (0.8 × 0.68 mm) with a concave ventral border and a protruding ventro-proximal tip (figure 4a,d) that is contiguous with the VAP. A small rise dorso-proximally with denser stereom and a poorly defined, weakly prominent spur on the proximal edge of the LAP (figure 4a,d). Six arm spine articulations on flattened distal surface of LAP (figure 4a,d), each with a large muscle opening bordered ventro-proximally by a small, thin ventral lobe, a much larger, swollen dorsal lobe dorsally bordering both the muscle opening and the slightly smaller nerve opening, widely separating both openings to form a distinctive ‘pig-snout'-shaped spine articulation (figure 4c). Distance between spine articulations strongly increasing dorsal-wards, with dorsal-most spine articulation widely separated from the other five. Interior LAP surface with two knob-shaped vertebral articular structures of dense stereom, one central and one near the distal margin and a ventro-proximal perforation (figure 4b). VAPs composed of open stereom, 2.0 mm long, 1.4 mm wide at widest point, 0.75 mm wide at narrowest point (figures 3f and 4g). Distally truncate, with distal-ward-divergent lateral sides, constricted mid-length around the large open tentacle pores, proximally flared, with a concave proximal border. Six arm spines. Upper arm spine (2.1 mm) with base slightly enlarged of open stereom, shaft with a shallow wide furrow along the dorsal side with a proximal pore and series of stereom pores along its length, some small denticulations laterally near the tip (figure 4e). Five ventral spines, subequal (1.2 mm), with an expanded base, and central series of stereom holes along a narrow groove, denticulations prominent laterally, long along dorsal side, including from the truncate apex (figure 4f). On distal arm segments, four to five dorsal arm spines replaced by distinctive, hyaline hooks (figure 4l,m), composed of a coarsely porous base with a nerve opening and a lateral muscle attachment extension, followed by two rows of thorns apically converging into a large, recurved terminal thorn. Some hooks of the same LAP seemingly opposable. No tentacle scales observed. Vertebrae with zygospondylous articulation, distal face with two lateral articulation surfaces and a median condyle with a concave face (figure 4h), proximal surface with median transverse bar and two ventral protuberances (figure 4i). Large dorsal and ventral muscle fossae, dorsal one with three to four areas of denser stereom. Deep dorsal furrow down the midline with several large apertures in the centre (figure 4j). Dorsal lateral wings relatively flat without pronounced curved grooves, nodulated external surface. Large podial basins (figure 4k). The unusual ‘pig-snout' form of the arm spine articulations is similar to Ophiobyrsa rudis Lyman, 1878, the type species of the type genus of Ophiobyrsidae [7]. The new species also shares other anatomical features with species within the Ophiobyrsidae as defined in O'Hara et al. [14], including the shape of the jaw, a series of pointed oral papillae along the oral plate ridge, a cluster of similar-shaped tooth papillae at the jaw apex, an adoral shield spine, reduced or absent DAPs, elongated contiguous VAPs, large tentacle pores that lack scales and reduced denticulate arm spines [7,30]. However, the vertebral articulation in the Ophiobyrsidae has been described as being streptospondylous (with hourglass-shaped articulation surfaces and no median condyle) [7,30]. The organization of the teeth can differ within that family, and there can be one to four vertical rows of teeth. Another important difference pertains to the ornamentation of the LAP external surface, the new species has areas of denser stereom and a spur on the proximal edge which are absent in the Ophiobyrsidae. The muscle and nerve openings of the spine articulations of the Ophiobyrsidae also tend to be separated by a thin, low ridge, whereas in the spine articulations of the new species, the muscle and nerve openings are separated by a much wider, swollen ventral-ward expansion of the dorsal lobe. The Ophiocamacidae can have a cluster of tooth papillae, but in general Ophiocamax species have a more obvious ophiacanthid morphology, especially in having a single column of teeth dorsally, oral papillae on the ventrolateral edge of the oral plate, small tentacle pores, multiple tentacle scales, robust arm spines, robust DAPs, an armoured disc and volute-shaped arm spine articulations. The taxonomic distinction of large- and small-pored Ophiacanthids (sensu [31]) is not congruent with the latest molecular trees [4,14], as large pores can also occur in the extant Ophiotomidae (e.g. Ophiotoma and Ophiopristis species). The other families within the Ophiacanthida have a disc either covered in plates or spines/granules, the tentacle pores on the arm are typically small and covered by tentacle scales (with a few exceptions, most notably Ophiotoma) and the teeth and oral papillae are generally solid, blunt or rounded. Furthermore, the Clarkcomidae and Ophiopteridae have their oral frame modified into a ‘buccal funnel' [28] to crush food before ingestion, with ventral teeth modified into multiple rows of blunt papillae. The new species also bears some resemblance to species in the Ophioscolecida (historically placed in a polyphyletic Ophiomyxidae), including the large tentacle pores, presence of tooth papillae, lack of flabelliform teeth, oral papillae that form along the oral plate ridge (rather than on a ventrolateral margin), the large 2nd oral pore with a podia that cannot enter the mouth slit when closed and the presence of an adoral shield spine [7,14,28]. The genus Ophioleptoplax in particular has a similar oral frame, VAPs, tentacle pores and reduced DAPs. However, it differs in having long smooth solid arm spines. It supposedly has no radial shields, but this was doubted by Martynov who discovered minute radial shields in similar species [7]. Multiple unbranched arms are unusual in the Ophiuroidea but their presence is not significant for higher-level classification. Of the almost 2100 accepted species of Ophiuroidea, only four species Astrochlamys sol Mortensen, 1936, A. timoharai Okaniski & Mah, 2020 (Order Euryalida), Soliophis bakeri Okanishi & Mah, 2020 (Ophioscolecida) and Ophiacantha decaactis Belyaev & Litvinova, 1976 (Ophiacanthida) have been recorded with 10 or more simple arms; one Ophiacantha enneactis H.L. Clark, 1911 with nine; and two species Ophiacantha opulenta Koehler, 1908; Ophiacantha vivipara kerguelensis Studer, 1876 can rarely have up to eight arms. Some fissiparous species, which usually have six arms, have also been recorded occasionally with seven, including Ophiocomella sexradia (Duncan, 1887) (Ophiacanthida) and various Ophiactis species (Amphilepidida). The presence of 10 or more simple arms has been associated with an arborescent habitat [32]. Wilkie & Brogger [29] have divided ophiuroid oral frames into two broad categories. Ophiojura conforms to type ‘A' with spine-like teeth, oral plates with distally projecting oral and aboral convexities, and large peristomial plates that cover the circum-oral water vascular ring. Since this is the likely ancestral state [29], it is not a good indicator of phylogenetic relationships. Of the few figured species, the triple peristomial plate arrangement is reminiscent of that in Ophioderma rubicunda ([28], fig. 12h) in the suborder Ophiodermatina, although in that case the jaw is much shorter and the central peristomial plate small and triangular. The LAPs of the new species are very similar to previously unpublished microfossils from Lower Jurassic (180 Myr) of Feuguerolles, France (figure 5). The distinctive morphology of these LAP microfossils, including the ‘pig-snout' form of the spine articulations with a thick ventral-ward extension of the dorsal lobe separating the muscle and nerve openings, and the presence of differentiated stereom and spurs on the external LAP surface, are consistent with the new family Ophiojuridae. These microfossil LAPs appear to be more heavily calcified than their extant counterparts, with thicker plates and more developed vertebral articular surface and tentacle notch. Figure 5. Fossil Ophiojuridae from the Lower Jurassic (Toarcian, Serpentinum and Bifrons ammonite zones, 180 Myr) of Feuguerolles, Normandy, France. SEMs of lateral arm plate microfossils with dorsal side up. (a–c) MHNLM 2020.1.4, proximal to median lateral arm plate in (a) external and (b) internal view and (c) with detail of spine articulations. (d) MHNLM 2020.1.5, proximal lateral arm plate in external view.
Previously published LAP microfossils with ‘pig-snout'-shaped spine articulations, all assigned to the extinct genus Ophiojagtus [31], differ in lacking external surface ornamentation and in having the muscle and nerve openings separated by a thin, low ridge, suggesting assignment to the Ophiobyrsidae. Furthermore, at least in the case of Ophiojagtus argoviensis (Hess, 1966), the associated vertebrae have a streptospondylous articulation corroborating the ophiobyrsid affinities. We have established a new family of brittle-stars based on a single damaged specimen. While this is not typical taxonomic practice, we justify publication on the basis of robust phylogenomic and detailed morphological data, and to facilitate the collection and identification of new specimens. The situation is similar to the first description of an extant glypheoid lobster, based on one old dry specimen found in the Smithsonian Museum [33], which subsequently inspired a French expedition to the Philippines to locate more material [34]. Our phylogenetic analyses estimate that the Ophiojura lineage arose in the Jurassic or late Triassic periods, most probably from 160 to 200 Myr. It is the longest branch for any single species on our extensive molecular phylogeny of the Ophiuroidea [5]. Despite the substantial survey effort that has occurred in New Caledonian waters, we are unaware of any other specimen that could plausibly be confamilial with this specimen. The fossil record indicates that the Jurassic was a period of elevated diversification of early ophiacanthid lineages, with many species occurring in the relatively shallow water of the continental shelves (0–200 m) [31]. The Ophiojuridae fossils presented herein confirm that the family was present by the Early Jurassic (180 Myr). However, several other derived ophiacanthid clades have a confirmed Early Jurassic fossil record [31], implying that the Ophiojuridae must have diverged before this period. Fossil and molecular evidence [35] has shown that the majority of the extant ophiacanthid clades originated in deeper waters (below 200 m), implying that their occurrence at shallower depths in the Jurassic represented a temporary expansion of their bathymetric range [31]. The fossil record is yet too sparse to confirm a similar pattern for the Ophiojuridae. It is noteworthy, however, that the only known fossils of the Ophiojuridae were found in strata that yielded one of the oldest diversified shallow-water cyrtocrinid faunas [36]. Cyrtocinids, in fact, are another clade with an assumed deep-water origin followed by a temporary expansion to the shallow epicontinental seas of the Jurassic [36]. We, therefore, speculate that additional fossil members of the Ophiojuridae which bridge the extensive temporal gap between the Jurassic and the present, may be found either in the extremely poorly sampled deep-sea fossil record, or in ancient shelf settings affected by increased bathymetric faunal exchange [31]. Although the Ophiojura lineage is ancient, and the taxa is not known to have any close extant relatives, this does not imply that Ophiojura is ‘primitive' [37] or that all its characters were present in the middle Jurassic. In particular, the presence of eight arms is not an indicator of deep taxonomic relationships. On the other hand, arm spine articulation morphologies do reflect the higher classification of the Ophiuroidea [6] and the shape outlined here may be an indicator of which of the many Jurassic fossils actually belong to the stem of the extant Ophiacanthina. The discovery of a new lineage on a seamount could be seen as evidence for the hypothesis that seamount faunas can be highly endemic [1], reflecting the age and geographical isolation of the seamounts. However, initial reports that endemism is a feature of seamounts have not been supported by more detailed studies [38]. Many seamount species are widespread and apparently good dispersers [38]. Nor do seamounts have consistently higher species richness than areas of continental slope [39]. However, our new discovery does reflect the deep and rich evolutionary diversity that is present in the upper tropical bathyal of the Indo-Pacific [40]. This biome, particularly between 200 and 1000 m, is a reservoir of palaeoendemism (formerly widespread but now restricted lineages) both for ophiuroids [40] and other taxa (see introduction). Banc Durand, the type location of Ophiojura, also has been noted as a hotspot of octocoral phylodiversity [41]. To protect the rich tropical deep-sea fauna, the entire New Caledonia Exclusive Economic Zone was placed in 2014 within a vast marine protected area, the Coral Sea Nature Park. Input and output files from the phylogenetic analyses and images from the micro-CT scan of the holotype, the latter including stacked vertical and horizontal images, rendered images of the dorsal and ventral surfaces, and a rotating animation, are available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.18931zcx3 [15]. Lists of materials, fossil calibration points and expanded phylogenetic trees are provided as electronic supplementary material. Requests to examine specimens can be directed to the relevant museum collection (MNHN and MHNLM). The data are provided in the electronic supplementary material [42]. T.D.O.: conceptualization, data curation, formal analysis, investigation, methodology, project administration, resources, supervision, visualization, writing-original draft, writing-review and editing; B.T.: formal analysis, investigation, methodology, resources, writing-original draft, writing-review and editing; A.F.H.: formal analysis, investigation, methodology, writing-original draft, writing-review and editing. All authors gave final approval for publication and agreed to be held accountable for the work performed therein. The authors declare no competing interests. Muséum national d'Histoire naturelle, Paris provided funding to identify specimens from expeditions to New Caledonia. The Australian National Environmental Science Program's (NESP) Marine Biodiversity Hub provided funding to sequence the taxonomic diversity of the Ophiuroidea around Australia and neighbouring nations. We thank the scientists and collection managers of the Muséum national d'Histoire naturelle, Paris who have collected, preserved and curated the extraordinary deep-sea fauna around New Caledonia; the Melbourne TrACEES (Trace Analysis for Chemical, Earth and Environmental Sciences) Platform for access to the micro-CT scanner and Dr Jay Black (University of Melbourne) for technical support; Caroline Harding and Mark Nikolic (Museums Victoria) with assistance with the macro photography, and Marc Chesnier (Cresserons, France) who collected the fossil material. FootnotesElectronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5448691. References
Page 3Globalization has contributed to a surge in the incidence, severity and spread of emerging infectious diseases (e.g. [1,2]). Emerging diseases of wildlife are particularly important to global biological diversity as they can cause devastating population declines and exacerbate other threats such as habitat loss, overharvesting, invasive species and climate change [3–6]. Recent advances in the study of disease emergence and spread integrate epidemiological and genetic data to test theoretical predictions about the ecological history of the pathogen given the underlying evolutionary signal [7,8]. However, most applications of this approach have been for quickly evolving pathogens (i.e. RNA viruses) and those that directly impact human health. There have been a handful of studies applying methodological advances in genetic epidemiology to emerging wildlife diseases (see recent reviews [9,10]), but such frameworks are still largely underutilized. Amphibians are declining worldwide [11,12]. One of the major drivers of amphibian declines is the global spread of the disease chytridiomycosis, caused by the fungal pathogen Batrachochytrium dendrobatidis (Bd) [13]. Bd infects the keratinized skin cells of susceptible host species, disrupts vital amphibian skin functions and can cause mortality [14]. In some cases, Bd infections can spread quickly across individuals, populations, and species leading to pizootic outbreaks and population and community collapses [15,16]. Since the earliest observations of Bd-related die-offs in the late 1990s, Bd has emerged as a global threat to amphibian biodiversity and now impacts amphibians on every continent where they are present [12]. Bd has a complex evolutionary history with multiple lineages found in different parts of the world. Phylogenetically, Bd is characterized by several early branching lineages endemic to different regions (BdCAPE, BdASIA1, BdBrazil/ASIA2 and BdASIA3) and one more recently derived hypervirulent panzootic lineage (BdGPL). BdGPL has been linked to declines of amphibian communities around the world and is the only Bd lineage with a truly global distribution [6,17,18]. Whole-genome data have been important for revealing the dynamics of BdGPL spread [6,19]. BdGPL typically exhibits little phylogenetic or spatial genetic structure (with the exception of two subclades BdGPL-1 and BdGPL-2) [20,21], suggesting that this lineage spread rapidly around the world [18,19]. Moreover, compared to other Bd lineages, BdGPL genomes have fewer pairwise genetic differences among them and highly variable genetic diversity values [6]. Observations of minimal pairwise genetic differences are consistent with rapid BdGPL spatial radiation, and variability in genetic diversity suggests episodes of population size fluctuation. However, we still lack a connection between our understanding of Bd evolutionary history at a global scale and regional Bd emergence and spread. Two of the most emblematic BdGPL-related declines occurred in the montane amphibian communities of the Sierra Nevada of California and Central Panama. In the Sierra Nevada of California, mountain yellow-legged frogs (Rana sierrae/muscosa), were historically one of the most abundant vertebrates [22]. Over the last century, these frogs vanished from more than 90% of their historic range, and Bd (along with invasive fish) was a significant factor in their decline [23]. Available information suggests that Bd has been spreading across the Sierra Nevada since at least the 1960s [24,25] and has caused epizootics and subsequent extirpations in hundreds of populations [16,26,27]. Some populations that experienced Bd-related declines are beginning to rebound, but remaining naive populations are still at risk for Bd epizootics [28]. Similarly, in Central America, amphibian population declines were first observed in the late 1980s [29,30]. As Bd spread southeast into Central Panama starting in the early 2000s [15], many susceptible amphibian host species declined—or even disappeared completely—across the region [15,31,32]. Although some species seem to be recovering [32], Bd-related declines have fundamentally reshaped these tropical communities [31,33,34]. From an epizoological perspective, amphibian declines in the Sierra Nevada and Central Panama appear quite similar. In both regions, initial detection of Bd was followed by devastating outbreaks and host mortality. Patterns of decline in both the Sierra Nevada and Central Panama also appear to provide evidence of a ‘wave’-like spread of Bd across the landscape [16,35]. Pathogen prevalence and population decline data in both systems suggest that new infections appear in a predictable spatial direction and that Bd outbreaks move a predictable distance each year [15,16,35]. Coupled with a global phylogenetic view of Bd, the prevailing hypothesis suggests that BdGPL is a recent invasive pathogen in these two regions [36]. However, epizoological data based on observed outbreaks and host outcomes may or may not reflect the true history of Bd arrival and spread. The Sierra Nevada and Central Panama differ dramatically in climate, habitat, and amphibian community composition. Ecological and environmental factors impact the physiological limits of Bd, transmission dynamics across the landscape and host immunity; extensive evidence suggests these factors influence prevalence and host disease outcomes [20,37–39]. Stark differences between regional environments likely contribute to different Bd dynamics in the Sierra Nevada and Panama, making their apparent similarities in disease outcomes all the more intriguing. Although it is often assumed that Bd arrived recently and spread in a wave-like fashion in both regions, it is possible that different evolutionary histories of Bd underlie these observed patterns. Molecular data can reveal nuances of a pathogen's history that cannot be obtained by field observations alone. Genetic and genomic approaches have previously been used to investigate the evolutionary history of Bd at regional and global scales [6,18,19,40]. However, most studies of the evolutionary history of Bd in emblematic systems like the Sierra Nevada and Central Panama have relied on a small number of Bd isolates for any one region. Live and pure Bd cultures have been the source of high-quality DNA for genomic sequencing (e.g. [6,19,41]) but are inherently challenging to obtain, isolate and maintain. Low sample sizes and poor spatial coverage has made it difficult to test fine-scale hypotheses about Bd emergence and spread. However, advances in sequencing technology now allow for leveraging fine-scale sampling of frog skin swabs, previously used to determine Bd presence/absence and load, to robustly characterize Bd genotypes across relevant spatial scales [42]. Thus, we can now test whether patterns of Bd emergence that appear similar across systems result from shared underlying processes. We used fine-scale genetic sampling to investigate assumptions about the history of BdGPL in the Sierra Nevada and Central Panama. Using non-invasive skin swabs collected across similar spatial and temporal scales, we targeted hundreds of loci across the Bd genome to examine the hypothesis of recent Bd emergence and unidirectional epizootic spread in these two emblematic systems. Our work provides an in-depth understanding of pathogen evolutionary dynamics in natural systems and an example for how researchers should not expect infectious diseases to emerge and spread similarly across the globe. We used skin swab DNA samples collected from the Sierra Nevada and Central Panama across similar timescales (2011–2017) and across equivalent spatial scales (approx. 130 km across Euclidean distance between furthest two sites) (figure 1a). Sites are defined as collections of lakes and streams that cluster together geographically within a region. We sampled 10 sites from both the Sierra Nevada (nsamples = 130; nspecies = 2) and Central Panama (nsamples = 80; nspecies = 17). Sierra Nevada samples comprised skin swabs from two sister species of frogs (R. sierrae/muscosa) [23] and Central Panama samples comprised skin swabs from 16 different frog species. Additionally, we included 120 previously published sequenced samples from a global Bd dataset, comprised of samples across 59 frog species from six continents (see electronic supplementary material). The global samples were all previously positioned within the BdGPL clade based on a comprehensive assessment of hundreds of Bd samples [17] and serve to add a global context to levels of genetic structure and diversity observed in the Sierra Nevada and Central Panama. Figure 1. Study system map and principal component analysis of within region genotypes. (a) Map of sites sampled in the study in the Sierra Nevada and Central Panama. (b) PCA within Sierra Nevada samples, coloured by the major site. Samples cluster by site, suggesting strong genetic structuring across the Sierra Nevada. (c) PCA within Central Panama samples, coloured by the site. Compared to samples from the Sierra Nevada, Central Panama samples exhibit a dramatically different pattern, i.e. panmixis, despite a similar spatial and temporal scale of sampling. Colours in (b,c) correspond to geographical locations in (a). (Online version in colour.)
We sequenced 240 regions (each 150–200 bp long) of the Bd genome from the Sierra Nevada and Central Panama skin swab samples. We first conducted a pre-amplification step (which improves performance for amplicon sequencing) and then used these pre-amplified products in a microfluidic PCR approach (see electronic supplementary materials for DNA extraction, preparation and PCR conditions) [42]. Pre-amplified products were loaded into a Fluidigm Access Array IFC, individually barcoded, then pooled for sequencing on ¼ of an Illumina MiSeq lane with 2 × 300 bp paired-end reads at the University of Idaho IBEST Genomics Resources Core. From raw sequence reads, we used the dbcAmplicons software (https://github.com/msettles/dbcAmplicons) to trim adapter, primer sequences and merged to continuous reads. We de-multiplexed and filtered sequences using the reduce_amplicons.R script within the dbcAmplicons repository into two sequence types: ambiguities and raw fastq for each sample. Ambiguities sequence files used IUPAC ambiguity codes to identify multiple alleles. Raw fastq files are all sequences for each sample. Ambiguity sequences were used for phylogenetic analyses and the fastq by the sample was used for alignment, variant calling and filtering VCF for downstream analyses. Specific parameters used for alignment, variant calling and variant filtering can be found in the electronic supplementary material. Post filtering our raw 4534 variants, we recovered 2268 variable sites across 235 amplicons. We applied PCA to examine genetic clustering and structuring among all samples. We estimated PCs using adegenet [43] and visualized in R (v. 3.6.1). We calculated summary diversity statistics using ANGSD [44]. For this analysis, diversity statistics were calculated based on genotype likelihoods, which is distinct from the variant calling approach by Freebayes described above. Given that sample sizes can impact diversity metrics, we randomly subsampled our the Sierra Nevada and global BdGPL samples to equal the number of Central Panama samples (n = 80). Additionally, 49 amplicons (previously developed as Central Panama-specific markers) were removed from the filtered 235 amplicons, leaving 186 amplicons for diversity statistics. Using filtered BAMs from our variant calls, we generated a folded site frequency spectrum given an unknown ancestral state. After estimating the site frequency spectrum for each region, we calculated per-site Watterson's θ and π for the Sierra Nevada (mean amplicon depth = 85.3; s.d. = 78.9; range = 7.0–464.8), Central Panama (mean amplicon depth = 207.8; s.d. = 206.2, range = 15.4–909.6) and global BdGPL (mean amplicon depth = 323.7; s.d. = 383.1; range = 5.5–1342.5) samples. We tested for significant differences in mean Watterson's θ and π and using analysis of variance followed by Tukey's HSD in R (v. 3.6.1), given that we had multiple pairwise comparisons of our global BdGPL reference, Sierra Nevada and Central Panama samples. Using ambiguity sequences by sample, we created a phylogeny including Sierra Nevada, Central Panama and our global BdGPL reference panel. We removed amplicons that had no data and included samples that had least 20 amplicons. We trimmed loci that had greater than 5 bp difference between the minimum and maximum sequence length to control for improper alignments near large indels. A final list of 206 loci were individually aligned using the MUSCLE package in R (v. 3.6.1, [45]) and concatenated (28 688 bp in length). We also included previously published sequences of UM142 BdBrazil as an outgroup [17]. Using this concatenated alignment, we built a phylogeny using IQ-Tree with 1000 ultrafast bootstrap replicates and chose the best model from AIC scores using ‘model-finder’ (GTR + F + I + G4) [46]. We assessed temporal signal in our phylogeny using TempEst (v. 1.5.3) and a date randomization test with TipDatingBeast (v. 1.1) (details available in electronic supplementary material, Methods) [47,48]. After confirming the temporal signal, we inferred time-measured phylogenies, with our concatenated alignment and recorded sampling years, using both BEAST2 [49] and Nextstrain [50]. We used time-measure parameters from a previously published whole-genome phylogenies for Bd to ensure comparability [6]. Briefly, for BEAST2 we used a GTR substitution model with estimated mutation rates 7.29 × 10−7 (lower; 3.41 × 10−7, upper; 1.14 × 10−6) and extended Bayesian skyline plot as demographic parameter [6]. Using this model, we ran a chain which drew samples every 3000 MCMC steps from a total of 575 000 000 steps, after a discarded burn-in of 57 500 000 steps. Convergence of distribution and effective sample size greater than 150 were checked through Tracer (v. 1.7.1) [51]. Our best-supported tree was estimated using maximum clade credibility through TreeAnnotator (v. 2.6) [49] and was visualized using FigTree (v. 1.4.4) (https://github.com/rambaut/figtree). Comparative methods for Nextstrain can be found in the electronic supplementary material. It is important to note that we used BEAST2 and Nextstrain as analytical frameworks to compare patterns between the Sierra Nevada and Central Panama but not to infer exact introduction dates. Applying the same evolutionary models across two geographical regions provides a powerful comparative tool and allows us to infer relative evolutionary rates and introduction timings. However, we interpret specific dates with caution given that patterns of Bd genome evolution may violate a number of model assumptions (e.g. variation across the genome in recombination and mutation rates, variation in chromosomal copy numbers, the potential for both meiotic and mitotic recombination) [19,41], and because our sampling dates do not necessarily correspond to first introduction dates. Given that any violation of basic model assumptions would be shared across study regions, comparisons between the Sierra Nevada and Central Panama can be used to draw conclusions about the relative invasion history in these regions. When comparing within regions, we found significant genetic clustering across the Sierra Nevada (figure 1b) but no genetic clustering across Central Panama (figure 1c). Samples collected from the same site in the Sierra Nevada clustered together, regardless of collection year or species (electronic supplementary material, figure S1). Starting with Unicorn Ponds at the north, samples generally follow a pattern of isolation by distance. LeConte Divide and Conness Pond are somewhat anomalous, however, because they overlap in PC space but are geographically separated by approximately 80 km (figure 1b). By contrast, Central Panama genotypes exhibited panmictic patterns, regardless of locality, collection year or species, indicating no genetic structuring across a similar spatio-temporal scale (figure 1c, electronic supplementary material, figure S2). We confirmed that Bd from the Sierra Nevada and Central Panama belong to the global BdGPL lineage. Both Sierra Nevada and Central Panama samples clustered with a panel of global samples that were previously identified as BdGPL [17]. Interestingly, the Sierra Nevada and Central Panama samples clustered separately from each other in PC space when compared to global BdGPL samples (figure 2a). Additionally, we found that overall genetic diversity was significantly higher in the Sierra Nevada as compared to Central Panama [Tukey HSD, p < 0.0001] (figure 2b,c). Remarkably, we also found that Sierra Nevada Bd samples have comparable and, in the case of Watterson's θ, higher diversity than the set of global BdGPL samples [Tukey HSD, p < 0.0001]. When comparing Central Panama and the Sierra Nevada using individual sites with similar samples sizes, we found that the majority of Sierra Nevada sites had higher mean diversity compared to Central Panama sites (both Watterson's θ and π) [Tukey HSD, p < 0.001], except in the lowest sample size pairing (N = 5) where El Valle S. had significantly higher mean diversity than LeConte Divide (electronic supplementary material, figure S3). Figure 2. Genetic differentiation and diversity among the Sierra Nevada, Central Panama, and global BdGPL samples. (a) PCA based on BdGPL genotypes from the Sierra Nevada (n = 130), Central Panama (n = 80) and global reference panel (n = 120). Colours indicate samples from each region. The global reference panel included samples from dozens of frog species across all continents with BdGPL. Samples from the Sierra Nevada and Central Panama are almost entirely separated in PC space with the Sierra Nevada samples showing greater genetic variation than Central Panama samples. (b) Distribution of mean genetic diversity (Watterson's θ) for all variable sites based on region. Samples from the Sierra Nevada and global panels were randomly subsampled to match Central Panama sample size (all regions n = 80). Mean genetic diversity was significantly higher for Sierra Nevada samples compared to Central Panama samples and to the global BdGPL panel [Tukey HSD, p < 0.0001]. (c) Distribution of mean nucleotide diversity (π) for all variable sites based on the region using the same samples as (b). Mean nucleotide diversity was significantly lower for Central Panama samples compared to Sierra Nevada samples and the global BdGPL panel [Tukey HSD, both p < 0.0001]. Each box plot shows the median (horizontal line), first and third quartiles (bottom and top of box, ‘hinges’), lowest and highest values within inter-quartile range of the lower and upper hinges (vertical lines) and outliers (points). (Online version in colour.)
Using a time-dated phylogenetic approach that included previously published global BdGPL samples for reference [17], we found branches from Sierra Nevada samples were comparatively older than those in Central Panama (figure 3, electronic supplementary material, figure S4). As discussed in the Material and methods, we do not assume the specific inferred dates are accurate given the likelihood that dynamics of Bd genome evolution violate several model assumptions. Our root-to-tip regression showed somewhat low temporal signal in our data (R2 = 0.02) (electronic supplementary material, figure S5). However, our date randomization tests showed no overlap between real and randomized datasets, indicating a sufficient level of the temporal signal (electronic supplementary material, figure S6). While these results may appear contradictory, root-to-tip regression is a conservative approach assuming a strict molecular clock [47], while date randomization provides a more statistically robust method of comparison [48]. Therefore, as discussed in the Methods section above, our time-dated approaches are appropriate for inferring relative invasion histories across regions rather than proposing specific divergence or invasion dates. Figure 3. BEAST2 timed dated phylogeny among the Sierra Nevada, Central Panama, and global BdGPL samples. Branch tips are colour coded by region. The tree is rooted by an outgroup from a more basal Bd lineage (BdBrazil isolate UM142). Sierra Nevada samples are found across the tree, in multiple clusters, and with longer branch lengths than Central Panama samples suggesting a longer history of Bd in this region. (Online version in colour.)
For BEAST2, the time to most recent common ancestor (tMRCA) for Sierra Nevada samples was estimated to be 474 years from present day (95%HPD 510–393 years from present day) and estimated tMRCA in Central Panama was 277 years from present day (95%HPD 389–60 years from present day) (figure 3). For Nextstrain, tMRCA for Sierra Nevada samples was estimated as 1407 years from present day (95% CI 4498–1151) and tMRCA for Central Panama was estimated as 666 years from present day (95% CI 1914–534) (electronic supplementary material, figure S4); dynamic Nextstrain visualizations are available at: https://nextstrain.org/community/andrew-rothstein/bd-gpl/auspice/viz. Therefore, even without ascribing weight to specific inferred dates, Bd in the Sierra Nevada appears to be much older than Bd in Panama. Confidence intervals for the inferred tMRCA do not overlap between regions with either analysis. The BEAST2 and Nexstrain time-dated phylogenetic approaches also corroborated PCA results (figure 1). Sierra Nevada samples largely clustered by the site while Central Panama samples had little to no structure based on-site location. (electronic supplementary material, figure S4) Finally, phylogenetic trees show an expected split within BdGPL, supported by high BEAST posterior node values (electronic supplementary material, figure S7). This correspond to a previously reported split separating BdGPL into two subclades: BdGPL-1 and BdGPL-2 [20,21]. Only GPL-2 is represented in Panama samples while GPL-1 and GPL-2 are both found in the Sierra Nevada samples. Bd has caused mass amphibian declines in many regions of the world and has been an example of the devasting effects of emerging wildlife diseases [12,15,16,35,52–54]. However, assessments of Bd emergence and spread have yet to incorporate genetically informed epizoology to examine pathogen dynamics at fine spatial scales. Our study used comparative population genetics to examine the genetic signatures of BdGPL across two emblematic regions with disease-related amphibian declines. The alpine lakes of the Sierra Nevada and the tropical forests of Central Panama have dramatically different climate, habitat and host communities. However, they have been described as having similar histories of recent Bd emergence and spread. We tested the assumption that BdGPL was recently introduced to these two regions and swept through each in a unidirectional epizootic wave. We found dramatic differences in Bd evolutionary history across regions, with an unexpectedly deep history of Bd in the Sierra Nevada. Here, we explore differences across regions providing a new perspective on these important historic declines. As wildlife disease rapidly continue to spread across the world, our framework is broadly applicable to interrogating observed patterns of pathogen emergence and spread to uncover important evolutionary pathogen histories. Our results from Central Panama support the hypothesis of a recent introduction, with Bd in this region lacking any spatial structure. All Bd genotypes from the Central Panama group tightly together, are generally distinct from Bd collected in the Sierra Nevada, and are all part of the GPL-2 subclade. This pattern supports previous studies reporting a single fast-moving outbreak of Bd through Central Panama [35]. Our samples from Central Panama were collected approximately 8 years after observed outbreaks (between 2012 and 2016), and the observed lack of genetic structure indicates that Bd did not diverge on a site-specific basis over this time period. Our findings support other recent studies showing a lack of genetic, phenotypic and functional shifts in Central Panama Bd across similar temporal scales [32]. BdGPL appears to have arrived in Panama much more recently than in the Sierra Nevada, maintained low levels of genetic diversity, and, over the last two decades, currently has no detectable genetic sub-structure. We observed a dramatically different pattern in the Sierra Nevada, where we found high levels of genetic variation between sampling sites and spatial structuring of Bd genotypes. Although Bd samples were collected across a similar spatial and temporal scale as those from Panama, our genetic data indicate that BdGPL has likely had a much longer historical presence in the Sierra Nevada than it has in Panama. This conclusion is supported by multiple lines of evidence. First, Sierra Nevada Bd contains more genetic variation and diversity than Central Panama (figure 2a). Measures of nucleotide diversity (π), are higher in Sierra Nevada Bd samples compared to Central Panama and Sierra Nevada Bd genetic diversity (Watterson's θ) is significantly higher than the entire global panel of BdGPL samples (figure 2b). This result is consistent with previous evidence that BdGPL in the Sierra Nevada has higher levels of genetic diversity than BdGPL from Arizona, Mexico, or Central Panama [55]. Second, we also observed a surprising pattern of spatially structured genetic diversity for BdGPL in the Sierra Nevada. Sierra Nevada BdGPL genotypes typically cluster by site and segregate by geographical distance in PC space and in the phylogeny (figures 1b and 3b). Much of the observed genetic structure in the Sierra Nevada is consistent with a pattern of isolation by distance, suggesting a much longer history of Bd on the landscape. Third, even the exceptions to the pattern of isolation by distance suggests a deeper and more complex history of Bd in the Sierra Nevada. Samples from LeConte Divide and Conness Pond are genetically distinct from all other samples in the Sierra Nevada and cluster in PC space (figure 1b). These samples belong to a separate, early branching clade referred to as GPL-1 (figure 3). The presence of both BdGPL-1 and BdGPL-2 subclades could represent multiple independent introductions or much deeper in situ divergence, possibilities we revisit below. One key factor that could contribute to radically different patterns of Bd genetic variation between Central Panama and the Sierra Nevada is invasion history (the timing and number of introductions). Our Nextstrain and BEAST2 analyses infer that Bd from the Sierra Nevada is older than Bd from Central Panama (figure 3, electronic supplementary material, figure S4). While our inference indicates that BdGPL has been in the Sierra Nevada longer than Central Panama, it is difficult to assert specific invasion dates. As discussed in the Material and methods section, patterns of Bd genome evolution may violate a number of model assumptions. Although our analyses used a species-specific mutation rate inferred from Bd whole-genome analyses [6] our assay targets regions of the Bd genome that are most informative for discriminating among Bd lineages [42] and therefore may not evolve with a shared background mutation rate. Even without specific introduction dates, studies using histology and qPCR to test for Bd presence in museum specimens have often shown Bd presence prior to field-observed die-offs [24,56,57], which could indicate older introduction timings than previously assumed. As such, Bd presence has been detected in samples as far back as 1932 in the Sierra Nevada [24] and 1964 in Costa Rica (adjacent to Panama) [56]. Moreover, field observations suggest that Bd may be present in the environment well before an outbreak is observed. In some lakes, Bd is present at almost undetectably low prevalence and load for years before Bd loads spike and die-offs occur [16,28,58]. In some systems, Bd can even be detected from eDNA surveys before die-offs occur [59]. Such dynamics challenge our a priori expectations that Bd die-offs occur immediately after the pathogen first arrives in an area. In some systems, such as the Sierra Nevada and parts of Costa Rica [24,56,57], it is possible that Bd had a more wide-spread presence earlier than perceived. Whether there actually were earlier Bd-caused die-offs remains an open question. Increased surveillance of Bd before and during early outbreaks is needed to decouple initial pathogen invasion from observed pathogen-induced declines. High levels of genetic variation, deep spatial genetic structure and the presence of both subclades of BdGPL in the Sierra Nevada suggest a longer evolutionary history of Bd in the region than previously appreciated. The presence of both BdGPL-1 and BdGPL-2 could represent multiple asynchronous invasions of BdGPL, a hypothesis raised by another recent spatial–temporal study of Bd presence in the Sierra Nevada [24]. An alternative explanation is that California is a potential source of Bd that has spread to other regions. As sampling resolution improves, it is possible that we will find other regions of the world with highly diverse and spatially structured BdGPL populations. However, it is also worth continuing to challenge our assumptions about the origin and spread of this lineage. While the most basal lineage of Bd is from Asia [6], the origin of BdGPL remains highly uncertain. Although we often assume that BdGPL presence results from recent invasions, the region from which BdGPL originated would be expected to have general characteristics similar to what we observe in the Sierra Nevada (i.e. relatively high genetic diversity and deep spatial structure). No such region other than the Sierra Nevada has yet been identified. Global sampling with greater spatial and temporal resolution will be needed to ultimately determine the origins of this highly virulent Bd lineage. Biotic and abiotic factors also likely influence patterns of Bd genetic variation in a consistent direction, with increased opportunity for pathogen mixing in Central Panama relative to the Sierra Nevada. Central Panama is home to a diverse amphibian assemblage, with dozens of sympatric species that use a variety of microhabitats and have different reproductive modes [15,31]. A diverse host community in Panama with year-round activity and some direct-developing species (i.e. those without an aquatic larval phase) could provide more opportunities for Bd spread [31,60]. Central Panama contains landscape features that may be barriers to dispersal for some amphibian species [61], but interconnected stream networks still allow for fairly high connectivity among sites. By contrast, in the Sierra Nevada, our samples are from the only common—and highly susceptible—amphibian species in the alpine lake habitats (R. sierrae/muscosa) [62]. Rana sierrae/muscosa have high site fidelity, limited overland movements, spend the majority of each year under ice and inhabit disjunct alpine lakes separated by high mountain passes [63]. These features all impede connectivity among host populations and provide fewer opportunities for Bd dispersal [64]. Therefore, landscape and host factors consistently provide decreased opportunities for Bd gene flow in the Sierra Nevada, which is reflected in the greater pathogen spatial structure in this region. In addition, Central Panama is significantly warmer and wetter than the Sierra Nevada. Temperature differences are particularly important because warmer temperatures (to a point) can lead to faster pathogen growth, an increased number of generations per year and greater opportunity for rapid evolutionary change. Slower Bd growth, generation time and evolutionary rates in the Sierra Nevada compared to Central Panama, make the patterns of higher genetic diversity and strong spatial genetic structure in the Sierra Nevada all the more interesting. Ultimately, integrating genetic, spatial and epizootic data within an evolutionary framework is a powerful way to understand the dynamics of emerging diseases of wildlife. Typically, studies of wildlife disease dynamics rely on a priori assumptions about pathogen introductions (i.e. based on earliest infection known from wild populations or museum records). However, our results, using Bd as an example for global pathogens, clearly demonstrates that outbreaks with similar epizoological signatures can still have radically different underlying pathogen histories. In our study, two regions with similar observed epizoological patterns in the field exhibit dramatically different pathogen evolutionary histories. In fact, one of the regions—the Sierra Nevada—has considerable pathogen diversity and genetic structure. Supporting evidence suggests that Bd in this region may persist in populations of highly susceptible host species at very low levels over many years without causing epizootics, opening the possibility that the pathogen has a much longer evolutionary history than previously appreciated. When we treat all population declines as the same, we overlook important nuances that could assist on-the-ground recovery and mitigation efforts. For example, if we incorporate Bd genotype data into choices of donor frog populations when planning translocations and reintroductions, we can mitigate human-induced mixing of Bd genotypes. Such actions could be an important component for species recovery efforts. By combining genetic and epizoological data, we can better understand differences in pathogen invasion history across regions and support more effective policies for biodiversity conservation and management. All sample collections were authorized by research permits provided by NPS, USFWS and the Institutional Animal Care and Use Committees of UC Berkeley, UC Santa Barbara, University of Nevada, Reno and the University of Pittsburgh. Processed data including BEAST2 tree XML, VCF file, genetic diversity values, tests for temporal signal and code for figures: https://figshare.com/s/0b3bcabff81fae2fb8e8. Raw sequences of new samples deposited in NCBI SRA BioProject ID PRJNA686993: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA686993. Nextstrain code: https://github.com/andrew-rothstein/bd-gpl.git. The data are provided in the electronic supplementary material [65]. A.P.R.: conceptualization, data curation, formal analysis, investigation, methodology, visualization, writing—original draft, writing—review and editing; A.Q.B.: conceptualization, data curation, formal analysis, investigation, methodology, visualization, writing—original draft, writing—review and editing; R.A.K.: conceptualization, data curation, funding acquisition, investigation, methodology, project administration, supervision, visualization, writing—original draft, writing—review and editing; C.J.B.: conceptualization, data curation, funding acquisition, investigation, project administration, supervision, writing—original draft, writing—review and editing; J.V.: conceptualization, data curation, funding acquisition, investigation, project administration, supervision, writing—original draft, writing—review and editing; C.L.R.: conceptualization, data curation, funding acquisition, investigation, project administration, supervision, writing—original draft, writing—review and editing; E.B.R.: conceptualization, data curation, funding acquisition, investigation, methodology, project administration, supervision, visualization, writing—original draft, writing—review and editing. All authors gave final approval for publication and agreed to be held accountable for the work performed therein. We declare we have no competing interests. This work was supported by DoD SERDP contract RC-2638 (to E.B.R., C.L.R.Z., J.V., C.J.B.), NSF DEB-1557190 (to E.B.R., C.J.B., R.A.K.), NSF DEB-1551488 (to E.B.R., C.L.R.Z., J.V.), NSF DEB CAREER - 1846403 (to J.V.), NSF DEB-166311 (to C.L.R.Z.), and the National Park Service. Sequencing done at IBEST Genomics Resources Core at the University of Idaho was supported in part by NIH COBRE grant P30GM103324. All sample collections were authorized by research permits provided by NPS, USFWS and the Institutional Animal Care and Use Committees of UC Berkeley, UC Santa Barbara, University of Nevada—Reno and University of Pittsburgh. We thank the Sierra Nevada Aquatic Research Laboratory field crew, Danny Boiano from the National Park Service, Matt Robak, Angie Estrada, Renwei Chen and Mary Toothman for assisting in field collection and pathogen qPCR. FootnotesElectronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5459508. © 2021 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 4Stony corals secrete calcium carbonate skeletons which create three-dimensional reef frameworks that support the most productive and biologically diverse marine ecosystems on Earth. Corals have a biphasic life cycle with swimming planktonic larvae and sessile adults [1]. Successful settlement of larvae and their subsequent survival (known as ‘recruitment’) is a key element for coral reef resilience and recovery following major disturbances [2]. Young coral settlers are highly vulnerable to mortality resulting from environmental changes, such as rising seawater temperature and ocean acidification (OA) associated with global climate change [3]. OA is expected to negatively affect mineral formation of marine calcifiers, such as corals, by making it more energetically difficult for these organisms to deposit their skeletons [4]. This circumstance has led to numerous studies on the effects of OA on coral calcification, that is predicted to substantially decline over the course of this century [5]. However, great uncertainty still remains about the extent of the threat of OA to coral persistence. This is partly due to the high variability of the coral response to OA and to the persistence of coral communities that thrive in naturally low pH environments [6,7], which has shed light on the potential for coral acclimation and adaptation to the projected acidification scenarios. Most of the studies conducted on primary polyps focused primarily on changes occurring in the skeletal development, which has been shown to be significantly delayed, with the potential for substantial deformities with decreases in seawater pH [8]. Studies on the functioning of the coral–algal symbiosis under OA are scarce, especially related to coral early-life stages. Coral symbiotic algae, which belong to the family Symbiodiniaceae, have a major role in the coral holobiont physiology and nutrition [9]. These dinoflagellates supply the coral host with photosynthetic products, thereby supporting the host metabolism, growth and reproduction [9]. A comprehensive understanding of how the earliest stages of coral, together with their algal partners, will respond to OA is critically needed to assess the capacity of these organisms to adapt and/or acclimatize to a future acidified ocean. To this end, we performed a controlled laboratory study by culturing larvae and primary polyps (metamorphosed and settled larvae) of the coral Stylophora pistillata under reduced pH values (high CO2 concentrations) predicted by the end of the century, specifically at pH 8.2, pH 7.8 and pH 7.6 [10]. Here, we investigate larval ability to survive and successfully recruit under OA conditions, and we examine key physiological processes and skeletal characteristics influenced by changes in seawater pH, at the level of both phenotype and gene expression. Our findings reveal that the response of young corals to OA involves an intricate network of interrelated changes in the coral holobiont. Coral larvae were collected from 20 randomly selected adult colonies of the stony coral S. pistillata on the reef adjacent to the Interuniversity Institute of Marine Sciences (IUI, 29°30′06.0″ N 34°54′58.3″ E) in the Gulf of Eilat (Israel), under a special permit from the Israeli Natural Parks Authority. All larvae were pooled together (approx. 50 larvae from each parental colony) and were transported to a controlled environment aquarium system at the Leon H. Charney School of Marine Science at the University of Haifa. The carbonate chemistry of seawater was manipulated in the experimental aquariums by injecting CO2 to reduce the ambient pH (pH 8.2) and obtain the target values of pH 7.8 and pH 7.6. After 9 days of experimental treatment exposure, the numbers of settled primary polyps, still swimming larvae (not settled) and dead larvae (larvae dissolved during the experiment) were recorded. The primary polyps were gently removed from the chambers and they were collected for subsequent physiological, skeletal and molecular analysis (for full details, see electronic supplementary material, Methods). Following the 9 days experimental exposure to decreased pH conditions, we observed a reduction in the percentage of settled larvae from 82% at the control pH 8.2, to 64% at pH 7.8, and to 46% at pH 7.6, which resulted to be significantly lower than the control (Fisher's exact test, p < 0.0001) (electronic supplementary material, figure S2). At the intermediate pH, we detected a significantly higher percentage of larvae that were still at the planktonic or swimming stage relative to pH 8.2 (Fisher's exact test, p < 0.0001), while at the lowest pH, we found significantly higher larval mortality compared to the control (Fisher's exact test, p < 0.0001). We then explored the physiological changes that allowed the surviving individuals to successfully settle and develop into a primary polyp. Measurements of the primary polyps dark-respiration rates, indicative of the coral metabolic rates, were significantly reduced at the lowest pH condition, with an average decrease of 64 and 69% at pH 7.6 relative to pH 8.2 and pH 7.8, respectively (one-way ANOVA, N = 6, F2,15 = 5.927, p = 0.01) (electronic supplementary material, figure S3A). In contrast with the reduction in respiration rate, coral host protein concentration, used as a proxy for tissue biomass, significantly increased in both pH 7.8 and pH 7.6 compared to the control pH 8.2 (electronic supplementary material, figure S3B) (one-way ANOVA, N = 3, F2,24 = 24.99, p < 0.0001). We then explored changes in the algal physiological attributes, and we found a significantly higher number of endosymbiotic algae per polyp surface area at both pH 7.8 and pH 7.6 relative to pH 8.2 (electronic supplementary material, figure S3C) (one-way ANOVA, N = 3, F2,24 = 99.25, p < 0.05), pointing to an enhancement of algal growth with decreasing pH. The number of symbiont cells per host protein did not show instead any significant change (electronic supplementary material, figure S3D). Furthermore, chlorophyll a concentration was higher in the polyps at the lowest pH (one-way ANOVA, N = 3, F2,24 = 73.41, p < 0.0001) (electronic supplementary material, figure S3E), even though its concentration within a single algal cell did not significantly vary relative to the control pH (electronic supplementary material, figure S3F), providing further evidence of the augmented algal density at pH 7.6. In addition, higher chlorophyll auto-fluorescence was detected in the live primary polyps at the lowest pH (one-way ANOVA, F2,6 = 73.38, p > 0.0001) (electronic supplementary material, figure S4A,B), matching the chlorophyll a concentration results (electronic supplementary material, figure S3E). The observed changes in the algal physiological properties are indicative of variations in photosynthetic traits. Indeed, we observed an increase of photochemical efficiency and activity in the endosymbiotic algae at pH 7.6. The maximal quantum yield (Fv/Fm) of photosystem II (PSII) was significantly higher in the primary polyps at pH 7.6 relative to the other two pH (electronic supplementary material, figure S4C) (Kruskal–Wallis test, H = 7.199, p = 0.027). In addition, although the measurements of the initial slope (α) and of the minimal photoinhibition point (Ek) showed no significant change (electronic supplementary material, table S5), the maximum values of the relative electron transport rate (rETR), were significantly elevated in the polyps at pH 7.6 relative to the control (electronic supplementary material, figure S4D) (Kruskal–Wallis test, H = 8.025, p = 0.018). Taken together, these results show that the observed increase of algal density under acidification conditions is paralleled by the enhancement of photosynthetic efficiency and activity. We also measured the maximum values of the non-photochemical quenching (NPQ) and found no significant change among pH treatments (electronic supplementary material, table S5). The existence of a linkage between tissue biomass and growth [11], and our observation of increased host tissue at reduced pH, prompted us to explore changes in the primary polyps skeletal characteristics. By employing the fluorescent dye calcein [12], we visualized patterns of incorporation of divalent ions (such as Ca2+) into the live primary polyps and subsequently into the skeleton. The in vivo imaging shows that in all three pH conditions, the highest calcein fluorescence intensity is found in the calyx area (see electronic supplementary material, figure S6 for additional details) and in the skeletal septa (electronic supplementary material, figure S5A). In corals, ions and other particles are ingested from seawater into the mouth cavity and are then transported to the site of calcification [13]. Our observations indicate that in all three pH conditions, corals are actively up-taking Ca2+ from the mouth, and that these ions are eventually incorporated into the skeleton. Measurements of the calcein fluorescence intensity of the primary polyps skeleton (N = 3 polyps per pH condition, only 1 polyp per pH is shown in electronic supplementary material, figure S5A) show a significantly lower calcein fluorescence at pH 7.8 and 7.6 compared to pH 8.2 (electronic supplementary material, figure S5B) (one-way ANOVA, F2,6 = 213.5, p < 0.0001). This indicates that a lower amount of Ca2+ was incorporated into the skeleton, thus pointing to a lower skeletal development at reduced pH conditions. Patterns of calcein incorporation correspond to the skeletal morphology imaged using micro-CT. Slices of the 3D reconstructions of the polyps skeleton (figure 1d–f) show a progressive reduction in the thickness of the skeletal septa as the pH decreases (figure 1g–l). Such reduction, that was significant at both pH 7.8 and 7.6 compared to the control (one-way ANOVA, F2,15 = 24.32, p > 0.001) (electronic supplementary material, figure S5C), indicates that at reduced pH, a lower amount of CaCO3 was deposited by the coral. Figure 1. Changes in skeletal growth patterns in coral primary polyps with decreasing pH. (a–c) Microscopy images showing the Calcein Blue fluorescence (expressed as percentage intensity) in the skeleton of the primary polyps at (a) pH 8.2, (b) pH 7.8 and (c) pH 7.6. Magnification: 4×. Scale bar: 400 µm. (d–f) Micro-CT images showing top views of three-dimensional reconstructions of the primary polyps' skeleton at (d) pH 8.2, (e) pH 7.8 and (f) pH 7.6. Scale bar: 400 µm. (g–i) Slices of the three-dimensional reconstructions showing thresholded two-dimensional slices of the base of the septa at (g) pH 8.2, (h) pH 7.8 and (i) pH 7.6. Scale bar: 200 µm. (j–l) Skeletal thickness distribution (μm) along the vertical axis of the septa of the primary polyps at (j) pH 8.2, (k) pH 7.8 and (l) pH 7.6. Scale bar: 200 µm. (Online version in colour.)
Further changes in skeletal growth patterns were additionally detected at a micro-scale. Measurements of the calyx area (figure 2a–c, cyan overlays) of the polyps imaged by scanning electron microscopy indicate that there was no significant difference between calyxes among the three pH conditions (figure 2m). Differently from the calyxes, the crown areas shown in figure 2a–c (yellow overlays), corresponding to the forming coenosteum (see electronic supplementary material, figure S6), were significantly smaller in the primary polyps at pH 7.8 and 7.6 relative to pH 8.2 (one-way ANOVA, N = 6, F2,15 = 8.17, p < 0.01). At higher magnification, the thickness of the septa (electronic supplementary material, figure S6) was significantly smaller in polyps grown at lower pH compared to the control (figure 2d–f,n) (one-way ANOVA, N = 27, F2,78 = 51.92, p < 0.0001). Overall, these observations show that the reduction of skeletal development detected at the macro-scale is also manifested at the microstructural level. Figure 2. Modifications of skeletal features of coral primary polyps at different pH conditions. (a–c) SEM images showing the micro-morphology of the calyx (cyan areas) and crown (orange areas) of the primary polyps at (a) pH 8.2, (b) pH 7.8 and (c) pH 7.6. Enlargements of the skeletal (d–f) septa and (g–i) spines at (d,g) pH 8.2, (e,h) pH 7.8 and (f,i) pH 7.6, respectively. (g) White arrows indicate the RADs on the spine. (j–l) Insets showing the RADs surface texture at (j) pH 8.2, (k) pH 7.8 and (l) pH 7.6. (m) Size of the calyx area (cyan bar) and the crown area (orange bar) calculated in the primary polyps (N = 6 per pH condition). (n) Thickness of the septa calculated in the primary polyps at (d) pH 8.2, (e) pH 7.8 and (f) pH 7.6 (N = 27 per pH condition). (o) RADs number measured per units of basal area of the spines (mm) at (g) pH 8.2, (h) pH 7.8 and (i) pH 7.6 (N = 9 per pH condition). (m–o) Asterisks (*) indicate statistical differences (p < 0.05, one-way ANOVA) relative to the control pH 8.2. Scale bars: (a–c) 200 µm, (d–f) 100 µm, (g–i) 10 µm and (j–l) 1 µm.
Moreover, the terminal portion of the skeletal spines located on the septa showed a significantly lower number of rapid accretion deposits (RADs; globular elements in figure 2g–i; see electronic supplementary material, figure S6) in the polyps at reduced pH, when compared with the control (figure 2o) (one-way ANOVA, N = 9, F2,24 = 19.72, p < 0.0001). RADs correspond to areas of the skeleton with rapid CaCO3 deposition [14]. Hence, the lower abundance of RADs in the polyps at pH 7.8 and 7.6 indicates a reduced development of these regions in acidified seawater. In addition, we found prominent differences in the RADs texture. At the control pH, the texture of the RADs was smoother and more compact in comparison to the acidified conditions (figure 2j–l). At pH 7.6, the aragonite bundles of fibres, corresponding to the microcrystalline features on the RADs surface (figure 2j–l), were less compact and were characterized by granular aggregates forming needle-like structures with differing shapes and orientations (figure 2l) compared to the texture at both pH 8.2 and 7.8. This microstructural pattern conferred a more porous appearance to the RADs surface. The exposure to acidified seawater conditions, therefore, generates changes at multiple levels of the skeletal hierarchy, from macro- to micro-scale. We examined the molecular controls underlying the observed modifications in the primary polyps physiological and skeletal characteristics. Analysis of the coral host's differentially expressed genes (DEGs) show that, compared to pH 8.2, a higher number of DEGs was found at pH 7.6 as opposed to pH 7.8 (figure 3a). This pattern, which reflects the logarithmic nature of the pH scale, indicates that the slope of the progression in DEGs with decreasing pH sharply steepens between pH 7.8 and 7.6, where larger transcriptional modifications are needed to adjust to acidified seawater. Figure 3. Coral DEGs, functional enrichment and biomineralization-related genes expression between experimental pH conditions. (a) Number of DEGs in coral primary polyps detected, respectively, in the pH 7.8 treatment relative to the control pH 8.2 (left bar), in the pH 7.6 treatment relative to pH 8.2 (middle bar) and in the pH 7.6 treatment relative to pH 7.8 (right bar). Red indicates the upregulated genes and blue indicates the downregulated genes within each comparison. (b) Dendrogram on the left showing enriched GO terms and KEGG pathways clustered according to the portion of identical genes shared. The heat map shows the positive (upregulated genes) or negative (downregulated genes) percentages of DEGs per GO term or KEGG pathway, detected, respectively, at pH 7.6 relative to pH 8.2 (right column) and to pH 7.8 (left column). Red numbers indicate the total number of genes within each per GO term or KEGG pathway. DE (%), percentage of DEGs. (c) Heat map showing the log2 fold change of significantly (adjusted p-value > 0.05) upregulated (red) and downregulated (blue) biomineralization-related genes at pH 7.8 and 7.6 when compared with pH 8.2. (Online version in colour.)
Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed to identify enriched genes groups and molecular pathways involved in coral response to pH reduction. Figure 3b shows that among all pH treatments comparisons, the largest number of enriched genes groups (both up- and downregulated) was found between pH 7.6 and pH 8.2. The upregulated gene groups at the lowest pH show a significant enrichment for genes associated with the activity of voltage-gated calcium channels, with cell differentiation and developmental processes, and with genes linked to sensory perception and detection of environmental stimuli (figure 3b). The downregulated genes groups at pH 7.6 include a large number of genes involved in lipid metabolism, mitochondria biosynthesis and activity, in energy production and respiration. The downregulation of these genes points to a reduced energy metabolism at acidified conditions, disclosing the molecular control behind the observed respiration rate decrease (electronic supplementary material, figure S3A). Looking into the molecular mechanisms underlying the observed changes in skeletal features, we evaluated the differential expression of putative coral biomineralization-related genes (electronic supplementary material, table S4). Among these 91 total genes, 5 were differentially expressed in the polyps at pH 7.8, and 23 were differentially expressed in the polyps at pH 7.6 (figure 3c), indicating a greater response of the biomineralization molecular machinery at the extreme pH. In particular at pH 7.6 relative to pH 8.2, a significantly higher number of biomineralization-related genes was upregulated, compared to the number of downregulated genes (Fisher's exact test, p < 0.05). Most of these upregulated genes code for proteins of the skeletal organic matrix with structural and adhesion functions. These proteins are embedded in the skeletal framework and control the growth of the aragonite crystal-like fibres within the coral skeleton [15]. Notably, we observed the downregulation of the coral acid-rich protein 4 (CARP4) gene at the lowest pH treatment (figure 3c). As member of the CARPs group, this biomineralization toolkit protein is expected to transport Ca2+ to the calcification site that are then concentrated in centres of calcification, also known as RADs [15]. Moreover, among the upregulated genes, we detected the S. pistillata carbonic anhydrase 2 (STPCA-2) and the L-type calcium channel (figure 3c), which is located in the calicoblastic epithelium (cell layer overlying the coral skeleton) [16]. These two genes are involved in regulating the pH and carbonate chemistry of the calcifying medium, the confined space where skeletal deposition occurs [17,18]. In particular, the increased expression of the L-type calcium channel, together with the enrichment of genes associated with calcium channels activity (figure 3b), indicates that the exposure to acidic seawater stimulates in the corals the transport of Ca2+. In S. pistillata algal endosymbiont, enriched genes groups (both up- and downregulated) were only found in pH 7.6 relative to the other two pH conditions, while no enrichment was detected in pH 7.8 relative to the control (figure 4). This indicates that, similar to the coral host, wider transcriptional changes take place at the lowest pH. Figure 4. Changes in expression levels of different functional groups of genes in the coral endosymbiotic algae. Dendrogram showing enriched functional groups clustered according to the portion of identical genes shared. The heat map shows the positive (upregulated genes) or negative (downregulated genes) normalized enrichment score (obtained using the functional enrichment analysis GSEA) per each functional group, detected, respectively, at pH 7.6 compared to pH 8.2 (right bar) and to pH 7.8 (left bar). No enrichment was found in pH 7.8 compared to pH 8.2. Genes groups are shown only for cases with significant (adjusted p-value < 0.05) enrichment in both comparisons. Red numbers indicate the total number of genes within each functional group. NES, normalized enrichment score. (Online version in colour.)
The downregulated groups at pH 7.6 include genes related to photosynthesis (e.g. ‘thylakoid’), gene expression (e.g. ‘ribosome’) and proton transport (e.g. ‘proton transmembrane transport’). The upregulated groups comprise a large number of genes involved in photosynthetic carbon fixation (‘ribulose-bisphosphate carboxylate activity’, also known as RuBisCO, ‘carboxy-lyase activity’, etc.), carbohydrates binding and transport (e.g. ‘polysaccharide binding’, ‘substrate-specific channel activity’, ‘active transmembrane transporter activity’). As shown by the clustering (figure 4), terms related to carbohydrates binding and transport share a portion of common genes, suggesting that symbiont genes associated with transmembrane transport are also involved in the transport of sugars. Taken together, our observations indicate that the symbiotic algae response to reduced pH is governed by the modulation of multiple photosynthesis-related genes. The current study on S. pistillata early-life stages reveals that the exposure to predicted future OA conditions leads to a reduction of larval survival and settlement (electronic supplementary material, figure S2). Nonetheless, the more acidic seawater did not entirely inhibit larval recruitment, suggesting potential acclimatory mechanisms that allowed the corals to successfully settle and start growing. Larval recognition of the settlement substrates occurs through environmental cues, that are exogenous factors capable of inducing larval settlement [19]. Seawater acidification was shown to reduce the induction of larval settlement in response to environmental cues [20]. Based on our transcriptomic analysis, genes shown to be related to larval selection of and attachment to the substrate (i.e. sensory perception and detection of environmental stimuli-related genes) [21,22] are significantly enriched at the lowest pH (figure 3b). This indicates that the modulation of genes with a role in sensory perception may account for the adjustment capacity and survival of the corals under OA conditions, as observed in other organisms exposed to a changing environment [23]. Diminished availability of settling juveniles could inhibit the replenishment of reefs after sporadic disturbances such as storms and bleaching events [2], with the potential to compromise coral reef resilience. However, coral populations can naturally persist in acidic environments, despite having a lower recruitment efficiency in comparison to corals living at ambient pH [24]. Moreover, improved adult performance could compensate for low recruitment rates, as observed for corals living around CO2 seeps [25]. In this study, larvae reared under acidified conditions that successfully metamorphosed into a primary polyp were characterized by reduced metabolic rates, as shown by the lower respiration rate (electronic supplementary material, figure S3A) and the downregulation of metabolism-related genes (figure 3b). A downward trend of energy metabolism seems a common response in young corals exposed to acidic seawater [26]. In an OA scenario, the formation of the skeleton is thought to be more energetically expensive for corals, requiring more energy to remove from the calcifying space the excess H+ produced during calcification [4]. The reduction of metabolic rates observed in the primary polyps may energetically limit the rates of H+ removal. This leads to a diminished ability to deposit CaCO3, and indeed our multifaceted examination of the skeleton reveals that under OA primary polyps were characterized by reduced skeletal development (figures 1a–c,j–l and 2a–i,m–o; electronic supplementary material, figure S5B,C). Furthermore, our observations of the reduction in primary polyp crown areas (figure 2a–c,m) and in the number of RADs (figure 2g–i,o) show that exposure to OA conditions involves changes within the skeleton that have never been investigated so far. In this context, the downregulation of the CARP4 gene detected at acidified conditions (figure 3c) reveals a potentially lower contribution of this protein to the formation of RADs, thereby contributing to the overall reduction in skeletal growth. In fact, this biomineralization toolkit protein was localized in S. pistillata skeleton in correspondence of RADs, and was suggested to form crystal binding substrates that lead to CaCO3 nucleation [15]. Given the reduction of skeletal development in acidified seawater, corals must put in place alternative strategies to sustain overall growth. Indeed, the greater amount of tissue biomass covering the skeleton, detected in this study (electronic supplementary material, figure S3B), could allow coral polyps to reach larger sizes, representing a critical strategy to counterbalance the lower skeletal development. Through a lower investment of energy into skeletal growth, polyps can instead allocate more of the available energy pool to tissue biomass [27,28]. For example, in the sea urchin Strongylocentrotus purpuratus, another marine calcifier, protein synthesis accounted for approximately 84% of the available energy pool consumption under exposure to seawater acidification [29]. With this energy-reallocation mechanism, corals could increase their tissue biomass despite having decreased metabolic rates. Similar to our study, a higher-biomass phenotype was observed under a long-term exposure to OA in the temperate coral Oculina patagonica [30]. When transferred back to ambient pH conditions, these corals returned to calcify normally despite being 12 months as soft-bodied polyps in low pH conditions [30], showing that corals possess a remarkable ability to adjust and survive to future acidification conditions. The potential ecological consequences of maintaining overall growth while reducing skeletal development are associated with reaching critical sizes for sexual maturity and assuring a high production of offspring, thereby perpetuating the species and sustaining the reef persistence [7]. The shift to a less-skeleton/more-tissue phenotype would not negatively impact the achievement of sexual maturity and gametogenesis, which appeared to develop similarly in both reduced-skeleton corals at acidified conditions and corals at ambient pH [30]. More porous skeletons with preserved rates of linear extension were observed in S. pistillata adults exposed to acidic conditions [31]. In addition in naturally acidic pH environments, adult corals build more porous skeletons to keep linear extension rate constant [7]. These observations indicate that, for corals, the achievement of larger sizes is indeed a critical need, and that the shift to a more porous skeleton is a common trade-off and an acclimatory strategy to cope with OA. In acidic seawater, our results show a more porous appearance of the primary polyps' skeleton at micro-scale (figure 2l). For corals, the accretion of a more porous skeleton under OA conditions is sustained by embedding more organic matrix proteins within the skeletal fibres [31]. Indeed, our molecular examination shows an enhanced expression of organic matrix proteins at acidic conditions (figure 3c), which points to an increased incorporation of these biomineralization toolkit proteins within the skeletal pores. Previous observations report changes in the expression of organic matrix genes [26] and in the skeletal fibres arrangement [32] of early coral stages reared under OA conditions. Nevertheless, these studies focused either on the molecular or on the morphological aspects, without integrating the coral response across multiple variables. It is noteworthy that energy investment into organic matrix synthesis is suggested to be three orders of magnitude less than the energy cost of active pumps regulating CaCO3 deposition [27]. This implies that increasing the organic matrix production represents a more cost-effective way to support overall polyp growth, when less favourable conditions for skeleton building occur. Despite the change in seawater carbonate chemistry, primary polyps still managed to favour CaCO3 accretion. With our work, we show that coral recruits finely tuned the expression of calcification-related genes. In particular, our molecular analysis implies a tight control of the calcification site chemistry exerted by the coral, through the modulation of Ca2+ transport (figure 3b,c) and of the activity of STPCA-2 (figure 3c). This latter enzyme catalyses the hydration of CO2 to HCO3̄ that is subsequently delivered to the calcification site [17]. The concentrations of both Ca2+ and HCO3̄ have a major role in the skeletal development under OA, as indicated in adults of S. pistillata and other coral species [33,34]. Furthermore, a recent study shows that genes related to Ca2+ management and to inorganic carbon regulation (i.e. carbonic anhydrases) have a major role in the persistence of coral populations adapted to naturally low pH environments [35]. In the light of the predicted decrease of oceanic pH, the ability of corals to modulate the Ca2+ transport and the STPCA-2 activity, especially during the delicate early-life stages, constitutes a crucial area of further research, that could disclose the susceptibility of different coral species to OA. In symbiotic corals, metabolic requirements for growth are critically supported by the photosynthetic activity of the symbiotic algae [9]. Through the carbon fixation process, these dinoflagellates incorporate CO2 into organic compounds, that are used as energy sources [9]. Growing evidence shows that enhanced CO2 levels in seawater could improve the symbiont photosynthetic performance, by alleviating the carbon limitation that is experienced by symbiotic algae inside the host tissues [36,37]. The predicted increase of [CO2] in seawater [10], simulated in our study, appears to boost the activity RuBisCO (the key photosynthetic enzyme that catalyses the first major step of carbon fixation), as suggested by the increased expression of carbon fixation genes in acidic seawater (figure 4). As a result, the algal photosynthetic activity and efficiency increase (electronic supplementary material, figure S4), ultimately enhancing algal growth (electronic supplementary material, figure S3C,E). It is noteworthy that endosymbiotic algae transfer a portion of the photosynthetically derived products, that are not ultimately used for algal growth, to the host [9]. The increased expression of genes linked to carbohydrate binding, that cluster with genes related to the activity of transmembrane transporters and channels (figure 4), suggests an enhanced transfer of photosynthates to the coral host. In corals, heterotrophic feeding requires additional energy to break down food particles and absorb nutrients, compared to sugars coming from the endosymbionts that are quickly metabolized [38]. Given the reduction of metabolic rates under OA, an increased transfer of sugars from the algal endosymbionts might constitute for corals a more cost-effective energy source, compared to heterotrophic sources. Notably, the energy derived from photosynthesis is hypothesized to have a 10–20-fold greater effect on tissue energetics compared to skeletal energetics [27]. Therefore, the augmented algal photosynthetic activity could substantially contribute to the tissue thickening response observed in the primary polyps. In addition, photosynthetic products have been suggested to be used as precursors for skeletal organic matrix biosynthesis [39]. An increased translocation of photosynthates from S. microadriaticum would thus considerably support the enhanced synthesis of organic matrix in the corals, and it would overall increase the host ability to cope with the effects of OA [40]. Earlier works on adult corals investigating the response to OA of the algal endosymbiont focused either on transcriptomic [41] or physiological changes [42,43], and reported negative or no physiological responses in the algal symbionts with the exposure to acidic seawater. However, it must be noted that broad physiological differences exist among different species across Symbiodiniaceae [44], and that gene expression changes of the endosymbiotic algae greatly vary among coral populations in response to acidification [45]. The possibility of symbiont-mediated changes in the coral response to seawater acidification, especially at the highly critical early-life stages, urges, therefore, more in-depth future examinations. In conclusion, we thoroughly describe the mechanisms underlying the response of S. pistillata early-life stages under OA. Our findings suggest that acclimatory mechanisms played a role in the observed coral response. However, we cannot definitively discern relative contributions of acclimation versus natural selection acting on corals genotypes, due to the differential survivorship among treatments and to genetic heterogeneity among coral larvae in the field. In any case, it is important to consider that coral recruits spawned from adults that are already experiencing the effects of decreasing pH may inherit from the parent colonies a greater tolerance to acidifying oceans [46]. These surviving individuals may be better able to cope with OA than prior generations, ultimately building coral reefs resilience through adaptive evolution (reviewed in [47]). As summarized in electronic supplementary material, figure S9, we show that in coral primary polyps, the increased energetic demands of the calcification process under OA resulted in reduced skeletal development and higher skeletal porosity. To counterbalance these changes and support overall growth, corals enhanced the incorporation of organic matric proteins into the skeletal framework and increased the production of tissue biomass. With more available CO2 dissolved in seawater, the growth and photosynthetic activity of the algal endosymbiont was greatly stimulated, leading to a potential enhancement of sugar translocation to the host. While extrapolating laboratory-based findings to projections of the effects of OA on corals in the field is possible, we caution that controlled aquarium systems often do not simulate the dynamic environment of the reef. Moreover, other environmental factors, such as elevated seawater temperature, might act synergistically with acidification, affecting coral fitness [48]. Although extensive, our study has analysed a single coral species and at only one time point, thus limiting the power of our predictions on the acclimatory nature of the observed coral response, and on the phenotypic changes that could occur at later coral life stages under long-term acidification exposure. Thus, future efforts should focus on assessing the degree of modifications in subsequent developmental stages of various coral species under acidification conditions, allowing us to make more accurate predictions on the vulnerability of corals to future OA scenarios. Our multidisciplinary approach reveals strong merit in investigating the response of both the host and the algal partner, as we show that the response to OA involves a wide and intricate net of interrelated parameters in both organisms. The overall correspondence of processes across different assays and biological scales showcases the robust nature of our work, and the importance of including interdisciplinary and complementary analyses towards understanding corals’ vulnerability to environmental change. All data needed to evaluate the conclusions in the paper are present in the paper and/or the electronic supplementary materials and/or are available from the Dryad Digital Repository, including the scripts employed to analyse the RNA-seq data and all the raw data of the study (i.e. raw data of physiological, morphological and photosynthetic measurements, and of the aquariums seawater chemistry: https://doi.org/10.5061/dryad.66t1g1k27 [49]. The RNA-Seq raw data were deposited in the National Center for Biotechnology Information (www.ncbi.nlm.nih.gov/bioproject/PRJNA640748). F.S.: conceptualization, data curation, formal analysis, investigation, methodology, project administration, writing—original draft, writing—review and editing; A.M.: data curation, methodology, software, visualization; P.Z.: data curation, methodology, resources, software, visualization; H.M.P.: conceptualization, data curation, funding acquisition, investigation, methodology, project administration, supervision, validation, writing—original draft, writing—review and editing; T.M.: conceptualization, data curation, funding acquisition, investigation, project administration, resources, supervision, validation, writing—original draft, writing—review and editing. All authors gave final approval for publication and agreed to be held accountable for the work performed therein. The authors declare no competing interests. This project has received funding from the Israeli Binational Science Foundation (BSF 2016321 to H.M.P. and T.M.), the Israel Science Foundation (312/15 to T.M.) and from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 755876 to T.M.). The experiment was performed in a controlled aquarium system which was funded by Institutional ISF grants 2288/16. Computations presented in this work were performed on the Hive computer cluster at the University of Haifa, which is partly funded by ISF grant 2155/15 to T.M. We would like to thank Maayan Neder for the help with collecting coral larvae from the wild. We also thank Jeana Drake for the valuable and critical comments on the manuscript. Lastly, we thank the Interuniversity Institute of Marine Sciences in Eilat for access to its infrastructure and services. FootnotesElectronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5463865. © 2021 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 5The adhesive capabilities of lizards have captivated naturalists since Aristotle, over 2000 years ago (e.g. [1–8]). The lamellae and scansors (sensu [5]) of gecko and Anolis adhesive toe pads are highly specialized. Adhesive ability is facilitated through hair-like, hypertrophied elaborations of the epidermis known as setae [9,10]. At a gross morphological level scansors possess tendinous connections to the digits and either reticular vascular networks or adipose pads to facilitate control while lamellae lack these structures [6,11]. Geckos exhibit a spectrum of digital morphologies, including toe pads with scansors, toe pads with a combination of scansors and lamellae or no adhesive structures at all [5,11,12], whereas Anolis only have toe pads with adhesive lamellae [13]. Excluding an analogous, yet poorly understood, evolutionary origin of digital adhesion in scincid lizards [14], adhesive toe pads are hypothesized to have evolved independently approximately 15 times (approx. 14 gains in gekkotans and one gain in Anolis; [8,12,15]). However, toe pads are not the only adhesive, setae-bearing structures of lizards. Several gecko lineages exhibit setae-bearing, adhesive scansors at the venterodistal tip of the tail. These lineages are geographically and phylogenetically disparate and comprise species in 21 genera in three of the seven gecko families (electronic supplementary material, S1 and S3). While no phylogenetic analyses of adhesive tail evolution have been done, it appears these structures evolved independently at least five times—once each in the families Sphaerodactylidae and Diplodactylidae and three times in the Gekkonidae [16]. The adhesive tail pads of diplodactylid geckos, henceforth called tail pads, are perhaps the most well-studied and are hypothesized to be serial homologues of adhesive toe pads [17]. Serial homologues are morphological structures that are present as multiple copies in the same organism and share a set of developmental constraints, such as fore- and hindlimbs of tetrapods [18–21]. Bauer's [17] hypothesis is based upon striking morphological similarities between adult adhesive toes and tails: reticular networks of blood vessels in addition to muscle fibres attaching directly to the dermal cores of scansors to provide control of the adhesive apparatus, adipose tissue to function as a cushion for the scansors, and of course, fields of setae covering the distal subcaudal tip [17]. Further evidence for the serial homology between tail and toe pads is their apparent evolutionary coupling. The absence of taxa exhibiting tail pads, but no toe pads, suggests that the evolution of toe pads is a prerequisite for evolving tail pads [17,22]. Corroborating the identity of a character as a serial homologue requires developmental data. The only developmental data available to Bauer [17] were a small post-natal series of Rhacodactylus auriculatus, preventing any further corroboration of serial homology. The development of lizard toe pads, in general, is poorly known. In some geckos, the first scansorial ridges form at the distal half of the digit and then develop along the entire length of the digit while becoming more asymmetrical in the proximal–distal direction (Tarentola, Ptyodactylus; [23–26]). In Anolis, the beginning of lamella development follows similar patterns of lepidosaurian scale development [27]; however, the epidermis subsequently undulates, giving rise to asymmetrical lamellae [28,29]. Alternatively, previous studies suggest that all other body scales, with the exception of tail scales, arise from individual, dome-like epidermal papillae (i.e. placodes; [30,31]). This suggests that the developmental programme that gives rise to adhesive toe pads is derived. With the exception of a handful of works [17,22,32–35], the evolutionary morphology of adhesive tail pads has been largely ignored and remains enigmatic and the function and development of these pads have not yet been investigated. The combination of developmental and functional data will provide a robust test of the serial homology hypothesis posited by Bauer [17]. Here we (i) characterize the anatomy and microanatomy of Correlophus ciliatus (Diplodactylidae) tail pads, (ii) characterize the functional capability of an adhesive tail pad in relation to its toe pads and (iii) reinvestigate Bauer's [17] serial homology hypothesis by using developmental data to identify potential developmental constraint in the evolution of adhesive digits and tails. Correlophus ciliatus are large-bodied (108 mm average snout–vent length, SVL), arboreal geckos native to New Caledonia [36–38]. Like all other New Caledonian diplodactylids, C. ciliatus exhibit not only adhesive toe pads, but also an adhesive tail pad at the tip of a robust, prehensile tail (figure 1a). This species was not included in Bauer's [17] investigation because it was thought to be extinct and only ‘re-discovered’ in the mid-1990s [41,42]. To compare tail development, we collected embryos of C. ciliatus and embryos of a digital pad-bearing, but non-adhesive-tailed gecko (Lepidodactylus lugubris) following protocols of Sanger et al. [43] and Griffing et al. [44], respectively. Embryos from both species were collected from captive colonies housed at Marquette University. Protocols for gecko husbandry are detailed elsewhere [37,44]. Using 286 collected embryos, we produced an embryonic staging series for C. ciliatus, the first staging series for any pygopodoid gecko, using published gecko staging series as a reference [45,46]. Figure 1. (a) A time-lapse of a subadult Correlophus ciliatus using its prehensile tail with adhesive tail pad to climb down a branch. (b) Scanning electron micrographs in ventral view of an adult C. ciliatus tail pad. Magnified images of different areas of the tail tip are framed by solid white boxes (distal), dashed white boxes (middle) and dotted white boxes (proximal). Each column is equally magnified relative to the three regions of the tail. Note that the distal tip of the pad exhibits dense fields of setae and more proximal regions have increasingly shorter spinules. (c) Caudal osteology of C. ciliatus depicted through μCT. Ventral view of caudal and sacral vertebrae and portion of pelvic girdle. Dashed white box illustrates autotomic vertebrae magnified in lateral view (1) and transverse view (2). (d) Hall–Brunt Quadruple stained sagittal section of the C. ciliatus adhesive tail pad. Scansor rows 7–12. a, Adipose tissue; d, dermis; e, epidermis; h, hypodermis; hm, hypaxial muscle; s, scansorial unit; sf, setal field. (e) Average heights of a non-pad-bearing gecko (Nactus) spinules, Anolis toe pad setae, Correlophus tail pad setae and Gekko toe pad setae. Figure adapted from Russell [39] and Peattie [40]. (Online version in colour.)
We examined tails and feet of pre- and post-natal specimens using scanning electron microscopy (SEM). Through a separate investigation of embryonic apoptosis, we also identified areas of substantial apoptotic activity within early tail pad development (stage 35 embryo) using Lysotracker Red DND-99 [47]. Setal densities of SEMs were estimated, following Bauer [17]. We investigated the relationship between setal density and maximum seta height in a phylogenetic context combining our C. ciliatus data with data from Bauer [17] and Schleich & Kästle [48] using the phylogeny from Skipwith et al. [49] and phylogenetic independent contrasts [50–53]. We also examined internal anatomy of post-natal tails using computed tomography [54], clearing and staining [55] and histology [56]. We measured frictional adhesive performance from the tails and forelimbs of 10 C. ciliatus over a range of body sizes (6.0–41.3 g) following the methods of Higham et al. [57], in which peak tensile force (Newtons) is obtained by placing the animal's adhesive pad(s) onto a pristine section of acrylic and slowly pulling the animal (or autotomized tail) in parallel opposition to an attached portable force gauge until pad slipping occurs; a single maximum force value was taken after multiple trials for each forelimb and tail. Obtaining reliable measurements often required the tail be autotomized. Once disconnected from the body, the adhesive strength of the tail was measured as above. Adhesive performance in geckos does not require active control [58], and removing the tail first avoided any variation due to behaviour or motivation (i.e. prehension). That said, we obtained tail adhesive performance before and after autotomy for most geckos. We quantified the scaling relationships between adhesive force and body mass using linear regressions. Variables were first log-transformed to linearize the data. The slope of the regression represents the scaling exponent, with a slope of one representing the expected relationship from previous studies of interspecific scaling [7,57]. Scaling relationships were obtained for both the manus and the tail. Additionally, we calculated the potential for the adhesive tail tip to support the entire mass of the animal in a vertical orientation, as might occur when hanging from a branch. To do this we calculated tail adhesive safety factor, the ratio of maximum adhesive force of the tail to the force due to gravity (body mass × acceleration due to gravity). A value greater than one indicates the tail alone could support the body. We investigated scaling of toe pad and tail pad area with relation to SVL using 24 formalin-fixed specimens, ranging from hatchlings to adult (36.8–108.3 mm). After log-transformation, we tested for differences between toe pad and tail pad area versus SVL scaling using linear regressions and analysis of variation [59]. The tail of C. ciliatus comprises 27 amphicoelous caudal vertebrae with reduced transverse vertebral processes, similar to other functionally prehensile-tailed geckos (e.g. Aleuroscalabotes felinus; [60]; figure 1c; electronic supplementary material, S4). Unlike most geckos, no autotomy planes are visible distal to the eighth vertebra. Hypaxial muscle bundles are larger than epaxial muscle bundles (electronic supplementary material, S4). Fields of long, branching setae cover the ventrodistal tip of the adult tail (figure 1b,d). This field occupies the distal most 13–16 scale/scansor rows, with proximal scale rows exhibiting shorter setae, and ultimately, non-branched spinules (figure 1b). Mean height of tail setae = 20.0 µm (N = 30 setae, measured from a ventrolateral tail scansor), the longest measuring 26.6 µm (figure 1b,e; electronic supplementary material, S2). Mean height of toe setae = 23.8 µm (N = 30), the longest measuring 32.7 µm (electronic supplementary material, S5). Setal densities between tail and toe pads are close to each other: tail = 32 950 setae mm−2; and toe = 30 000 setae mm−2. There is an inverse relationship between maximum setal height and setal density in toe pads (r2 = 0.2889, F = 5.468, p = 0.0414), but not tail pads (r2 = 0.2909, F = 3.282, p = 0.1076; electronic supplementary material, S6). Distal scansors/scales are more imbricate and asymmetrical than proximal scales or the scales on the dorsal surface of the tail (figure 1d; electronic supplementary material, S4). Unlike the dorsal scales, the dermal core of scansors sits above a thick hypodermal layer of adipose tissue (figure 1d; electronic supplementary material, S4). The adhesive performance of the manus ranged from 1.33 N to 8.12 N, with larger geckos clinging with greater force (r2 = 0.95, p < 0.001; figure 2a). Adhesive performance of the manus scaled with mass1.04. The adhesive performance of the tail pad ranged from 0.22 N to 1.19 N, and also increased with body size (r2 = 0.88, p < 0.001; figure 2a). Adhesive performance of the tail scaled with mass0.97. The safety factors for the tail pad ranged from 2.45 (largest gecko) to 5.57 (moderately sized gecko). Given an expected scaling exponent of 2 under isometry, the toepad area scaled with isometry (scaled with SVL1.88±0.183), whereas the tail pad area scaled with negative allometry (scaled with SVL1.54±0.298; figure 2b). Figure 2. (a) Adhesive performance of Correlophus ciliatus tail pads and forelimbs. Points to the right of the dashed arrow denotes adults from juveniles/subadults. (b) Isometric scaling of toe pad area and negative allometric scaling of tail pad area with respect to SVL.
Correlophus ciliatus toe pad development begins shortly after interdigital webbing recession (stage 36; figure 3; electronic supplementary material, S7). Four subdigital scansorial ridges initially form in the widest, distal portion of the digit. Shortly after, a small number of new ridges form distally and many more ridges form proximally, all while simultaneously expanding laterally (figure 3). Individual scansors become more imbricate with one another until toe pad development is complete at stage 42 (figure 3). Figure 3. Scanning electron micrographs of digital and caudal development in Correlophus ciliatus and caudal development in Lepidodactylus lugubris. Embryonic stages 36–42. Plantar views of left manus, digit IV and ventral views of the distal tail tip. In C. ciliatus, both toe pads and tail pads exhibit a distal-toproximal development of scansor rows, which subsequently subdivide into what will become the adult scansors. In L. lugubris, scale annuli form synchronously along the entire length of the tail and wrap around the circumference of the structure (stage 38, white arrows). Initial annuli become more distinct (stage 39, black arrow) with annuli appearing near the distal tail tip (stage 39, white arrows). Eventually additional annuli form in between the initial annuli (stage 40, white arrows) and subsequently become granular (stages 41–42). Scale bars = 500 µm.
Tail tip shape changes drastically during C. ciliatus embryonic development. The tail tip is initially pointed and subsequently sculpted away into a wide, blunt end through apoptosis (electronic supplementary material, S8; figure 3). Tail pad development occurs immediately after toe pad development begins. Prior to any signs of scansor development, a subcaudal sulcus forms along nearly the entire length of the tail (stage 36). This sulcus is associated with enlarged hypaxial muscle bundles found in other prehensile-tailed geckos [17]. The distal-most portion of the tail exhibits lateral outgrowths, creating a distal pad which is somewhat wider than the rest of the tail (figure 3; early stage 36). Shortly after (mid stage 36), large lateral scansorial ridges form in the distal portion of the pad before forming new ridges in a distoproximal direction (stage early–mid stage 37). At stage 38, the distal scansorial ridges begin subdividing into numerous, raised units within a scansorial ridge (figure 3). We henceforth refer to this process as granularization. Subsequently (stage 39), the distal portion of the tail expands further laterally, forming a spatulate pad. The pad expands laterally and scales granularize in a distoproximal direction until tail pad development is complete by stage 42 (figure 3). Embryonic development of the non-padded tail of L. lugubris is notably different from C. ciliatus. Lepidodactylus lugubris does not exhibit a subcaudal sulcus and the first signs of tail scale development occur at stage 39 when approximately eight evenly spaced annular scale rows form simultaneously along the length of the tail (figure 3). Further annular scale rows form in between these initial rows by stage 40 and then granularize simultaneously by stage 41 (figure 3). Correlophus ciliatus has been considered unique among adhesive-tailed geckos in having a paddle-shaped distal tail tip, thus expanding the adhesive field compared to tail pads of other species [17,61]. However, despite this unique paddle-shape, the surface morphology and histomorphology of the C. ciliatus adhesive pad is largely similar to other gecko tail pads [17,62,63]. Setae are branched and of comparable height to toe setae, hypaxial muscle bundles are enlarged to presumably assist with prehension, and subdermal adipose tissue likely plays a role in cushioning scansors against surfaces (figure 1b,e; electronic supplementary material, S2 and S4; [5,17,64–66]). Setal heights and densities are largely similar between toe and tail pads (electronic supplementary material, S2). We demonstrate a significant inverse relationship between setal height and density in diplodactylid toe pads. The lack of this relationship in tail pads is likely due to a small sample size. In the gekkonid gecko genus Gekko, Bauer & Good [67] hypothesized that as body size increases between species, setal height and density increase and decrease, respectively. This relationship generally holds true for New Caledonian diplodactylids, but not New Zealand diplodactylids (electronic supplementary material, S2; [17]). Small-bodied Bavayia exhibit moderate setal heights (29–32 µm) with extremely dense fields of setae (35 600–42 900 setae mm−2) while the large-bodied Rhacodactylus exhibit larger setal heights (34–38 µm) with less dense fields of setae (13 700–18 000 setae mm−2). The microanatomy of C. ciliatus tail pads deviates from these trends. Although they exhibit setae of a comparable height to other large-bodied New Caledonian diplodactylids, setal densities of C. ciliatus tail and toe pads are much higher, on par with densities exhibited by Bavayia sauvagii toes and tails (electronic supplementary material, S4 and S2). Further, the setal tip width (i.e. the amount of branching) is larger than any other studied diplodactylids (electronic supplementary material, S2), with both toe and tail setae having a setal tip that is two to three times wider than the closely related Correlophus sarasinorum. The increased density of tail setae, coupled with large setal tip width, may provide C. ciliatus with adhesive ability which exceeds that of other diplodactylid geckos. However, it should be noted that setal densities can vary on gecko toe pads depending on location of measurement [68–70]. Digital adhesion in C. ciliatus is similar to some of the highest absolute adhesive forces recorded for geckos (see [7,57] for comparative values). Gekko gecko is the only species with higher recorded forces [7,58]. Correlophus ciliatus tail adhesion is substantial; forces frequently exceeded 1 N, which, in those cases, represented up to 80.2% of the force estimated to be produced by a single digit. Correlophus ciliatus tail adhesion values far exceed digital adhesion of Anolis, another pad-bearing group of lizards often studied in the context of digital adhesion. With both forelimbs engaged, A. carolinensis and A. sagrei generate 1.5 and 1.3 N of adhesive force, respectively [7]. Thus, crested gecko tail adhesive force often exceeds the forces that these anoles can generate with a single manus. Safety factors for tail adhesion are, in all cases, sufficient to support the mass of the entire body in a vertical orientation, which means that C. ciliatus could potentially hang from a branch using only their tail. In fact, the maximum value of safety factor exceeds five, indicating that a single tail could theoretically hold up to five C. ciliatus without losing grip. Additionally, the tails are capable of prehension, adding yet another component to their clinging ability. Although it is unclear how much grasping force the tail could exert on a perch, some lizards (e.g. chameleons) can exert up to 35 times their own body mass in grip force with their tail [71]. The scaling factors of digital and tail adhesion with respect to body mass in C. ciliatus were not significantly different from 1, which is comparable with previous research on geckos and other pad-bearing lizards [7,57], as well as leaf-cutting ants [72]. This ‘functional similarity’ is not found when examining the relationship between toepad area and adhesive force [7], where it largely follows the predicted scaling factor of 0.67 with respect to body mass. Therefore, other aspects of adhesive morphology, such as setal dimensions or density, are likely driving the functional similarity. Although the scaling exponents of adhesive pad area of the tail and digits (with respect to SVL) are not significantly different from each other, only toepad area scaled with isometry. By contrast, the tail pad scaled with negative allometry, possibly indicating that tails are more important for clinging in smaller geckos. This could be attributed to a shift in habitat use through ontogeny, although little is known about substrate use in nature. Correlophus ciliatus toe and tail pads exhibit strikingly similar patterns of pad subdivision and extension during development (figure 3). Additionally, the onset of pad formation and subdivision occurs near synchronously in both structures. The development of adhesive toe pads in C. ciliatus is similar to Tarentola geckos [24–26] suggesting comparable developmental mechanisms underly the formation of adhesive structures in both species. Consequently, our data provide further support to Bauer's [17] hypothesis that the adhesive tail pads of some diplodactylid geckos are serially homologous to their toe pads. Our results only support this conclusion for the single origin of adhesive tail pads which includes C. ciliatus (clade comprising Pseudothecadactylus + New Caledonian geckos; [49]). A comparative developmental investigation of other groups (e.g. New Zealand geckos, Lygodactylus, Euleptes, etc.) is required to corroborate serial homology in other tail pad-bearing taxa. Further evidence of serial homology between tail and toe pads comes from similarities of setal development on the digits and regenerating tail of the New Zealand diplodactylid, Woodworthia maculata [62]. Such an experiment is impossible with C. ciliatus as they do not regenerate a full tail after autotomy, nor amputation [73]. Although there are underlying differences between toe and tail pads in skeletal, muscular and tendinous morphologies, these likely reflect the distinct ancestry upon which the adhesive pads evolved [17,41] and do not detract from the serial nature of the tail and toe adhesive apparatus. The degree to which tail pads and toe pads are evolving independently (paramorphs; sensu [19]) or in tandem (homomorphs; sensu [19]) is unclear. Although the presence of tail pads appears linked to the presence of adhesive toe pads, anecdotal evidence suggests toe pad shape and size does not predict tail pad shape and size [17,74], suggesting some degree of independent evolution. Further investigations into other tail pad-bearing taxa may determine whether they exhibit different ranges of covariation between tail and toe pad shape and size [75]. Following pad subdivision and extension, tail pads deviate from the developmental pattern seen in toe pads and begin to exhibit granularization of the individual scansors (figure 3). We posit the adhesive scansor developmental programme, via homeosis, was supplanted onto the tail tip, resulting in markedly similar development of scansorial ridges. Soon after, the scansors granularize, creating numerous placode-like structures, which resemble typical reptile body scale development [76]. Unlike C. ciliatus, other adult gecko tail pads, like those of Pseudothecadactylus and Lygodactylus, exhibit mediolaterally broad scansors with few to no granular scansors [17,77]. To our knowledge, this derived, granularizing pattern has not been documented in other studies of amniote integumentary development [76,78]. Identifying the patterns of activator and inhibitor morphogens in tail pad development is required to determine whether this is a two-step process of lateral inhibition to form tail scansors [79]. These derivations further demonstrate that our current understanding of epidermal development is incredibly simplified and requires further descriptive embryology to fully characterize the diversity of epidermal developmental patterns [78]. Further investigation into the molecular patterns and processes that produce both digital and tail adhesive pads is necessary to definitively determine the degree of homology the two structures share. Tail pads of C. ciliatus, but not the non-adhesive-tailed L. lugubris, appear to pass through the three main stages exhibited by developing toe pads (figure 3; electronic supplementary material, S8). Adhesive toe pads of C. ciliatus and Tarentola [24,25] exhibit these stages following digital webbing reduction: (i) pad formation, (ii) distal scansorial/lamellar ridge formation and (iii) distal-to-proximal and lateral ridge extension. By contrast, non-padded lizards of the genus Pogona develop all plantar scales synchronously across the length of the digit [78]. Following these stages, the scansorial rows of developing tail pads begin to granularize, presumably being released from previous constraint (figure 3). Further developmental research is required to determine: (i) if the pattern exhibited by Pogona is the ancestral state, and (ii) if there are biases in the production of morphological variation during pad morphogenesis (i.e. developmental constraint; [80]). Our in-depth investigation into the structure, function and development of C. ciliatus reveals their tail pads are largely similar to other diplodactylid lizards [17,62], with the exception of their extraordinarily dense fields of setae for their body size and large branching setal tips. The adhesive tail pad of C. ciliatus is highly functional, with adhesive capabilities on par with an entire Anolis manus. Paradoxically, the highly functional C. ciliatus tails do not regenerate, unlike nearly all other gecko species which autotomize their tails, including their close relatives, C. sarasinorum. Finally, we add evidence that toe pads and tail pads are serial homologues. Investigation into the molecular underpinnings of toe and tail pad development are required to definitively corroborate serial homology and identify the degree to which developmental constraint has affected the evolution of these enigmatic structures. All Correlophus ciliatus were housed at Marquette University (IACUC protocol AR-279) or University of California, Riverside (IACUC protocol 20170039) following standard husbandry protocols. Morphological and adhesion measurements: Figshare doi:10.6084/m9.figshare.14220098; raw μCT data: Figshare doi:10.6084/m9.figshare.14220134, doi:10.6084/m9.figshare.14222834, doi:10.6084/m9.figshare.14222993; 3D μCT models: Sketchfab https://sketchfab.com/SangerLab. A.H.G.: conceptualization, data curation, formal analysis, investigation, methodology, resources, visualization, writing—original draft, writing—review and editing; T.J.S.: conceptualization, data curation, investigation, methodology, resources, supervision, visualization, writing—original draft, writing—review and editing; L.E.: data curation, formal analysis, visualization, writing—review and editing; A.M.B.: investigation, writing—original draft, writing—review and editing; A.C.: formal analysis, investigation, methodology, writing—original draft, writing—review and editing; T.E.H.: conceptualization, data curation, formal analysis, methodology, visualization, writing—original draft, writing—review and editing; E.N.: formal analysis, investigation, methodology, writing—original draft, writing—review and editing; T.G.: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, resources, supervision, visualization, writing—original draft, writing—review and editing. All authors gave final approval for publication and agreed to be held accountable for the work performed therein. We declare we have no competing interests. This study was funded by the National Science Foundation (DEB 1657662 to T.G.; MRI 1726994 to Jacob Ciszek (Scanning Electron Microscope facility, Loyola University in Chicago)). A.M.B. was supported by the Gerald M. Lemole Endowed Chair Funds through Villanova University. We thank student employees of the Marquette University animal facility; Ryan Kerney for advice about Hall–Bunt Quadruple stain; Tony Russell for comments concerning scansorial identity and translation of [32]. FootnotesElectronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5448688. References
Page 6Although sleep constitutes a period of reduced brain activity and quiescence, specific neural circuitry may be activated at times, and this activity may serve important functions. A well-studied example is neural reactivation resembling daytime patterns in circuitry involved in the formation of motor memory. This so-called replay is thought to aid the consolidation of procedural memory in the acquisition of motor skills [1,2]. Experimental evidence supports this model, as the performance of newly learned tasks is improved by this reactivation [3–7]. Reactivation-related solidification of motor memory has been postulated for a number of different tasks [3,8–11] all of which require the involvement of forebrain and thalamic areas for sensorimotor learning and motor execution. In addition, reactivation of established and stereotyped motor programmes also occurs during sleep [12–15]. This replay can lead to the execution of the motor instructions at the peripheral muscle systems [14,15]. An important example is the activation of the song production circuitry during sleep in oscine songbirds. In the best-studied species, the zebra finch (Taeniopygia guttata), song in adult males is highly stereotyped and consists of 4–8 different syllables, each of which is accompanied by distinct activation patterns of the muscles of the vocal organ, the syrinx [16–19]. Song-like activation during sleep involves contractions of the syringeal muscles without concurrent activation of the respiratory system, thus no sound is generated [15]. The specific muscle activation patterns facilitated the identification of the silent song sequences during sleep and showed that birds only rarely generate the stereotyped syllable sequence of the daytime song. It is not clear whether and how this form of sleep activation is tied to learning, but it may function in the maintenance of stereotyped motor programmes by enabling error processing during sleep [15,20]. Sensory input during sleep can trigger reactivation events [12,13,20–25]. Thus, irrespective of whether sleep activation occurs in the context of newly acquired motor skills or established motor programmes, elaborate and distributed forebrain circuitry is involved in its motor planning, execution and sensorimotor processing. The involvement of this complex circuitry raises the question about the evolutionary origins of these memory-related mechanisms [26]. During sleep, there is a general lack of mobility, which appears to be promoted by ‘sleep neurons’ comprising a network across sensory and motor pathways throughout the brain [27]. These mechanisms for the prevention of movement are not universal, as some motor activity facilitated by skeletal muscles must occur during sleep. Eye movements and spontaneous body movements occur regularly or occasionally during sleep. Another example relevant for birdsong is the activation of syringeal muscles during the expiratory phase of each respiratory cycle in awake and sleeping songbirds [15]. The mechanisms for blocking movement during sleep must, therefore, be circumvented to facilitate activation of peripheral skeletal muscles, whether this activation is linked to life-sustaining activity, spontaneous movement or replay of existing motor programmes [14,15]. Combining these neural observations, motor memory-related reactivation during sleep may constitute neural processes that have evolved from pre-existing mechanisms. Therefore, it is pertinent to test whether song-like motor behaviour occurs during sleep in suboscine birds, which generate song innately [28,29] and do not possess the elaborate forebrain circuitry for acquisition and generation of learned vocal behaviour [30,31]. Experiments were performed on four adult great kiskadees (Pitangus sulphuratus) and on two adult western wood pewees (Contopus sordidulus). Both species belong to the Tyrannidae family. Kiskadees were sexed via polymerase chain reaction amplification of a size-different intron within the highly conserved chromo-helicase-DNA binding protein gene located on the avian sex chromosomes [32]. Three individuals were found to be females, while the sample from the remaining individual could not be amplified for analysis. Pewees were sexed post-mortem by visual inspection of the gonads and determined to be males. Birds were housed in acoustic chambers and subjected to a 14 : 10 h light : dark cycle. In kiskadees, both the male and female produce song. They sing through the year, but most predominantly during the breeding season (October–March in South America). They are most active and vocal in the early morning after sunrise, and in the afternoon, before sunset, although they continue to sing throughout the day. In the case of western wood pewees, only the males sing during the breeding season [33]. Birds were captured using mist nets, they were identified using coloured rings and later transported to the laboratory for experimentation. For kiskadees, no longer than two weeks after capture, birds were released in the same area of the capture and monitored for a few days after release, as described in the permits from the state of Buenos Aires (DI-331-2018-GCDEBA-DFYMAGP). All experiments were conducted according to the regulation of the animal care committees: Institutional Animal Care and Use Committee of the University of Utah (protocol no. 16-03014) and University of Buenos Aires (protocol no. 113, 2019). We measured the electromyographic (EMG) activity of the obliquus ventralis muscle (ovm) a large syringeal muscle [34], which is situated on the ventral surface of the syrinx. To perform this measurement, we used custom-built bipolar electrodes (25 µm diameter, stainless steel 304, heavy polyimide HML insulated, annealed, California Fine Wire). Pairs of electrodes were inserted in the muscle as previously described [35]. For kiskadee recordings, the signal from these electrodes was filtered (150 Hz high-pass RC filter) and differentially amplified (225×) using a custom-built electronic board, which was mounted on a backpack previously fitted on the birds. In wood pewees, electrode signals were amplified (1000×, a larger amplification to ensure signal integrity, since in this case, the amplifier was farther from the measurement site) and band-pass filtered (Brownlee 440; 100–3000 Hz). The sound was recorded using either a condenser microphone (Takstar SGC568) connected to an audio amplifier (Behringer MIC100) or an Audio-Technica (AT 8356) microphone. All conditioned signals were acquired using a National Instruments acquisition board (NIDAQ-USB-6212) controlled using a MATLAB script or using Avisoft Recorder software at 44.15 kHz sample rate. We recorded both signals (sound and EMG) during the day and during the night. EMG signals were then band-pass filtered (100–3000 Hz). Daytime recordings were triggered by a sound level, while night-time recordings were triggered by EMG level (trigger was set at twice the standard deviation of the quiescent signal). In both cases, a 1 s pre-trigger window was recorded. To control that the birds were actually sleeping we performed additional measurements. In pewees, we measured air sac pressure as previously described [15]. In kiskadees, we were not able to perform both measurements simultaneously, so we decided to record the EMGs while inspecting the bird's state using an infrared camera (Sony HDR-SR7) to check that their eyes were closed. During the night, and in addition to EMG, in all cases we recorded sound. EMG signals during song production were recorded throughout the day in two kiskadees. The song of the great kiskadee consists of three syllables, composing a motif. This motif (as in the case of zebra finches) can be repeated a variable number of times. In figure 1a,b, we display an example of a song, in which the three-syllable ‘kis-ka-dee’ is repeated two times. This song structure is stereotyped for each individual, but there is context-dependent variability in the number of repetitions and in the duration and timing of each syllable. Some individuals may also skip the first syllable. The song is accompanied by characteristic patterns of strong EMG activation in the ovm during the first two syllables (figure 1a–c). Figure 1. Song and EMG activity of the kiskadee. (a) Recorded sound of two consecutive songs. (b) Spectrogram of the sound. (c) EMG activity of the ovm muscle. (d) Detail of the EMG activity during the production of the first two syllables of kiskadee song. (e) Frequency spectrum of the shaded segment of the signal in (d), corresponding to the first syllable of the song. (f) Distribution of the mean EMG burst frequencies for the activity during song production in the great kiskadee. The bimodality in the frequency distribution of bursts from daytime vocal behaviour arises from the two types of vocalizations that involve the ovm. Typically, the burst frequency produced during the first syllable of the song was 170 ± 8 Hz, while bursts produced during the isolated ‘kis’ syllables have a frequency of 176 ± 4 Hz. This, together with the fact that the number of events is not balanced between the two vocalizations (n = 25 songs, m = 70 isolated ‘kis’ syllables), explains the asymmetry in the distribution. (Online version in colour.)
The ovm muscle has been shown to be involved in the control and amplification of an amplitude modulation present in the first syllable of the kiskadee song [36], via a burst-like pattern of activation. This muscle also presents a minor activation during the second syllable of the song and is not activated during the third syllable. There are other syringeal muscles, but their activation patterns and roles in sound production are still unknown. Kiskadees use the air sac pressure to control the fundamental frequency and timing of their vocalizations. A detailed view of the EMG signal during the first syllable shows burst-like activation with a characteristic burst frequency (figure 1d,e). This activity is involved in sound amplitude modulation and gives rise to a characteristic spectral enrichment of the sound [37]. The ‘kis’ can also be produced in a call-like manner by itself. To quantify the burst frequency in the EMG accompanying the first syllable, we computed the power spectra and obtained the frequency of its first peak. An example is displayed for the shaded region in figure 1d (corresponding to the first syllable of the song on figure 1a–c), with its spectrum depicted in figure 1e. The burst frequency distribution for 95 vocalizations produced by two birds is shown in figure 1f. This dataset was manually curated to include only vocalization-related patterns of activity observed during sound production, and it consisted of isolated ‘kis’ syllables (n = 75) and song (n = 20) events. During the day, we did not observe vocalization-like patterns of activity without sound production, either during waking hours or during daytime sleep (12 h were measured). In four kiskadees, we recorded EMG activity for at least two nights. We found that two different patterns of strong EMG activation occurred during sleep. During these patterns of activation, no sound was produced. We also checked that the bird's eyes were closed in the video recordings, whenever possible (i.e. when the bird was facing the camera). First, we found song-like activity, consisting of bursts (figure 2a—note the similarity to the EMG activity shown in figure 1d). Figure 2. (a,b) Examples of two replays produced during sleep. (a) Song-like activity, corresponding to the activity of the first syllable of the song (note the similarity with figure 1d). (b) Trill-like activity: long intervals (0.7–1 s) of activity with a slow (10–25 Hz) modulation in amplitude were observed, consistent with the production of a trill. (c,d) Trill produced during daytime. (c) Sound amplitude (blue line) and envelope (orange line) (d) Spectrogram of the sound. (Online version in colour.)
Second, another qualitatively different type of EMG pattern occurred at night, which was typically much longer than a song-like event (about 5×) and is characterized by a slow modulation of the EMG activity (15–20 Hz). Although our birds did not produce a vocalization corresponding to this EMG activity during the day, based on its temporal characteristics and the role of the ovm in sound amplitude modulation, we infer that it is consistent with a ‘trill’ vocalization produced by kiskadees in the context of territorial disputes (figure 2c,d). This vocalization is part of a multi-modal display during which the bird beats its wings and raises its head feathers to display its bright yellow crown (see electronic supplementary material, movie S1). To systematically find all events of song-like activity during the night, we first identified intervals of the relevant activity. This was done by evaluating the envelope of the EMG signal and searching for segments in which it exceeded a threshold. To calculate the EMG envelope, we calculated the signal peaks considering a minimum separation of 0.02 s and then performed a cubic interpolation between the peaks. The considered threshold was 0.04. This threshold was selected for two reasons. First, song EMG activity was found to have peak values an order of magnitude greater than this value. Second, by analysing one hundred 30 s segments of EMG activity, in which the birds were not singing, we were able to estimate the range of resting EMG activity. We found an average standard deviation σ = 0.006 (σ ≤ 0.02 for all cases). We used this low threshold in the first stage to avoid missing any relevant events, and to focus on the subsequent classification (either as noise or song-like activity) of the events. Recently, it was shown that the EMG activity of the ovm consists mainly of bursts of a characteristic frequency, which are produced during the first syllable of the song [37]. The stepwise extraction is exemplified in figure 3 for the data of one of the birds. In order to find song-like replays during sleep, for each segment over the threshold we computed its autocorrelation and the most prominent autocorrelation peak (other than the one at zero-time lag). The inverse of the time lag corresponding to this peak is then the fundamental frequency of the segment. The prominence of the peak, a measure of its harmonicity, was calculated as the height difference between the maxima and its closest minima. In figure 3a, we show a typical song-like pattern of activity. Figure 3b shows its autocorrelation as a function of time lag. From this result, we extracted the most prominent peak (black dot in figure 3b, located near 7 ms) and its prominence (vertical black bar). A full analysis for all the data of one bird is displayed in figure 3c. Each dot corresponds to a segment of above threshold activity, with its size proportional to the prominence of the autocorrelation peak. In this figure, we can see that, while most of the segments aggregate to the left (indicating the absence of a harmonic component), some data points cluster to the right (high autocorrelation), at a specific frequency range (approx. 140 Hz), and with a relatively high value of prominence. The dashed vertical line indicates the minimum autocorrelation value that we considered significant. Figure 3. Detection of activity events in kiskadees. (a) An example of a song-like replay during sleep. (b) Autocorrelation of the signal as a function of lag. The black dot indicates the most prominent peak, and the vertical black bar its prominence. The inverse of the lag of maximal autocorrelation is the burst frequency. (c) Frequency and maximal autocorrelation for all events of one bird. Point size is proportional to the prominence of the autocorrelation peak. Orange dots indicate those manually identified as burst-like activity. The vertical dashed line indicates the autocorrelation threshold (see electronic supplementary material, figure S1) considered for manual classification. (d) An example of a trill-like replay during sleep. (e) Autocorrelation of the signal's envelope as a function of lag. The black dot indicates the most prominent peak, and the vertical black bar, its prominence. The inverse of the lag of maximal autocorrelation is the frequency. (f) Frequency and maximal autocorrelation for all events of one bird. Point size is proportional to the prominence of the autocorrelation peak. Red dots indicate events manually identified as trill-like activity. Note how events identified as trill-like replays are characterized by a relatively higher autocorrelation, a frequency of 14–20 Hz and a high prominence of the autocorrelation maxima (encoded as marker size). We proceeded in a similar manner to select the threshold (dashed vertical line) as in the previous case. (Online version in colour.)
To calculate a threshold autocorrelation, we first computed the average logarithmic prominence of all events above each possible threshold. This curve shows two different regimes (electronic supplementary material, figure S1). For low threshold values, it grows rapidly due to the decreasing number of low-prominence points for increasing thresholds. For high threshold values, it still grows, but at a slower rate. We used this curve to select the autocorrelation threshold as the transition point between the two growth rates. To estimate this point, we computed linear regressions for each regime (see orange and green lines in electronic supplementary material, figure S1) and used their intersection point as the threshold. After establishing a threshold, we manually verified each data point above it to determine its nature. In electronic supplementary material, figure S2, we present three examples of events that lie above the autocorrelation threshold but are qualitatively different. While the three cases presented show similar autocorrelation maxima, prominence and frequency, only the first one was classified as song-like activity. To determine if they were in fact song-like activity, we again used the autocorrelation function. In the case of song-like activity, which has a very defined frequency, the autocorrelation function is smooth and close to zero between the two maxima. In the other cases (both in trill-like activity and in other events, such as EMG activity accompanying defaecation, or movement artefacts caused by the pecking of the surgery suture) the autocorrelation has many maxima and minima, and changes sign multiple times. Neither the duration of the event nor its frequency was used as an aid for this classification. The events that were confirmed as song-like activity are represented in figure 3c by orange dots, while blue dots indicate other types of events (note the absence of song-like events near and above the threshold). To find events of trill-like activity we proceeded in a similar manner, but since we were interested in events with a much slower modulation, we computed the autocorrelation of the signal's envelope. The results of these analyses are displayed in figure 3d–f. In this case, the classification of over the threshold events was aided by the autocorrelation function of the signal's envelope. After event detection, all song-like replays were analysed in the same manner as daytime songs to calculate their frequency and duration. We found that in both quantifications night-time activity differed from song activity. The frequency distributions for song-like activity during sleep and during song production show that sleep replays are characterized by a significantly decreased characteristic frequency (figure 4a, Welch's t-test, p<0.001). While daytime bursts have a frequency of 174 ± 6 Hz, sleep replays have a frequency of 146 ± 4 Hz, representing a 16% decrease in mean frequency. Furthermore, all the sleep replay events found (n = 105 from four birds) were of a lower frequency than any of those produced during daytime song (n = 95; maximum fnight = 155 Hz; minimum fday = 158 Hz). Figure 4. Comparison of awake and asleep song activity. (a) Distribution of EMG frequency for both singing activity (blue) and night-time replays (orange). (b) Duration of EMG bursts produced during awake singing (left, blue), and replayed during sleep (right, orange). Night-time replays show greater variability in duration. (Online version in colour.)
EMG burst duration was more stereotyped and consistent during daytime song production than during night-time replay (figure 4b, Levene test, p<0.001). Daytime bursts are clustered in three groups, which correspond to the specific types of vocalizations included in the dataset. The longest duration bursts stem from two atypical vocalizations shortly after the surgery. The middle cluster was produced exclusively by one of the birds which generated a longer first syllable during the first song of a song bout. During sleep, not only the dispersion of the burst durations increased (from 0.03 s to 0.06 s), but also bursts of shorter duration were produced (down to 0.06 s, while the shortest daytime burst had a duration of 0.09 s). Another difference was found when we analysed the production of the activity corresponding to the second syllable of the song. During the day, the EMG activity was produced either during a call or during the first syllable of the song, in which case it is followed by the activity of the second syllable (after a time gap Δt = 180 ± 40 ms, see the short period of activity after the initial burst in figure 1c) During song production, the motif can be played once, or repeated, in which case the activity of the second syllable is followed (after a gap of Δt = 360 ± 20 ms, n = 10 events) once again by a first syllable burst. During the replays, we found that, whenever the second syllable activity was produced (n = 10), it was followed by a first syllable activity. By quantifying the time difference between these two events we found that it was significantly shorter than the observed time difference during song production (Δt = 150 ± 90 ms), suggesting that during the night-time activation the third syllable of the song may have been omitted, something that we did not observe in any case during daytime song production. We then compared the night-time EMG burst frequency of the presumed trills to the trill rate of sound recordings obtained in the field. We quantified the amplitude modulation frequency of the trill by computing the envelope of the sound and determining modulation frequency as the inverse of the lag of maximum autocorrelation of the envelope. The trill rate (15 ± 4 Hz; n = 16, from two birds) matched the night-time EMG burst frequency very closely (15 ± 4 Hz; n = 93 from four birds). Since the ovm activation displays a dynamics in the same frequency as the modulation of the sound amplitude during the production of a trill, we infer that this EMG activity represents its silent replay. These patterns were also qualitatively different from any other caused by movement artefacts (caused by flight, eating, drinking or pecking of the surgery suture, see electronic supplementary material, figure S2). Interestingly, the recordings of the infrared camera showed that the production of trill-like replays during sleep is consistently accompanied by a partial behavioural display, in which the birds raised their crest feathers to expose their yellow crown (see electronic supplementary material, movies S2 and S3). This behaviour accompanied every trill-like EMG activation event produced during the night. The song of the wood pewee (figure 5a) is accompanied by strong, burst-like activation of the ovm (figure 5c). We did not observe any event of EMG activation without sound production during the day (120 h were analysed). During the four nights of recording in each individual, song-like EMG activity occurred during sleep. To control that birds were sleeping we measured the respiratory rate for sequences of 10 breaths. While the rate during wake hours was (2.12±0.4)Hz (n=180, range 1.74−2.66 Hz), the rate during the night was (1.0±0.1)Hz, (n=380, range 0.85−1.38 Hz). The rate during song-like activation was (1.0±0.1)Hz, (n=60, range 0.82−1.1 Hz). Night-time activation events varied in burst frequency and duration from daytime song. As in the kiskadee, bursting frequency decreased by 15% (Welch's t-test, p<0.001) from song in the daytime (142 ± 3 Hz) to night-time events (120 ± 9 Hz) (figure 5f). As in kiskadees, we also found a significant difference in the duration of the song-like activation. However, while in kiskadees it differed in its dispersion, in this species, bursting activity produced during night-time events (0.07 ± 0.04 s) was significantly shorter (Welch's t-test p<0.001) than the EMG burst produced during song production (0.28 ± 0.02 s) (figure 5g). Figure 5. Analysis of the western wood pewee. (a) Song of the western wood pewee. (b) Spectrogram of the song. (c) EMG activity of the ovm muscle during song production. (d) Detail of the EMG activity during song production. (e) Frequency spectrum of the shaded segment of the signal in (d). (f) Distribution of the EMG frequencies for the song-like activity replayed during sleep and during singing for the western wood pewees. Burst frequency is significantly lower during sleep replays. (g) Duration of EMG bursts produced during awake singing (left, blue), and replayed during sleep (right, orange). Night-time replays are significantly shorter. (Online version in colour.)
The song-like activation patterns of syringeal muscles during sleep in the two suboscine species extend our records of motor reactivation to motor behaviour, which (i) is innate (i.e. not learned), (ii) might occur in males and females, (iii) includes vocalizations that are not considered song and (iv) is executed by circuitry residing in the mid- and hindbrain, rather than the forebrain [31]. Suboscines in the family Tyrannidae develop song innately without the need for acoustic models or auditory feedback [28,29]. Unlike oscine songbirds, they do not exhibit the clearly delineated forebrain and thalamic circuitry necessary for learning from tutor models and for generating learned vocal behaviour [30], although some connections from the arcopallium to mid- and hindbrain vocal and respiratory areas appear to exist [38,39]. It is thought that the vocal motor programme arises in the mesencephalic dorsomedial nucleus of the intercollicularis (DM, equivalent to the mammalian site for vocal control, the periaqueductal grey, PAG). Electrical stimulation in this region can elicit unlearned, species-specific vocalizations [40–42]. Thus, unlike in oscines, the execution of the song motor programme does not require contributions from the forebrain. Telencephalic and hypothalamic input to the DM must play a role in initiating song production [31], but it is not required for the control of normal, species-typical vocalizations. The reproduction of song-like activation of the syringeal muscles in both suboscine species is similar to that seen in the zebra finch in several regards [15]. Its occurrence is distributed throughout the night, and its timing in the zebra finch suggests that it occurs during slow-wave episodes in NREM sleep [43]. This is consistent with the finding that auditory responses in the sensorimotor nucleus HVC preferentially occur during this sleep state [44]. A further parallel is lower stereotypy in the reactivation patterns, which results in different durations of the burst trains and different burst frequencies. In oscine birds, sleep has been shown to be involved in the development of learned vocal behaviour [8,45]. In species of this group, sleep consists of both rapid-eye-movement (REM) and non-REM sleep [46,47] states that are comparable to those observed in mammals [48]. These sleep patterns have been found in many species [49], but to the best of our knowledge, specific data for suboscines are not available. Whereas our indirect evidence (IR videos and respiratory rate) indicates that birds were asleep during the song-like activation events, we could not determine whether these events occurred during REM or non-REM phases. In addition, the kiskadee replay of the motor programme of another vocalization that is not considered to be song is remarkable. In the case of the kiskadee trill, even part of the visual component of the multi-modal display is replayed during sleep. This adds a motor programme for a complex multi-modal display to the known cases of replay events. As with the suppression of the respiratory contribution in the vocal display component to prevent sound production, the wing movements of the visual display are not executed, while head feathers are raised to display the yellow crest. In both cases, display components that might alert potential predators to the sleeping bird are suppressed during the replay of displays. The burst repetition rate during sleep tends to be slower during song in kiskadees and wood pewees, but not during the trill in kiskadees. The latter observation makes it less likely that the slower rhythm is caused by lower brain temperature during sleep [50,51]. Alternatively, or in addition, intrinsic sleep-induced changes in circuit dynamics may differentially affect the rhythms. Melatonin may play a critical role in these timing changes. The slowing of the rhythm and the shorter duration of song-like activation patterns compared to daytime song agree well with observations following reduced melatonin production in the vocal learning zebra finch and the non-learning Japanese quail (Coturnix japonica). Song and crowing were shorter and had a faster rhythm when night-time melatonin production was prevented either through exposure to constant light or pinealectomy [52]. In birds with innate vocal behaviour, melatonin receptors are present in the intercollicular vocal nucleus [53], suggesting that melatonin could directly affect the rhythm of vocalizations where the vocal motor execution presumably originates. Although the songs of the two species in this study are comparatively simple, the changed timing between second syllables and the following first syllable of the next song suggests that deviations from stereotyped song sequence also occur, as was observed in the zebra finch [15]. The parallels to song-like reactivation in the oscine zebra finch indicate that the observed activity in suboscines also constitutes ‘replay’ of song and even multi-modal displays. The fact that suboscines show these remarkable similarities warrants questions about the evolution of motor replay during sleep. It is currently unclear whether tyrannid suboscines secondarily lost the vocal learning ability or whether vocal learning may have arisen independently multiple times in the clade including parrots and passerines [38,54]. In either case, the fact that the midbrain centre for unlearned vocal behaviour is activated during sleep suggests an ancestral mechanism that may have predated vocal learning. The possible functions of this reactivation in innate vocal behaviour remain to be investigated. Among the exciting possibilities are a neutral hypothesis, peripheral motor system maintenance and/or general consolidation of motor programme circuitry [15]. The latter mechanism fits nicely with the variable reactivation patterns of the stereotyped daytime song in oscines [15] and may, therefore, be a precursor to the more sophisticated telencephalic sensorimotor processes for learned vocal behaviour. All experiments were conducted according to the regulation of the animal care committees: Institutional Animal Care and Use Committee of the University of Utah (protocol no. 16-03014) and University of Buenos Aires (protocol no. 113, 2019). Kiskadees were released in the same area of capture and monitored for a few days after release, as described in the permits from the state of Buenos Aires (DI-331-2018-GCDEBA-DFYMAGP). Code and processed data available from the repository http://doi.org/10.5281/zenodo.4599099 [55]. We declare we have no competing interests. This work was supported by ANCYT, CONICET and UBA. We thank Alejandro Martínez for providing the movie of a great kiskadee performing a trill in the wild, and Melina Atencio (UBA) for her help sexing the kiskadees in our study. Two anonymous referees provided helpful suggestions, which improved the manuscript. FootnotesElectronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5470465. References
Page 7One of the most extreme ecological transitions has been the colonization of land by fishes. Water and air differ dramatically in their physical properties (e.g. density and viscosity), which pose critical challenges for processes such as respiration, osmo- and ionoregulation, nitrogen excretion, feeding and locomotion [1,2]. These challenges suggest that traversing the air–water interface should be difficult for fishes and consequently quite rare, yet the ability to leave water has evolved approximately 87 times among bony fishes in the past 65 million years [3]. There are more than 200 extant species of amphibious fish that transition between water and land as part of their natural history [4]. While there are many notable challenges associated with an amphibious lifestyle, spatial navigation is one difficulty that has been largely overlooked. Some amphibious fishes must learn to navigate one habitat in both the presence and absence of water, such as mudskippers and blennies that become air-exposed due to receding tides [5,6]. Others must learn to navigate dramatically different aquatic and terrestrial habitats, such as killifishes that move overland from pond to pond [7,8]. Fish that traverse the air–water interface arguably experience a more complex suite of environmental conditions than those that remain solely in water. Terrestrial environments provide a number of novel biotic (e.g. predators) and abiotic (e.g. wind and rain) stimuli that are not necessarily encountered in aquatic environments. Moreover, fish on land must negotiate obstacles and terrains that could otherwise be avoided at appreciable water depths. Previous studies have described how amphibious fishes may use vision [9–12], chemoreception [13] and the otolith–vestibular system [14,15] to navigate unfamiliar aquatic and/or terrestrial environments. However, learning the location of important environmental features (e.g. predators, food, refuge and mates) in both aquatic and terrestrial settings can be critical for survival and reproduction [16]. There is mounting evidence that spatial cognition is a highly plastic trait in fishes. Exposure to dynamic and structurally complex environments can improve spatial learning [17,18], and this effect is strongly associated with neural plasticity in the dorsolateral telencephalic pallium (i.e. homologue to the mammalian hippocampus [19]). For example, juvenile Atlantic salmon (Salmo salar) raised in tanks with a rocky substrate and aquatic plants exhibited increased neurogenesis in the telencephalon, and were faster at finding their way out of a four-armed maze than those from barren tanks [20]. Physical exercise can similarly have positive effects on neurogenesis and spatial learning, although the vast majority of studies have focused on mammalian models [21]. Over the past few years, however, studies have emerged, suggesting that swimming can have similar effects in fishes [22,23]. Locomotor movement in air is far more difficult than movement in water [24]. Even small bouts of terrestrial movement may therefore constitute exercise for amphibious fishes. Indeed, just 3 min of terrestrial exercise on alternate days induced dramatic skeletal muscle remodelling in the amphibious mangrove rivulus (Kryptolebias marmoratus) after only two weeks [25]. If amphibious fishes face more complex environmental conditions and more challenging locomotion when they leave water, then does leaving water make fish smarter? We used an amphibious killifish (K. marmoratus) to test the hypothesis that the spatial learning ability of amphibious fishes would be altered by exposure to terrestrial environments because of altered neurogenesis in their dorsolateral pallium. We predicted that fish subjected to air–water fluctuations and terrestrial exercise would exhibit enhanced neurogenesis in their dorsolateral pallium, as well as better maze learning abilities, than fish maintained solely in water. Kryptolebias marmoratus inhabits ephemeral ponds and crab burrows throughout mangrove swamps of the tropical western Atlantic [26] and frequently leaves water to escape unfavourable aquatic conditions (e.g. hypoxia [27]), and also to disperse, forage or reproduce [26,28]. During seasonal droughts, K. marmoratus may also be forced out of water for several weeks [26,29]. Once out of water, K. marmoratus encounters complex terrestrial terrains (electronic supplementary material, figure S1) which they move across using a terrestrial ‘tail-flip’ jumping behaviour [30]. Although previous studies have investigated the sensory cues that allow K. marmoratus to navigate unfamiliar environments [15], whether fish learn the location of ecologically important features in the environment remains ambiguous. We subjected fish to either eight weeks of fluctuating air–water conditions or terrestrial jumping exercise before assessing spatial learning and neurogenesis in the dorsolateral pallium. All experimental fish (n = 85, 0.094 ± 0.001 g, 19.0 ± 0.1 mm; mean ± s.e.m.) were adult hermaphrodites of the self-fertilizing K. marmoratus (strain HON9, originating from the Bay Islands, Utila, Honduras [31]). Prior to experiments, fish were individually maintained in 120 ml plastic holding cups (approx. 60 ml water, 15‰ salinity, 25°C) in the Hagen Aqualab at the University of Guelph on a 12 light : 12 dark cycle. Fish were fed live Artemia sp. nauplii three times weekly. All experimental procedures were approved by the University of Guelph Animal Care Committee (AUP 3891). Fish were randomly assigned to one of three 8-week experimental acclimations: control (n = 29), air–water (n = 28) and terrestrial exercise (n = 28) (electronic supplementary material, figure 2). At the start of the acclimation period, we transferred each fish from their holding cup into a 750 ml container (Rubbermaid TakeAlong Deep Squares; 140 × 140 × 80 mm) with a rugose paraffin wax bottom (electronic supplementary material, figure 3). The purpose of the rugose wax bottom was to mimic the complex terrestrial terrain that K. marmoratus must navigate when air-exposed under natural conditions. The air–water containers were drained and refilled with the same water every 1–3 days (randomly assigned) such that water was absent for half of the acclimation period, except within the deepest crevasses of the paraffin wax bottom. We maintained constant water levels in the control and exercise containers (approx. 500 ml water; 15‰ salinity, 25°C) throughout the acclimation period. On the same days that water was absent from the air–water containers (3–4 d week−1), fish from the exercise treatment were removed and placed on moist filter paper in a terrarium (30 × 60 cm). After a 2 min adjustment period, fish were exercised as previously described [18]. Briefly, fish were encouraged to jump repeatedly via gentle prodding with a clicker ballpoint pen for 3 min (approx. 50% exhaustion [25]). Fish in the air–water treatment were fed Artemia sp. nauplii each day they were in water (3–4 d week−1), and the control and exercise groups were fed on the exact same days. All fish were also fed bloodworms once per week on the same day. Water changes were performed in all groups every second week. At the end of the eight-week acclimation, spatial learning was assessed in a subset of fish from each treatment group (control, n = 20; air–water, n = 19; exercise, n = 20) using a bifurcating T-maze modified from the existing designs [32] (figure 1). The T-maze was made of opaque grey plastic; a down piece (start arm; 24 × 6 × 6 cm) was connected to the midpoint of a cross piece (30 × 6 × 6 cm) to create the two goal arms. We familiarized fish to the T-maze by allowing them to swim freely within the maze in the 3 days prior to the learning trial (1 h d−1). To assess spatial learning, we suspended a bloodworm in each goal arm, which was blocked from view using two staggered opaque partitions (figure 1). The bloodworm in one arm was accessible to the fish (reward), while the bloodworm in the opposite arm was made inaccessible by placing it inside a mesh bag. The location of the reward bloodworm (i.e. right or left goal arm) was randomly assigned for each fish. At the start of each trial, an individual fish was collected in a hand-held dip-net and placed behind a removable opaque partition approximately 4 cm from the end of the start arm. After a 5 min adjustment period, the partition was lifted and fish were given 30 min to find (and eat) the reward bloodworm before they were returned to their acclimation container. Preliminary trials revealed that fish routinely found the reward bloodworm within a 30 min period. All fish were fasted for 3 days prior to the learning trials to ensure a similar motivation to find the reward bloodworm. We tested fish once per day over 10 consecutive days between 08.00 and 12.00. Each trial was video recorded using a webcam (Logitech Quickcam Pro, Fremont, CA, USA) mounted above the maze, so that fish performance was not affected by an observer. Videos were analysed to determine (i) the latency to find the reward bloodworm, (ii) the distance fish travelled before finding the reward bloodworm (standardized to body length), and (iii) the number of times fish entered a non-reward arm before finding the reward bloodworm, hereafter referred to as the number of errors. We used ImageJ (http://imagej.nih.gov/ij) to calculate the distance travelled by marking the location of the fish in the maze each second until the reward bloodworm was eaten and then summing the distance between each point. Finally, we also measured the time fish spent attempting to eat the inaccessible bloodworm on the first day of the learning trial. Figure 1. Schematic of the T-maze used for spatial learning trials.
We assessed cell proliferation as a proxy for neurogenesis in the dorsolateral pallium of each fish using the S-phase marker, proliferating cell nuclear antigen (PCNA). The dorsolateral pallium is thought to be homologous with the hippocampus of mammals and birds—the brain region responsible for spatial learning [33,34]. Following the eight-week acclimation period, the subset of fish that were not assigned to the learning trials (control, n = 9; air–water, n = 9; exercise, n = 8) were euthanized via cold-water immersion and immediately decapitated. Whole heads were fixed in 10% buffered formalin for 24 h and stored in 70% ethanol (4°C) until routine paraffin embedding was performed. The paraffin-embedded heads were serially cross-sectioned at 5 µm increments through the telencephalon (approx. 200 µm range) and then stained for PCNA [35]. Briefly, paraffin sections were deparaffinized in xylene and rehydrated with a graded ethanol series. Antigen retrieval was then performed by submerging the sections in citrate buffer (10 mM sodium citrate, 0.05% v/v Tween20; pH 6) at 95°C for 12 min and cooled at room temperature for 20 min. Following antigen retrieval, the sections were blocked (Immobilon Block, Sigma Millipore) for 1 h at room temperature, then incubated overnight at 4°C in primary antibody (1:200 PCNA; Proteintech Group) and diluted in phosphate-buffered saline (PBS; 137 mM NaCl, 15.2 mM Na2HPO4, 2.7 mM KCl, 1.5 mM KH2PO4; pH 7.8). The sections were then rinsed in PBS (3 × 5 min), incubated for 2 h at room temperature with secondary antibody diluted in PBS (1:500 Alexa Fluor 488 goat anti-rabbit IgG; Invitrogen), counterstained using DAPI and coverslipped. We photographed three sections per individual in which the dorsolateral pallium was intact using an Eclipse Ti2 Series inverted epifluorescent microscope (Nikon Instruments, Melville, NY, USA) and averaged the total number of PCNA + cells in the dorsolateral pallium (right and left side of the brain combined) across the three replicates. We used one-way analyses of variance (ANOVAs), followed by Tukey's post hoc tests, to compare baseline (day 1) maze performance (i.e. latency, distance and errors) between groups. We performed a linear regression to determine the relationship between the latency to find the reward bloodworm and the time fish spend attempting to eat the inaccessible bloodworm on day 1. To assess spatial learning, we first calculated the cumulative time, the cumulative distance travelled and the cumulative number of errors made for the successive experimental days, correcting for baseline performance. The resultant ‘learning curves’ were fit with a second-order polynomial model to allow for a curved relationship. We then tested whether the slopes of the ‘learning curves’ differed between the experimental treatments. Finally, we used a one-way ANOVA to compare the number of PCNA + cells in the lateral pallium between groups. All data were assessed for normality of residuals (Shapiro–Wilk) and homogeneity of variance (Bartlett's test), and appropriately transformed when necessary. All statistical analyses were performed using RStudio (v. 1.1.447), and all graphs were generated using GraphPad Prism (v. 8). The experimental treatments influenced baseline maze performance in K. marmoratus. The latency to find the reward bloodworm differed between groups on day 1 (ANOVA; p = 0.01; figure 2a). Control fish took significantly longer to find the reward bloodworm than both the air–water (Tukey's; p = 0.02) and the exercise group (Tukey's; p = 0.01). We found a strong positive correlation between the latency to find the reward bloodworm and the time fish spend attempting to eat the inaccessible bloodworm on day 1 (linear regression; p < 0.001, R2 = 0.43; electronic supplementary material, figure S4). We found no significant difference in the distance travelled (ANOVA; p = 0.20) or the number of errors made (ANOVA; p = 0.24) between groups on day 1. The experimental treatments also influenced spatial learning in K. marmoratus. Compared with control fish, the air–water and exercise groups travelled an increasingly shorter distance before finding the reward bloodworm as the trial progressed (curve comparison; p = 0.01; figure 2b). Although we found no such statistical differences in the time it took to find the reward bloodworm (curve comparison; p = 0.14; figure 2c), fish from air–water and exercise groups spent 10% and 26% less cumulative time, respectively, searching for the reward compared with control fish. The experimental treatment had no effect on the number of errors made before finding the reward bloodworm (curve comparison; p = 0.60). Finally, we found that the experimental treatments influenced cell proliferation in the dorsolateral pallium of K. marmoratus (ANOVA; p = 0.026; figure 3). The exercise group had significantly more PCNA + cells in their dorsolateral pallium than the control group (Tukey's; p = 0.03). The number of PCNA+ cells in the dorsolateral pallium of air–water fish did not significantly differ from that of the control (Tukey's; p = 0.11) or exercised fish (Tukey's; p = 0.75). Figure 2. (a) Latency to eat the reward bloodworm on day 1 of the learning trial. Different uppercase letters denote significant differences between the experimental treatments. (b) The cumulative distance travelled (in body lengths), and (c) the cumulative time elapsed (in seconds) before eating the reward bloodworm. The lines indicate the second-order polynomial fit for each treatment group. An asterisks indicates that the ‘learning curves’ differed significantly across the experimental treatments. Error bars represent s.e.m. n = 19–20 per group.
Figure 3. (a) A representative transverse section through the telencephalon of K. marmoratus. The left half is a schematic that indicates the different telencephalic regions in the fish brain, and the right half is a representative photo of a stained paraffin section. The red box indicated the location of (b). Abbreviations: dD, dorsodorsal; DMd, dorsomedial dorsalis; DMv, dorsomedial ventralis; DLd, dorsolateral dorsalis; DLv, dorsolateral ventralis; DC, doralis centralis; DP, dorsalis posterior; VT, ventral telencephalon. (b) Representative region of the dorsolateral pallium indicating PCNA+ cells in green and DAPI-stained nuclei in blue. Scale bar = 50 µm. (c) The number of PCNA+ cells in the dorsolateral telencephalic pallium (DLv) of each treatment group. Different uppercase letters denote significant differences between the experimental treatments. Error bars represent s.e.m. n = 8–9 per group. Means are shown as crosses within the box.
We hypothesized that neurogenesis and the spatial learning ability of amphibious fishes are plastic and influenced by terrestrial exposure. Indeed, we found that K. marmoratus in the air–water and terrestrial exercise groups were taking a more direct route to the food reward compared with control fish towards the end of the learning trials, suggesting that they had better spatial learning abilities. Moreover, air–water and exercised fish had 39% and 46% more proliferating cells in their dorsolateral pallium relative to control fish, respectively, which likely reflect higher rates of neurogenesis [36,37]. We also found that baseline maze performance (day 1) differed between treatment groups. Air–water and exercised fish found the food reward significantly faster than control fish, in part because they spent less time trying to eat the inaccessible bloodworm. Taken together, our findings suggest that fish that spend more time on land may have a cognitive advantage over those that remain in water, and in turn may be more successful at navigating both aquatic and terrestrial environments with potential fitness advantages. The eight weeks of experimental acclimation dramatically altered baseline maze performance in K. marmoratus. The air–water and exercised fish, which found the reward bloodworm significantly faster on day 1, were possibly employing a ‘lose-shift foraging strategy’, whereby fish quickly moved on in the maze after realizing that the non-reward bloodworm was inaccessible [38,39]. Alternatively, enriched environments and physical exercise have both been shown to reduce stress in fishes and result in bolder behaviours, such as being more exploratory [40,41]. If fish subjected to air–water fluctuations and exercise were bolder and more exploratory than control fish, it may explain why they found the reward bloodworm faster on day 1 of the learning trial. Different baseline maze performances between groups may also relate to alterations in the ‘latent’ spatial learning ability of K. marmoratus. Latent spatial learning involves the gradual creation of a cognitive map of the environment in the absence of a reinforcement (e.g. food reward) [42]. We allowed fish from all groups to become familiar with the maze in the 3 days prior to the learning trial in the absence of bloodworms. If air–water and exercised fish had superior latent learning abilities, or simply explored the maze more thoroughly than controls during the familiarization period, it might explain their better maze performance on day 1. Other amphibious fishes are known to build spatial maps of their environment, such as the frillfin goby (Bathygobius soporator) that memorizes the topographical features of the environment when submerged at high tide in order to efficiently move from pool to pool during low tide [43]. Although there are a number of possible reasons why baseline maze performance may have differed between groups, fish that were exposed to the terrestrial environment were faster foragers than those maintained solely in water. Exposure to fluctuating air–water conditions significantly improved some markers of spatial learning in K. marmoratus. Notably, air–water fish travelled an increasingly shorter distance as the learning trials progressed relative to controls, suggesting that their route to the food reward was becoming more efficient. We also found a trend towards increasingly faster maze performance in the air–water group, but the number of errors fish made in the maze did not differ from the control. Fish from the air–water group experienced the structurally complex terrain of their acclimation container under both aquatic and terrestrial conditions, and were forced to physically interact with the terrain when air-exposed. Numerous studies have demonstrated that acclimation to structurally complex environments can have positive effects on spatial cognition and neurogenesis in vertebrates, likely because manoeuvring around obstacles and/or over tortuous terrains pose an increased demand on cognitive functions related to spatial navigation [18,20,44]. Air exposure may have acted as a form of environmental enrichment by providing novel stimuli (e.g. puddles in wax crevasses) that are not necessarily encountered under aquatic settings and by forcing fish to negotiate obstacles that could have been avoided at appreciable water depths (e.g. wax ridges). On the other hand, exposure to terrestrial environments involves a number of morphological and physiological changes that could have also impacted brain function. For example, air-exposed K. marmoratus proliferates epidermal capillaries in the skin to enhance their capacity for aerial O2 uptake [45] and increase their blood O2-carrying capacity [46]. Increased circulating blood O2 has been linked with improved cognitive performance in mammals [47], but it is unclear whether similar mechanisms can improve cognition in fishes. Finally, it is important to note that fish remained relatively inactive during periods of air exposure compared with periods in water, although we observed them periodically moving overland within their acclimation containers (G.S.R. 2020, personal observation). Terrestrial movement may therefore have contributed to the improved spatial learning of fish from the air–water group. Regardless, it is clear that air–water fluctuations in general can have dramatic effects on neural and cognitive processes in K. marmoratus. Exercise training is well known to enhance neurogenesis and spatial cognition in vertebrates [48], although relatively few studies have focused on fishes. Remarkably, we exercised fish for only 84 min out of the 80 640 total minutes (0.1%) during the eight-week experimental period, yet this was sufficient to induce significant neural (increased PCNA + cells) and cognitive (better maze performance) changes in this species. Since we were unable to exercise fish in water (K. marmoratus refuse to swim against a current), it is difficult to disentangle the effects of exercise from air exposure. Although brief, a few minutes of air exposure three to four times per week may have resulted in some of the neural and behavioural changes observed in our study. It is also possible that exercise may induce neural and cognitive plasticity in K. marmoratus regardless of whether it occurs in water or on land. Swim training was found to increase the expression of neurogenesis-associated genes in the dorsolateral pallium of Atlantic salmon (S. salar), including the expression of pcna [23]. However, wild K. marmoratus have been reported to spend up to 90% of their time out of water [49], and thus, terrestrial locomotion is probably a more ecologically relevant form of exercise in this species. Finally, it is unknown how the positive effects of exercise on neurogenesis are mediated in fishes, but in mammals it is thought that hippocampal angiogenesis plays an important role [50,51]. Perhaps similar mechanisms are responsible for enhancing neurogenesis in the dorsolateral pallium of K. marmoratus as in the mammalian hippocampus, although this hypothesis warrants further investigation. In the wild, K. marmoratus are found in ephemeral pools in which they experience air–water fluctuations, but they also leave water voluntarily and traverse terrestrial landscapes in search of prey and/or new aquatic environments [26,28]. Although our study was performed in captive fish, our findings indicate that terrestrial episodes promote neurogenesis in the telencephalon and enhance the spatial learning ability of this species. Thus, spending more time out of water on a daily basis during the wet season may improve navigation. However, when seasonal droughts occur, K. marmoratus seeks moist terrestrial habitats (e.g. rotting logs) and may remain quiescent for weeks at a time in complete darkness [52]. The lack of movement and the deprivation of environmental stimuli may reverse the terrestrial enhancement on neurogenesis and spatial learning, an avenue for future investigation. Finally, the capacity for spatial learning and its neural mechanisms are highly conserved across vertebrate groups, suggesting inheritance from ancient fishes that gave rise to land-dwelling tetrapods [19,33,53]. Is it possible that ancient fishes similarly exhibited neural and behavioural plasticity from repeatedly traversing the air–water interface? If so, could it have helped them navigate and successfully colonize novel terrestrial environments? MacIver et al. [54] recently demonstrated that the vertebrate invasion of land was preceded by the evolution of an increased visual range above water, which is thought to have promoted an increased capacity for spatial planning in ancient fishes. Since the colonization of land by extant amphibious fishes has obvious parallels with the origin of all land vertebrates, understanding the factors that facilitate successful land invasions in extant fishes can provide important insights into the evolutionary path to terrestrial life. All experimental procedures were approved by the University of Guelph Animal Care Committee (AUP 3891). The datasets supporting this article have been uploaded as part of the electronic supplementary material [55]. G.S.R.: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, validation, visualization, writing—original draft, writing—review and editing; P.A.W.: conceptualization, funding acquisition, supervision, validation, writing—review and editing. All authors gave final approval for publication and agreed to be held accountable for the work performed therein. We declare we have no competing interests. This work was provided by the Natural Sciences and Engineering Research Council of Canada (NSERC) Canada Graduate Scholarship to G.S.R. and an NSERC Discovery grant to P.A.W. We thank Mike Davies and Matt Cornish for assistance with animal care, Dr Andreas Heyland for the use of his microscope and Dr Fred Laberge for helpful discussions regarding data interpretation. FootnotesElectronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5448685. References
Page 8The biological pump is the process by which organic nutrients are transported from shallow ocean to deep sea [1]. Today this consists of two major vectors: passive sinking of organic matter aggregates and vertical movement of animals [1] (figure 1). The biological pump, the main driver of the marine carbon cycle, is responsible for approximately two-thirds of the vertical gradient of carbon in the ocean [2]. While transport of carbon to the deep ocean through vertical mixing of dissolved organic carbon (DOC) carries significant amounts of carbon to depth, 95% of DOC cannot be used as food by marine organisms [1,3]. By contrast, aggregates and vertical migrants concentrate organic matter in a form that can be used by demersal animals [1] and are therefore crucial for the establishment and sustenance of deep-water communities. Figure 1. Simplified architecture of the modern biological pump. Arrows trace pathway of carbon (C) and nitrogen (N) and other nutrients. Numbers indicate two main biologically mediated vectors that transport nutrients fixed by phytoplankton from the surface ocean to demersal communities. (1) Aggregation vector of phytoplankton, faecal pellets and other organic matter which sinks passively through the water column. (2) Vertical migration vector driven by active two-way migration by metazoans. (Online version in colour.)
The biological pump was very different before the Cryogenian Period (more than 720 Ma). Cyanobacterial picoplankton (0.2–2.0 µm) domination led to a stratified and turbid water column (e.g. [4–6]). Cells were too small to sink in significant numbers, so most nutrients were recycled within the surface waters, with little export to the deeper ocean [4–6]. The flux of the aggregation vector (figure 1) would have increased when primary productivity shifted to be eukaryote-dominated during the Cryogenian ∼ 650 Ma [7], and further when Cambrian suspension feeding sponges and pelagic phytoplanktivores applied a size-selective pressure for larger primary producers [8–12]. The carcasses, sloppy feeding and faecal pellets of macrozooplankton, including phytoplanktivorous crown-group crustaceans (e.g. [9]), would have further increased the formation and size of aggregates, and thus the flux along this vector. However, despite the importance of vertical migration in the transport of nutrients necessary to sustain mesopelagic and deep-water ecosystems in the modern ocean (e.g. [13]), this vector is poorly understood in the Cambrian. Comparisons of the morphology of Cambrian fossil organisms with modern vertically mobile pelagic animals provides the opportunity to infer whether vertical migration occurred in oceans over 500 Ma. Qualitative comparisons between the pelagic crustacean Gnathophausia and the nektonic stem group euarthropod [14] Isoxys (e.g. [15,16]) suggest the latter is a promising candidate. The 20 Isoxys species so far formally described are united by the presence of a bivalved carapace which bears both anterior and posterior spines [17] (figure 2), and the genus has a cosmopolitan distribution (e.g. [18]). The presence of eyes that can comprise approximately 10% of the body length, a digestive tract with paired serial midgut glands, and a pair of anteriorly positioned raptorial appendages (figure 2) support a predatory habit for Isoxys, which would have been well suited for capturing small soft-bodied invertebrates [19]. Isoxys is unusual for Cambrian animals as a pelagic lifestyle has been proposed (e.g. [18]), although some recent workers have suggested a potential hyperbenthic lifestyle (1–10 m above the bottom), with individuals capable of moving small distances vertically [19,20]. However, while Isoxys carapaces appear to show adaptations for hydrodynamic streamlining, interspecific differences in both carapace asymmetry and spine lengths (e.g. figure 2), as well as soft parts, suggest that different species may have occupied distinct niches, including some much closer to the seafloor [21]. Figure 2. Morphology of Isoxys. (a) line drawing of idealized Isoxys illustrating known soft parts. (b) YPM IP 005804, Isoxys acutangulus from the Burgess Shale, British Columbia, Canada (Cambrian: Wuliuan) (credit: W. K. Sacco). (c) USNM PAL 189170, Isoxys longissimus from the Burgess Shale, British Columbia, Canada (Cambrian: Wuliuan). Image courtesy of the Smithsonian Institution (EZID: http://n2t.net/ark:/65665/m372f2a644-97c3-441c-87e2-24b1dccb2e8c, credit: Xingliang Zhang). Abbreviations: as, anterior spine; en, endopod; ex, flap like exopod; ey, eye; g, gut with paired diverticulae; ps, posterior spine; ra, raptorial appendage. (Online version in colour.)
Here, we provide the first quantitative assessment of carapace shape variation across the group and compare Isoxys to other Cambrian ‘bivalved’ euarthropods and Gnathophausia. Subsequently, through computational fluid dynamics (CFD) simulations, we test the importance of the spines and carapace outline for generating lift and reducing drag, and thus the ability of different taxa to move vertically through the water column. These analyses support the hypothesis that Isoxys taxa occupied a variety of distinct niches in Cambrian oceans, including vertical migrants. Two-dimensional outlines of 20 Isoxys species, its sister taxon Surusicaris elegans [22], 6 Tuzoia species and 11 Gnathophausia species were constructed in Inkscape from the literature sources by S.P. (electronic supplementary material, table S1) and imported into R [23] for elliptical Fourier analysis (EFA). D.A.L. independently constructed 20 Isoxys and 1 Surusicaris outlines directly from photographs of fossil specimens to allow assessment of the error introduced in the creation of outlines (electronic supplementary material, table S2; and S1). Tuzoia was chosen for comparison as it is also common in Cambrian communities, and has been suggested to be closely related to Isoxys based on similarities in the structure of the eyes and carapace shape (e.g. [24]). Gnathophausia was selected because similarities in the carapace morphology of one species (G. zoea) and Isoxys have been repeatedly noted (e.g. [16,19]), the carapaces of these animals are similar in size (10–30 mm), and multiple species of Gnathophausia are known to be vertically mobile in the modern ocean, having been sampled from surface waters and at depths of over 3000 m [25]. Outlines were sampled at the same resolution (64 points provided sufficient detail to distinguish taxa), centred, scaled by centroid size and subjected to EFA using the Momocs package [26]. Harmonics describing 99.9% of the variation were retained. EFA results were visualized with a principal components analysis. A hierarchical clustering analysis (cluster package; [23]) quantitatively grouped similarly shaped carapaces together using all principal components. Isoxys species reflecting the variation in carapace shape over both PC1 and PC2 were chosen for inclusion in CFD simulations, and Gnathophausia zoea was analysed for comparison. We chose to undertake a two-dimensional analysis to conserve computational resources and minimize errors in the modelled geometries (the exact three-dimensional shape is unknown for most taxa). This is justified because undeformed Isoxys specimens preserved in dorsal view show a narrow profile [16,21,27,28]. Propulsion during swimming derived from movement of the ventral appendages, and not from flapping of the bivalve carapace [15,27]. A lack of adductor muscles means that Isoxys was unable to alter the size of the gape during swimming [15], and the numerous specimens preserved in ‘butterfly’ orientation are considered exuviae [27]. In addition, the full variation of the carapaces considered in this study can be visualized in two dimensions. Two-dimensional analyses are commonly performed on analyses of wing outlines to assess aerodynamic performance by both biologists and engineers, (e.g. [29,30]). A two-dimensional analysis is suitable for this study as the wake behind the Isoxys carapaces is steady at the Reynolds numbers considered, so can be modelled with two-dimensional simulations [31,32]. The impact of soft parts such as eyes protruding from the carapace was analysed for two taxa: Isoxys acutangulus and I. longissimus. Isoxys was assumed to be negatively buoyant, like Gnathophausia [33]. All Isoxys species were assumed to have the same carapace composition and density. The cuticle ornamentation in some Isoxys species is not expected to impact the drag at the low Reynolds numbers considered in this study, as the roughened surface falls within the slowly moving fluid near the carapace surface [31]. Following validation and verification of the model and set-up for low Reynolds numbers, and mesh quality assessments using ANSYS Mesh and ANSYS Fluent (Ansys Academic, Release 2020 R2; electronic supplementary material, S2), outlines of the selected Isoxys species and G. zoea were exported from R as .txt files readable by ANSYS DesignModeler (Ansys Academic, Release 2020 R2; electronic supplementary material, S2), and standardized to a dorsal length (chord length) of 25 mm. This allowed size-independent comparisons of hydrodynamic performance of shapes. While some Isoxys taxa (e.g. I. communis, I. longissimus) can reach up to 50 mm, other adult forms only reach approximately 20 mm (e.g. I. glaessneri, I. volucris) [18,27]. A size of 25 mm represents a compromise size for comparison between these larger and smaller forms, with the influence of larger size able to be assessed by simulating higher Re (as Re depends on both size and swimming speed). CFD simulations were conducted using the steady-state laminar solver in ANSYS Fluent (Ansys Academic, Release 2020 R2). The laminar solver performed best of the three considered during validation and verification (laminar, SST, k-epsilon; electronic supplementary material, S2), as expected for the low Reynolds numbers in this study [31]. Coefficients of drag (Cd) and lift (Cl) (electronic supplementary material, S2) were calculated under three flow speeds equating to 0.75, 1.00 and 1.18 body lengths per second (chord Reynolds numbers, Re, 255, 340 and 400 respectively for saltwater conditions at 0°C; electronic supplementary material, S2). These Re were chosen as swimming speeds of between 75% and 100% of the body length per second have been observed in adult Gnathophausia ingens (carapace mean length approx. 25 mm [34]), and as the laminar model was validated against published drag and lift data for NACA aerofoils at exactly Re = 400 [29]. The chord length (25 mm) was taken as the reference area for both coefficients. Solutions were considered converged when residuals were less than 10–6. Simulations were run at numerous angles of attack, to evaluate the hydrodynamic performance of carapaces at multiple orientations. In tank experiments, Gnathophausia ingens has been observed to change angle of attack to generate more lift or less drag at different swimming speeds [34]. The angle of attack was increased from 0 to 8° at all Re, until the stall angle could be identified and/or unsteady flow was observed. If the stall angle was not reached, further experiments were run until the maximum lift coefficient was obtained. Negative angles of attack were also simulated to assess the negative lift generated by the different outlines. In all cases, the absolute value of the negative angle of attack was increased until the drag coefficient was equal to or greater than that of the stall angle. When unsteady flow was suspected to be the reason that steady-state simulations did not converge, the steadiness of the flow field was determined by carrying out a time-dependent analysis of 100 time steps, with each time step equal to the flow speed (so, for an inlet velocity of 0.01875 ms–1, the time step = 0.01875 s). In the outline analysis, a total of 18 harmonics were retained. Principal coordinates 1 and 2 described 79.3% of the total variation. Carapace asymmetry, narrowness and spine length increased as PC1 became more positive, while the length of the anterior spine relative to the posterior spine corresponded to an increase in PC2 (figure 3). Isoxys occupied the largest area in the morphospace. Visual overlap of the areas occupied by the genera demonstrated that some Isoxys taxa were more similar in shape to Gnathophausia than their Cambrian relatives Surusicaris and Tuzoia. Confirmation was provided by the cluster analysis (figure 3). All six Tuzoia species formed a single cluster, with Isoxys taxa spread over the three remaining groups (Surusicaris group, Gnathophausia groups 1 and 2). Species clustered with Surusicaris have symmetric and deep carapaces and relatively short spines. Species in the Gnathophausia groups displayed narrower outlines whose narrowness, asymmetry and spine length increased with PC1. The single species of Isoxys in Gnathophausia group 1, I. paradoxus, displayed an anterior spine much more elongate than its posterior one, that contrasted with the seven species in Gnathophausia group 2 with their spines of approximately equal length. Figure 3. Principal component analysis of results of EFA conducted on the outlines of 11 species of Gnathophausia, 20 Isoxys, 1 Surusicaris and 6 Tuzoia. Convex hulls indicate optimum four groupings as recovered by clustering analysis. Labelled Isoxys species chosen for subsequent hydrodynamic analysis. (Online version in colour.)
CFD simulations assessed the impact of increasing asymmetry, spine length and relative lengths of anterior and posterior spines on hydrodynamic performance. Inclusion of eyes did not significantly impact the hydrodynamic performance of carapaces (electronic supplementary material, S3). Isoxys zhurensis, the most symmetric species chosen for analysis, created an unsteady wake with Kármán vortex street at the lowest Re considered, and so no drag or lift coefficients were obtained (electronic supplementary material video). The flow around remaining Isoxys outlines was laminar with a steady wake, and there was no evidence for three dimensionality (electronic supplementary material, S4). Greater asymmetry and narrowness (more positive in PC1, figure 3) resulted in lower drag coefficients, as demonstrated by a comparison of the short-spined taxa I. chilhoweanus, I. acutangulus and I. mackenziensis. The most asymmetric of these forms, I. mackenziensis, produced lower drag coefficients than the other two, but ranges of lift coefficients were similar for all three (figure 4a). More elongate spines increased the range of lift coefficients (vertical bars, figure 4a) and, significantly, negative lift coefficients at negative angles of attack (e.g. compare I. mackenziensis, I. communis and I. longissimus). In I. paradoxus, where the anterior spine is much longer than the posterior one (more positive in PC2, figure 3), the range of lift coefficients further increased (figure 4a). Similar drag coefficients and ranges of lift coefficients were obtained in an analysis of the hydrodynamics of Gnathophausia zoea, when either the carapace alone or both the carapace and abdomen were considered (figure 4b). Figure 4. Drag polars (plot of Cd against Cl) of taxa analysed at Re = 255 (0.75 body lengths per second for an animal 25 mm long). Each point corresponds to a single simulation at a different angle of attack. Vertical bars show range of lift coefficients. Note that flow was unsteady for Isoxys zhurensis at Re = 255, and so no quantitative lift or drag coefficients were recorded. Drag polars at faster flow speeds and raw data presented in electronic supplementary material, S3. (Online version in colour.)
Functional morphology of Isoxys fossil specimens supports an off-bottom (hyperbenthic or pelagic) life habit for this animal [15,16,18,19,32], based on the eye orientation (forwards, slightly ventral) and the elongate slender carapace shape of Isoxys. Our study combines outline morphometric and CFD analyses and suggests that Isoxys species occupied a variety of niches, including some as pelagic vertically mobile predators. Lift and drag coefficients of Isoxys carapaces indicate variation in the depth range and swimming speeds of these species. Isoxys taxa clustering with Surusicaris (figure 3) generate positive lift, but do not generate negative lift (figure 4). This supports suggestions of previous workers that these Isoxys species may have occupied a hyperbenthic (1–10 m above the seafloor) and/or nektobenthic [21] niche, perhaps moving vertically short distances in the water column [15,19,35]. Vertical movement would be achieved by altering the angle of attack to produce lift force greater than (ascent), equal to (horizontal swimming), or less than (descent) the impact of their negative buoyancy. Drag reduction associated with streamlining would have allowed some taxa (e.g. I. mackenziensis) to capture faster-moving prey animals, as faster swimming speeds could have been maintained over longer distances for the same metabolic cost. Isoxys species clustering with Gnathophausia show convergent adaptations to moving vertically in the water column, such as asymmetric carapaces with elongate anterior and posterior spines (figures 3 and 4). This does not preclude elongate spines from also acting as anti-predatory deterrents, as suggested by Vannier & Chen [15]. These adaptations provided hydrodynamic benefits that would have allowed the Isoxys species to operate over a wider bathymetric range. A streamlined carapace facilitates not only faster movement, but also more efficient swimming, beneficial for migrations over a long distance. The carapace shapes of I. longissimus and I. paradoxus generate lift coefficient ranges and minimum drag coefficients comparable to the modern crustacean Gnathophausia zoea (figure 4), which has been recovered at depths ranging from surface waters down to approximately 3000 m in the modern ocean [25]. These results also suggest that an elongate abdomen in G. zoea reduces the drag experienced by the animal slightly but does not greatly impact on the range of lift coefficients (figure 4b), though the abdomen may also play a physical role. Animals that move vertically in the water column do not have to cover the entire distance from the surface ocean to demersal communities, and instead sometimes migrate across only a shorter vertical distance. Thus, Isoxys taxa with the broadest ranges of lift coefficients (I. longissimus and I. paradoxus) probably covered a wider depth range than those with smaller ranges of lift coefficients (e.g. I. communis). Corroborating evidence for variation in bathymetric range for different Isoxys species comes from the fossil record itself. Members of different groups as resolved in the cluster analysis (convex hulls, figure 3) co-occur with different relative abundances in Cambrian deposits preserving soft-bodied fossils. In general, species with inferred vertically migrating lifestyles are much rarer than those that lived close to the seafloor. In the Chengjiang Biota, Isoxys auritus (Surusicaris group) greatly outnumbers both I. paradoxus and I. curvirostratus (Gnathophausia groups 1 and 2, respectively) (figure 3; [15,21,36–38]). A similar pattern of relative abundances can be observed in the two Burgess Shale taxa (figure 2): in the Walcott Quarry, I. acutangulus (Surusicaris group) comprises nearly 0.5% of the total community, vastly outnumbering the extremely rare inferred vertical migrant I. longissimus (Gnathophausia group 2 [16,39]). The relative abundances of these Isoxys species can be partly explained by the differences in lifestyle predicted by the carapace outline and soft anatomy. The less hydrodynamic taxa (those with higher drag coefficients and narrower ranges of lift coefficients) probably lived near the seafloor, with the more streamlined species living in the water column and occupying a broader bathymetric range. This broader bathymetric range would have included more open water settings, beyond the maximum depth of the shelf where Cambrian deposits preserving soft-bodied fossils occur—Gnathophausia zoea for example has been found at depths of up to 3000 m [25]. As modern euarthropod carapaces disarticulate quickly after death (e.g. [40]), the preservation potentials for pelagic euarthropods living high in the water column are lower than for those living closer to the seafloor. In addition, animals which occupy an ecological niche in the open water are less likely to find themselves over shelf environments such as those which preserve soft-bodied fossils or be trapped and transported by an obrution event responsible for the preservation of soft-bodied communities in these settings. The small numbers of vertically mobile Isoxys individuals observed may have been at the bottom of their vertical migrations and/or been transported horizontally by currents. Isoxys species are not globally distributed [41]. Many species (for example those clustering with Surusicaris) appear suited to hyperbenthic habits, and so would be expected to have provincial distributions. The limited geographical distribution of I. longissimus and I. paradoxus is most likely to be due to a combination of factors. Firstly, deposits where Isoxys is expected to be preserved are not evenly distributed in time and space—Stage 3 deposits are mostly in South China, while Wuliuan and younger are mostly in Laurentia [42], though the absence of the Chengjiang species I. paradoxus from Sirius Passet is notable. Secondly, the lower preservation potential of pelagic (compared to hyperbenthic) species means that they are rare even in Tier 1 Burgess Shale-type Lagerstätten (sensu [42]). However, despite its rarity, the Burgess Shale species I. longissimus has a wider known geographical range than the co-occurring I. acutangulus. The former has also been reported from the Wheeler Formation, House Range, UT, USA [43]. The situation appears more complex in the Emu Bay Shale, where the more hydrodynamic species I. communis greatly outnumbers the less streamlined I. glaessneri [27]. However, the Emu Bay Shale is not a traditional Burgess Shale-type deposit, instead it represents a localized deep-water micro basin on the inner shelf [44]. Here fluctuating oxygen levels may have periodically deoxygenated the water column, possibly killing pelagic taxa like I. communis in great numbers and creating a taphonomic bias that preferentially preserves pelagic taxa. Further support for the Chengjiang taxon Isoxys auritus occupying a niche closer to the seafloor than I. curvirostratus comes from a comparison of the soft anatomy (soft parts are unknown for I. paradoxus) [21]. The stout endopods of I. auritus appear well suited for interacting with the substrate, while exopods with broad fringing lamellae and a sophisticated vascular system in the more streamlined I. curvirostratus suggest it was a more powerful swimmer, providing additional support for a pelagic habit [21]. Lastly, a compendium of fossil, geochemical and phylogenetic data show that vertically mobile Isoxys species would have had access to a variety of pelagic prey items and an oxygenated water column to travel in. Cambrian oceans were not stratified, instead displaying wedge-shaped oxygen minimum zones broadly comparable to modern oceans [45]. Isoxys prey size range (approx. 5–20 mm [19]) includes Cambrian phytoplanktivores, such as crustaceans with setae and filter plates from Sirius Passet and Mount Cap (15–50 mm) [4,8,9,12], as well as crown-group branchiopods and copepods and total group ostracods [46], as well as bradoriids, some of which were also likely pelagic [47,48]. The presence of planktonic larvae, another possible prey item, can be inferred from tip-dated phylogenetic analyses that support the evolution of metamorphosis in euarthropods by the Cambrian [49]. Thus, a range of different data sources suggest multiple Isoxys taxa were vertically mobile, and that the Cambrian ocean could support such an ecology. The presence of the likely vertically mobile Isoxys paradoxus in the Chengjiang Biota makes it the oldest confidently identified euarthropod vertical migrant, and probably among the first animals to employ this life habit. For this vector to be significant by the Cambrian Stage 3, a large biomass of Isoxys would need to move vertically. While pelagic animals have a lower preservation potential than benthic ones (for example, very few fossil copepods are known [50]), Isoxys species with inferred (hyper)benthic habits are extremely abundant in both the Chengjiang and Burgess Shale [15,16,21,36–38]), and their pelagic counterparts Isoxys longissimus and I. paradoxus may have been similarly numerous. The Chengjiang and coeval Qingjiang biotas also preserve the earliest evidence for gelatinous zooplankton, which move vertically small distances in the modern ocean [51,52]. However, an active swimming Isoxys would have covered greater distances more rapidly. Furthermore, the presence of a through gut would have increased processing time for food, vital for transporting nutrients consumed in surface waters to the deep ocean as faecal pellets. Vertical migration was one of many important eukaryotic and metazoan innovations key to establishing a modern-style biological pump. A series of metazoan innovations which appear in the fossil record in quick succession during the early Cambrian gave the biological pump a modern-looking structure (figure 5), which was strengthened during the Phanerozoic. Figure 5. First known appearance in the fossil record of metazoans that impacted the biological pump. Circle colour denotes preservation style. Fossils with more than one preservation style indicated with split circles, with the preservation style that provided the oldest evidence on the left of the circle. Abbreviations: Drum., Drumian; Guz., Guzhangian; Wul., Wuliuan. (Online version in colour.)
The shift to primary production dominated by eukaryotes, and the innovation of active suspension feeding, would have ventilated the oceans, cleared organic matter in the water column, and increased transfer of oxygen and plankton of increasing sizes from the surface to the sediment–water interface (e.g. [6]). The first of these events occurred during the Cryogenian [7], but while benthic passive suspension feeders are known from the Ediacaran (e.g. [53]), early Cambrian stem and crown-group sponges represent the oldest benthic active suspension feeders [10,11]. As active suspension feeders, sponges were able to transport large volumes of water between benthic and pelagic realms (e.g. [54]), and the most abundant sponges from the Cambrian, the reef building archaeocyaths, display pore size differentiation within reef systems [55], presumed to be evidence of prey size-selectivity. This illustrates that predator–prey feedbacks were present in the Cambrian Stage 2 (figure 2), presenting a potential driver towards a larger size of plankton, increasing the sinking speed and efficiency of this part of the biological pump, and ventilating the water column. The invasion of the pelagic realm by eumetazoan zooplankton provides the next step towards a modern-style biological pump. These zooplankton would have further cleared surface waters and contributed faecal pellets to organic aggregates sinking to the deep ocean, and also increased oxygen levels at depth (e.g. [6,54]). The small shelly fossil (SSF) record provides a source of evidence for the invasion of the pelagic realm by eumetazoans. Terreneuvian SSFs include the possible chaetognath Protoherzina and the molluscs Watsonella, Aldanella and Oelandiella, whose widespread distributions are suggestive of a pelagic lifestyle, or at least a planktonic larval stage (figure 5) [56]. Euarthropods, likely early occupants of the plankton [35], are represented in the SSF record from the Cambrian Stage 3 by millimetre-scale bradoriids among others (figure 5; e.g. [47,57,58]). The first macroscopic nektonic suspension feeders, such as the radiodont Tamisiocaris, also appear at this time [59], while the first centimetre-scale phytoplanktivores are identified close to the Stage 3–4 boundary (figure 5) [12]. These data suggest that there was an increase in the diversity of millimetre-scale zooplankton at or close to the base of Stage 3, very close in time to the appearance of the first vertical migrants. Most bradoriids are considered benthic, however Anabarochilina increased its distribution in three phases, providing complementary evidence for a steady strengthening of the pump during the Cambrian. In Epoch 2 Anabarochilina was coupled with benthic assemblages, by the Wuliuan it spread to a wider spectrum of lithofacies, and by the Guzhangian two species became widely distributed [48]. The metazoan innovation of vertical migration would have impacted both demersal and pelagic communities. The strength of the impact depends on the amount of biomass undertaking vertical migration. Models based on Cambrian environmental parameters predict that vertical migration would have increased the efficiency of the carbon pump by around 7% [60], however, more significantly, vertical migrants transport organic nutrients to the deep sea more quickly than aggregates, with a different nutritional balance, and repackage decaying sinking organic matter (figure 1) [1,13,61–64]. In addition, vertical migrants are major contributors to ocean mixing and ventilation, spreading oxygen and nutrients throughout the water column (e.g. [54,65]). These effects probably played a role in contributing to the rapid rate of diversification during the Cambrian explosion, interwoven with numerous evolutionary and ecological feedbacks. For example, the higher metabolic needs and nutrient requirements of large biomineralizing animals and motile predators [66] may have been facilitated by increased quality and quantity of nutrient transport in the biological pump, and resulted in increased oxygenation of bottom waters. In turn, the increasing size and motility of predators would have provided a further ecological pressure for animals to ‘escape’ into the pelagic realm. The establishment of a biological pump with a modern-style architecture by the Cambrian Stage 3 does not mean that the fluxes along the aggregation and vertical migration vectors (figure 1) remained constant to the modern day. Indeed they likely strengthened through the Palaeozoic with an increase in biomass (from an increased number of taxa, individuals and size of individuals). Fossil evidence points to additional metazoan innovations during the Palaeozoic that would have affected the fluxes along these vectors and strengthened the pump. For example, the aggregation vector would have been strengthened following the evolution of centimetre-scale and decimetre-scale phytoplanktivores later in the Cambrian [12,67], and the major radiation of plankton and the evolution of metre-scale nektonic suspension feeders during the Great Ordovician Biodiversification Event [68–70]. The flux of nutrients along the vertical migration vector would have increased as pelagic and vertically migrating animals diversified and increased in size—an innovation that would also have increased the mixing of waters and ocean ventilation. For example, the evolution of large, fast moving fish as part of the Devonian nekton revolution is expected to have been especially significant [54,71]. In summary, the innovation of vertical migration in some Isoxys species was one of several interwoven and coevolutionary feedbacks during the early Cambrian that increased the efficiency and altered the architecture of the biological pump, likely contributing to the rapid expansion in metazoan diversity at this time. All data and supplementary materials are available through the Open Science Framework: https://doi.org/10.17605/OSF.IO/2JDRS. S.P.: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, validation, visualization, writing–original draft, writing-review & editing; A.C.D.: writing–review & editing; D.A.L.: data curation, writing–review and editing; I.A.R.: methodology, validation, writing–review & editing. All authors gave final approval for publication and agreed to be held accountable for the work performed therein. We declare we have no competing interests. S.P. was supported by an Alexander Agassiz Postdoctoral Fellowship (Harvard University) and a Herchel Smith Postdoctoral Fellowship (University of Cambridge), D.A.L. by a Dame Kathleen Ollerenshaw Research Fellowship (University of Manchester), and I.A.R. by a Museum Research Fellowship (Oxford University Museum of Natural History). We thank the associate editor, two anonymous referees, and Christian Klug, who provided helpful reviewer comments. We thank members of the Ortega-Hernández Lab for Invertebrate Paleobiology (Harvard University) for fruitful discussions. S. Butts (Yale Peabody Museum) and M. Florence (Smithsonian National Museum of Natural History) provided curatorial assistance. Footnotes© 2021 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 9Currently, the oceans provide about half of the global net primary production (NPP), approximately 48.5 Pg C yr−1 [1]. The fate of most primary production in the ocean is remineralization and rapid recycling of nutrients and CO2 in the upper ocean [2]. On average, only 20% [3] is exported in the form of particulate organic matter into the ocean interior. There, approximately 97–99% [4] is remineralized, but now spatially separated from the ocean surface and atmosphere. A return to the surface from the deep ocean depends on the much slower action of ocean circulation and upwelling. A small fraction of particulate organic matter and hence carbon also escapes the ocean-atmosphere system and is buried in marine sediments, helping regulate atmospheric pCO2 on geological time scales. This dynamical biogeochemical partitioning between ocean surface and interior, particularly of carbon, is known as the marine ‘biological pump’. This process influences not only nutrient availability and primary production, but also atmospheric pCO2, and hence climate, among the fundamental functions of the marine pelagic ecosystems. The strength (magnitude) and the efficiency (remineralisation depth) of organic matter transfer by the biological pump is highly variable [5]. This variability is largely a result of differences in phytoplankton, zooplankton and microbial communities and food web structures [4,5]. For example, the proportion of large cells within the phytoplankton community influences the flux of particulate organic carbon (POC) [6], while the type and nature of zooplankton feeding can lead to an increase or decrease in the sinking rate of organic matter [4]. Understanding the controls on both primary production as well as the efficiency of transfer of carbon to the ocean floor is hence essential for an understanding of the regulation of pCO2 and climate. Model projections show marine primary production decreasing globally due to climate change [7,8] particularly in the tropics and the North Atlantic, though uncertainties are high [7,9]. Environmental changes in temperature and CO2 concentration impact the ratio of calcifiers to non-calcifiers (potentially changing the composition and hence sinking rate of particles), as well as the proportion and size of plankton and grazers [10–13]. Modelling studies reveal that under projected future climate change, the decrease in phytoplankton size exerts the largest influence on decreasing POC flux out of the surface ocean, whereas changes in zooplankton communities are important in the subsurface ocean [14]. Assessing these interactions is essential for understanding the critical elements in, and the resilience of, the biological pump to perturbations. Observational time series are not long enough to reject or support these projections [15]. The longer perspective of the fossil record, however, could fill this gap as well as providing crucial information to address potential impacts of marine extinctions, given the anthropogenic effect that humans are having on Earth's processes in the wake of a potential 6th mass extinction [16,17]. As isotopically light carbon (12C) -CO2 is preferentially taken up during primary productivity in the surface ocean, the carbonate shells mineralized by planktic foraminifera are relatively enriched in 13C compared to the organisms that precipitate carbonate in the deep ocean (e.g. benthic foraminifera), where organic matter is remineralized releasing 12C, resulting in a vertical carbon isotope (δ13C) gradient [18]. The greater the surface-deep δ13C difference the more efficient the biological pump. This signal of biological pump efficiency (which differs from strength, see [19] for details) is recorded in the shells of foraminifera and combined in the geological record with information on the ecological composition of planktic ecosystems. The Cretaceous/Palaeogene (K/Pg, 66.02 Ma) mass extinction [20] was the most important extinction in the evolutionary history of modern plankton [21] and provides an excellent opportunity to understand how a severe perturbation to the planktic ecosystem affects the biological pump. Based on interpretations of a surface to deep ocean δ13C gradient of near zero [22], it has been suggested that the extinction led to a near complete shutdown of the biological pump, with a recovery time of several millions of years [22–24]. However, a lack of extinction within benthic communities, which depend on surface-exported food supplies, challenges this interpretation (see [22]). While some proxies, such as biogenic barium [25], suggest spatial heterogeneity in productivity between the major ocean basins and open ocean versus shelf environments, the δ13C signal is a global phenomenon [23,26]. Previous studies have focused on either the recovery of the marine biological pump [22,24,27] or the pelagic biota [28,29]. Few studies have tried to establish the link between the two [30–32]. While ideally the whole ecosystem would be interrogated, most pelagic organisms do not preserve in the fossil record [33]. Fortunately, two important autotrophic and heterotrophic components have excellent fossil records i.e. calcareous nannoplankton (haptophyte algae) and planktic foraminifera (shell-building micro-zooplankton) that provide representative signals of the pelagic ecosystem. Differences in rates of recovery of the biological pump, plankton diversity and size are apparent in the fossil record [30,31,34–37]. Recently, post-K/Pg ecosystem recovery, assessed via community stability history of nannofossils, has been linked to a return of the biological pump in the Pacific, approximately 1.8 Myr years after the event [31]. Here, we assess whether the marine biological pump drove ecological changes, or, whether the recovery of the biological pump was itself contingent on pelagic community recovery (diversity, size) to fulfil their ecological function. We explore whether ecosystem function (i.e. the contribution of marine plankton in the regulation of the global carbon cycle) depends on diversity to be re-established. Specifically, we question how important diversity recovery at higher trophic levels is for the efficiency of the biological pump. Moreover, did the evolution of certain traits in plankton, such as body size and photosynthesis, drive the restoration in biological pump efficiency, or do certain environmental and ecological thresholds need to be met for these traits to be established? Here, we focus on planktic foraminifera, which benefit from a well-established understanding of both modern and paleo-diversity and ecological preferences [38], and several methods for documenting the evolution of body size [35,39]. The K/Pg event is captured in Ocean Drilling Program (ODP) Site 1262, Walvis Ridge (27°11.15' S and 1°34.62′ E; electronic supplementary material, figure S1). The K/Pg boundary occurs at approximately 216.6 m composite depth (mcd), calibrated to 66.02 Ma on an astronomically tuned time scale [40]. Consistent preservation of calcitic microfossils suggests deposition above the carbonate compensation depth throughout the K/Pg interval [41]. Core samples were washed over a 38 µm sieve and dried in an oven at 40°C. Planktic foraminifera species abundance counts were made on 49 samples. Taxonomy follows the Palaeocene Atlas [42] and Cretaceous chapters of Plankton stratigraphy [43]. The summed coefficient of variation metric (∑CV), which quantifies the level of stability, was calculated per sample for the five designated ecogroups, clusters of species sharing similar ecologies, following the method of [44], without the SiZer smoother step, as in [31]. ∑CV was chosen as it is independent of taxonomic composition. Separated sieved size fractions (electronic supplementary material, table S1) were weighed in 94 samples. Samples were weighed using an A & D semi-microbalance (standard deviation of 0.1 mg). The species/genera counts were assigned to one of five ecogroup categories for the Palaeocene and one of four categories for the Cretaceous (electronic supplementary material, table S2). A representative split of 356 washed samples (greater than 38 µm) was analysed for foraminifera size. Size parameters of randomly oriented foraminifers were measured using a Malvern Mastersizer laser granulameter. The maximum diameter of the object was chosen as the most suitable size estimator because it is least affected by random orientation. The 90th percentile of the maximum diameter (or D90) was used to describe these strongly skewed distributions. Our data show that overall planktic foraminifera diversity remained low compared to pre-extinction assemblages for approximately 4.8 Myr post-extinction (figure 1b). Test size dropped dramatically from approximately 400 to 150 µm at the K/Pg and small size persists throughout the studied interval (figure 1e). Linked to this, the relative proportion of particles greater than 38 µm increased indicating a higher contribution of planktic foraminifera to bulk carbonate (figure 1f). Dividing the greater than 38 µm foraminiferal fraction further (electronic supplementary material, table S1) shows that most size classes contributed roughly equally to the assemblage in the Cretaceous (approx. 5.5–15%). After the K/Pg boundary, approximately 75% of the total foraminiferal fraction were in the size fraction below 106 µm (figure 1g). Consequently, carbonate accumulation rates dropped strongly, enhanced by the mass extinction of calcareous nannofossils, which typically contributed 80–90% of the bulk dry weight before the extinction (figure 1h,i, respectively). However, the abundance of small opportunistic planktic foraminifera increased (figure 1g) and resulted in an increased contribution of foraminifers to bulk carbonate (figure 1j). Between the partial and full recovery in the marine biological pump (300 kyr to 1.8 Myr later; line 1 to 2) many of the sedimentary parameters (% carbonate, F/N ratio and foram fraction) began to stabilize. Alpha diversity increased after approximately 1.8 Ma and approximately 4 Ma, indicating that originations exceeded extinctions for a brief period (figure 1b). These diversity increases were associated with minor increases in foraminifera size (figure 1e). Figure 1. Multiproxy K/Pg datasets from Site 1262. (a) Adjusted stable carbon isotope (δ13C) differences between planktic and benthic foraminifera [22]. (b) Shannon diversity of planktic foraminifera greater than 106 µm (black line) and total number of species greater than 106 µm (red line). Blue and red bars indicate intervals of enhanced planktic foraminifera species origination and extinction[34], respectfully, and vertical dashed lines show minimum K/Pg value. (c) Relative abundance (%) of Palaeocene planktic foraminifer ecogroups. (d) Summed coefficient of variation metric (∑CV). (e) 90th percentile foraminifera size (µm) [35]. (f) weight contribution (g) of the foraminiferal fraction (i.e. coarse fraction greater than 38 µm) to the bulk sediment. (g) Planktic foraminifera test size spectrum shown as cumulative weight % (g) of size fractions. (h–j) The % carbonate, carbonate accumulation rate and foraminifera/nannofossil accumulation ratio [45]. Ages are based on the tuning of the records to La2011 [40]. Horizontal black dashed line represents the K/Pg boundary and coloured numbered lines relate to significant steps in the recovery process; 1—initial partial recovery of the surface-to-deep carbon gradients, 2—‘full’ carbon gradient recovery, 3—appearance of photosymbiotic ecologies [22,46] and 4—‘full’ planktic biotic recovery. (Online version in colour.)
The survivor and opportunistic species dominated the initial recovery interval (approx. 100 ka). These were followed by ‘transitional’ taxa (figure 1c) whose early representatives were deep dwellers while descendent species migrated to surface waters and/or became symbiotic. The transitional taxa declined and were replaced by thermocline species around 1 Ma (figure 1c). A surface/symbiotic group appeared approximately 3 Myr later (figure 1c) when size spectra widened and small taxa and individuals lost their dominance (to approx. 50%, figure 1g). The balance between the relative contribution of the two carbonate producing groups (figure 1f,j) returned to pre-extinction levels also by approximately 3 Ma. ∑CV shows a small peak (approx. 3.5 Ma) above background variation (electronic supplementary material, figure S3), which coincides with the decline and extinction of the transitional taxa and the early radiation of the surface/symbiotic group. This group only became an important component of the assemblages approximately 4.3 Myr after the extinction, when assemblage ecological characteristics were restored but with completely different species. While size classes diversified, they did not reach the overall sizes of the Cretaceous. The impact of an asteroid at the K/Pg had devastating effects on Earth's fauna and flora and created environmental instability [20]. Extreme temperature changes, reduced pH, heavy metal loading, stratification and increased nutrients are all thought to have played a part in the initial extinction and dictated what organisms survived [20,32,47]. Environmental conditions feedback on the diversity and abundance of marine organisms [48–50] and thereby change the fixation and export of carbon and utilization of nutrients. As marine ecospace was re-created [28,36,51], rates of evolutionary turnover were far above typical background rates [34]. Small, opportunistic planktic foraminifera dominate the early Danian [29] and a high frequency of morphological abnormalities in Tunisian planktic foraminifera has been linked to severe environmental instability [52]. Successive acmes of opportunistic, eutrophic calcareous nannofossils are attributed to high-nutrient loads [31,37], as are increases in high-nutrient, opportunistic dinoflagellates, indicative of eutrophication in marginal settings [53]. Blooms of neritic opportunistic dinoflagellates declined by the end of biozone P1a [53], approximately 50 kyr after the extinction. Planktic foraminiferal abnormalities reduced to background levels 200 kyr after the K/Pg [52]. Disappearance of these bloom forming, high-nutrient taxa suggests that the reduction of nutrient levels in marginal surface waters predates the initial return of export productivity observed in the open ocean (line 1, figure 1a) [22] supporting suggestions that neritic and marginal marine environments recover quicker (<10 kyr) than open ocean environments [54]. The microperforate group (which is similar to our opportunistic/survivor group) at the shallow water impact Site M0077 [37] declined in abundance from approximately 1.2 Ma ending near the return of the biological pump (figure 2). However, no notable biotic change, neither in traits such as size nor diversity, is observed in our open ocean region with the full return of the biological pump at approximately 1.8 Ma (line 2, figure 1a). The succession of acmes of nannofossil species in the Pacific [31] ends at this time (figure 2), yet, pre-extinction nannofossil size and ecogroups are not re-established until 3.5 Myr later, and diversity not until 10 Myr later. Figure 2. Multiproxy K/Pg datasets. (a) Adjusted stable carbon isotope (δ13C) differences between planktic and benthic foraminifera for Site 1262 (see [22]). (b) Relative abundance (%) of Palaeocene planktic foraminifer ecogroups for Site 1262. (c) Nannoplankton community variance; Simpson's Diversity Index ODP Site 1209, Global species richness from [31]. (d) Ichthyolith accumulation rates; Tethys Sea (Gubbio, Italy; Icth. >38 µm/cm−2 Myr−1) from [55] and Walvis Ridge, DSDP Site 527 in the South Atlantic (Icth. >63 µm cm−2 Myr−1) from [56]. (e) Percentage planktic groups from Site M0077 [37]. Ages are based on the tuning of the records to the latest orbital solution La2011, [40]. Horizontal black dashed line represents the K/Pg boundary and coloured numbered lines relate to significant steps in the recovery process; 1—initial partial recovery of the surface-to-deep carbon gradients, 2—‘full’ carbon gradient recovery, 3—appearance of photosymbiotic ecologies (see [22] for details) and 4—‘full’ planktic biotic recovery. (Online version in colour.)
Symbiosis is widespread in modern planktic foraminifera and allows populations to thrive in low-nutrient, oligotrophic environments [38]. The reacquisition of photosymbiosis [46] in planktic foraminifera did not lead to rapid diversification (figure 1d; electronic supplementary material, figure S2) [30,57]. We speculate that surface waters during the earliest recovery may have been rich in nutrients due to lower consumption, such that this ecology was not selected for. This speculation is supported by independent records of calcareous nannofossils which also do not begin to show a return to more oligotrophic taxa until 3.5 Myr later [31,37]. The contribution of the different size classes of foraminifera may give insight into the complexity of the food webs, as organisms tend to eat prey approximately a tenth of their size [58]. The equal contributions of the size classes present before the extinction suggest a complex food web. The dramatic shift to smaller shell sizes after the boundary suggests that larger organisms could not thrive [59,60], as food web complexity reduced, perhaps as mean food size decreased and food chains shortened [20,61,62] an interpretation supported by a decline and long recovery in mid food chain predator fish, based on ichthyoliths (fossil fish teeth and shark dermal scales) (figure 2) [55]. In contrast with the recovery of the marine biological pump, plankton diversity and traits point to a much more drawn out return to pre-extinction states. Ecogroups recover around 4.3 Myr after the extinction, while size and diversity do not return to pre-extinction levels within the studied interval. Calcareous nannofossil diversity does not reach pre-extinction values for more than 10 Myr [34,63] and foraminiferal size even longer [35]. Similar long delays in diversity recovery are also documented for calcareous red algae (approx. 6 Myr to recover [64]), Neoselachian sharks (7 and 10 Myr [65]) and corroborated by other studies showing a macrofaunal genus level diversity lag of approximately 10 Myr [66,67]. The biological pump recovery (especially its strength) is likely driven by aspects of the marine community which do not fossilize, such as many microbes [68]. However, unexpectedly, larger size classes of organisms and specialists appear unnecessary for the recovery of the marine biological pump. Concurrently, the recovery of the biological pump does not result in a recovery of diversity and pre-extinction organism size [31]. The process of biotic recovery may be governed by other factors, such as niche creation, which are suggested to work on much longer time scales [36,67]. This interpretation supports the notion that rather than there being predefined ‘niche’ space (i.e. fixed ecological real estate) available to be filled or vacated at times of diversification or extinction, organisms themselves, and their interactions, ‘construct’ the environment they inhabit [69]. In this study, we speculate that increasing diversity created oligotrophic conditions, that in turn increased niche spaces which selected for novel traits such as symbiosis creating feedbacks between environment (and the biological pump) and evolution. The K/Pg extinction impacted the global biogeochemistry of the oceans and marine life for millions of years. Our data show that biotic recovery, as measured by trophic levels, microplankton size classes, and diversity of planktic foraminifera occurred much later than the re-establishment of a marine biological pump. Diversification of ecology in planktic foraminifers was linked to a return in the dominance of the surface symbiotic species adapted to lower nutrient conditions. These data suggest that a large range of nutrient conditions, including oligotrophic conditions, is necessary for high diversity independent of the marine biological pump recovery. Our findings highlight the need to link climate projections, models of primary production and ecology in both coastal and open ocean environments to improve our ability to project the repercussions of climate-induced extinctions or reorganization into novel environments on marine ecosystems and their services to people [70,71]. Raw data used in this study can be found in the supplementary information file and Pangea. We declare we have no competing interests. The material studied was provided by the Integrated Ocean Drilling Program. The research was funded by a Natural Environment Research Council (NERC) studentship to H.B. and Royal Society University Fellowships to D.N.S. and H.K.C. The work was also partly funded by an ERC grant PAst Links in the Evolution of Ocean's Global ENvIronment and Ecology (PALEOGENiE) awarded to A.R. D.N.S. would like to acknowledge the Royal Society for a Wolfson Merit Award. We would also like to thank two anonymous reviewers for their time and constructive comments. FootnotesElectronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5459523. References
Page 10In order to predict how biodiversity patterns on today's Earth will respond to climate change, the factors that cause biodiversity distributions must be understood [1,2]. Deep-time perspectives can provide novel insights into the controls on biodiversity distribution. By examining biodiversity distributions at times in Earth's history when climate, continental arrangement, and oceanic currents were different than today, ecological hypotheses about the causative mechanisms behind biodiversity distribution and the establishment of modern patterns can be tested [3–7]. However, if palaeontologists wish their data to be informative to those working on the causative mechanisms of modern-day biodiversity patterns, we must first demonstrate that actual biodiversity patterns are preserved in our reconstructions of past ecosystems, and that we are able to overcome the many sampling biases that affect the fossil record (e.g. [8–10]). Dinosaurs are an exceptional model system for studying biodiversity and macroevolution in terrestrial vertebrates. For greater than 150 million years, from the Late Triassic to the end of the Cretaceous, they dominated terrestrial ecosystems, occupied every continent, and radiated into a wide variety of ecological niches. Because of public interest, they are the best-sampled Mesozoic terrestrial vertebrate group and their fossils have been collected for well over 150 years [11–14]. Arguably, the best-sampled part of the dinosaurian fossil record is the Late Cretaceous of the Western Interior region of North America [15–17]. During the Late Cretaceous, North America was divided into two landmasses, Laramidia to the west and Appalachia to the east, by the epicontinental Western Interior Seaway. In the latter stages of the Late Cretaceous (Campanian and Maastrichtian), large-bodied herbivorous niches in Laramidia were dominated by two groups of ornithischian dinosaurs, the hadrosaurs and ceratopsids. The body fossil record of the latter is entirely restricted to North America at this time, with the exception of a single taxon [18,19]. Study of these herbivorous dinosaurs has provided major insight into dinosaur behaviour, palaeoecology, and biogeographic patterns (e.g. [20–32]). Numerous workers have argued for strong faunal provinciality in Laramidia throughout the Campanian and have divided the landmass into northern and southern faunal provinces (e.g. [33–35]). This signal is particularly clear in chasmosaurine ceratopsids, where virtually all species are recognized from either northern or southern Laramidia, but not both [34–36]. This endemicity in ceratopsids is thought to have driven high levels of diversity, underpinning their radiation [34]. Since no geological or geographical barrier has thus far been identified between northern and southern Laramidia [37], the boundary between the northern and southern provinces has been suggested to be related to latitudinal climate, temperature, or rainfall patterns [34,35,38] or was maintained due to competition between local populations [39]. The patterns of apparent provincialism decrease in the Maastrichtian, coincident with overall regression of the Western Interior Seaway [15,17,33,36]. These hypotheses of biogeographic provincialism, however, remain controversial. With very few exceptions [35,36], studies that advocate for provincialism are based on qualitative observations (e.g. [33,34]) and arise from comparisons of the fauna of specific geological formations (e.g. [33–35]). Recent research has, however, suggested that some of the formations used in such studies are not contemporaneous [17,40] and that the length of time intervals used results in the amalgamation of multiple successive faunas [17,37]. Many studies advocating faunal endemism are based on taxonomic decisions that have proven controversial and the conclusions have been called into question as a result (e.g. [37,39]). Additionally, it remains a possibility that faunal provinces within the Campanian are an artefact of sampling: most Campanian dinosaur occurrences are known from Alberta, Montana, southern Utah, and northern New Mexico, with far less sampling having occurred in northern Utah and Wyoming [16]. Three quantitative studies have investigated the provincialism hypothesis in dinosaurs of the Late Cretaceous Western Interior. Gates et al. [35] used a variety of statistical techniques to assess the similarity between Campanian northern and southern faunas, and found evidence for either two distinct provinces with a broad area of overlap between them, or a latitudinal diversity gradient. The statistical techniques employed were unable to distinguish between these two hypotheses, and their results regarding dinosaurs were inconclusive. They suggested a further investigation into the causes of dinosaur distribution in the Western Interior. Berry [36] used a phylogenetic approach to assess biogeography within the Campanian and found no evidence for endemic sub-clades of chasmosaurine ceratopsids, arguing that this would be expected if there was a major barrier to dispersal, or niche conservatism related to climate. Vavrek & Larsson [15] investigated faunal endemism in the Maastrichtian of Laramidia using measures of beta diversity. After correcting for sampling biases, they found little evidence of provincialism, instead suggesting a homogeneous dinosaurian fauna across the Western Interior region at the very end of the Cretaceous; however, they did not test to see whether apparent biogeographic patterns within the Campanian were also caused by sampling. Herein, we quantitatively test hypotheses of faunal endemism in both the Campanian and Maastrichtian using biogeographic and multivariate statistical approaches. We focus our study on ceratopsid and hadrosaurian dinosaurs, as these megaherbivores have well-understood phylogenies and taxonomies and have been at the centre of previous discussions of faunal provinciality in this region. The distinctiveness of northern and southern Laramidian provinces are tested using phylogenetically corrected Biogeographic Connectedness (pBC). This quantitative method uses a network approach to assess phylogenetic distances between taxa in different geographic areas, resulting in a metric that quantifies the degree of faunal provinciality versus cosmopolitanism. It has been used successfully to understand changes in faunal compositions through the Carboniferous–Permian transition [41], and the Permian–Triassic and Triassic–Jurassic extinction events [42]. We also introduce additions to the methodology that address concerns regarding variation in sampling through time. To investigate the impact of sampling bias on our results, we examine correlations between occurrences (records of specimens) and collections (sites where specimens have been collected) with latitude, and use multivariate regression to examine which sources of sampling bias best explain sampling patterns. We use our results to determine whether it is possible to identify true geographic patterns of biodiversity on a continental scale in this very well-sampled area. Since no complete phylogenetic analysis of all ceratopsids is available, we built an informal supertree of all ceratopsid taxa considered valid in recent phylogenetic analyses by combining the phylogenetic results of [18] for chasmosaurines and [19] for centrosaurines. We resolved polytomies by removing Nedoceratops, a taxon some workers consider to be invalid ([43], but see [44]), from the data matrix in [18] and re-analysing the dataset. This resolved polytomies in the clade containing the common ancestor of Eotriceratops, Triceratops, and all of its descendants. The resulting supertree includes 67 taxa and represents a consensus of current views on ceratopsian phylogeny (electronic supplementary material, figure S1a). The structure of the hadrosaurid tree is based on several key recent analyses [45–47]. We resolved polytomies and added taxa considered valid but not included in those references using other recent phylogenetic analyses [48–51]. The resulting supertree includes 55 taxa and represents a current reasonable estimate of hadrosaur phylogeny (electronic supplementary material, figure S1b). Age for North American hadrosaur and ceratopsid species was obtained from the primary literature. The formations in which taxa occurred were found from the Paleobiology Database (PBDB; www.paleobiodb.org), and the most recent absolute age estimate of those formations was obtained from the primary literature (see electronic supplementary material, OSM, for sources). The age and geographic data for taxa outside of North America were obtained from the PBDB. pBC requires a priori assignment of geographic regions to test hypotheses of biogeographic connectedness, so we assigned dinosaurs to either northern Laramidia or southern Laramidia. Northern Laramidia includes taxa found in Wyoming and further north; southern Laramidia includes taxa found in Utah and further south, following previous studies. Age data were used to time-calibrate the phylogenetic trees using the ‘timePaleoPhy’ function of the Strap package [52] in R v. 3.5.2 [53] with the minimum branch length option specified (type=’mbl’). While it would be ideal to use high-resolution bins to test patterns of biogeography through the Late Cretaceous [17], too few taxa would be present in each bin to permit the use of pBC. Consequently, we divided taxa into Campanian and Maastrichtian time bins, which also has the benefit of allowing for comparison between previous studies of faunal provincialism in this area. Where a taxon's stratigraphic range/uncertainty crossed the Campanian–Maastrichtian boundary, it was included in both time intervals. We calculated pBC for Campanian ceratopsids, Maastrichtian ceratopsids, and Campanian hadrosaurs. The sampling of Maastrichtian hadrosaurs was too sparse, particularly in southern Laramidia, to calculate meaningful pBC values. Trees were pruned to exclude taxa from timeslices other than the one being analysed, and were made ultrametric prior to analysis. pBC was calculated using the function BC of the package ‘dispeRse’ (available at github.com/laurasoul/disperse). We initially varied the constant µ (see [41]) from 1 to 15 million years; subsequent analyses used a constant µ of 10. Data were jack-knifed 1000 times to produce a distribution of possible pBC values. To address concerns about the potential for a relationship between pBC and taxon sample size [54], we calculated rarefaction curves for pBC for the ceratopsian data (to facilitate comparisons between the Campanian and Maastrichtian). Sample sizes were rarefied down to a minimum number of five taxa. Ninety-five per cent confidence intervals for the rarefaction curves were generated using 1000 replicates at each sampling level. In order to determine whether pBC for each clade and time interval was significantly different from random, we randomly permuted the geographic areas in which taxa are found. We generated 1000 permutations of the data for each clade and time interval and calculated pBC for each permutation. The pBC for the unpermuted data was compared to the distribution of permuted pBC values to establish statistical significance (p < 0.05). To investigate whether biogeographic patterns we observed in the pBC results were influenced by sampling bias, we downloaded raw occurrence data for ceratopsids and hadrosaurs for the Campanian and Maastrichtian from the PBDB. We then downloaded North American dinosaur-bearing collections and North American tetrapod-bearing collections for each timeslice, and plotted occurrences and collections with latitude. We compared the curves using Spearman's rho and Kendall's tau. To investigate the possible causes of sampling bias we identified, we statistically examined correlations between occurrences and outcrop area, depositional environment, and proxies for exposure. First, we imported publicly available United States Geological Survey (USGS) state-level and Canadian Province digital geological maps (www.ngmbd.usgs.gov; https://ags.aer.ca/publication/map-600; https://geohub.saskatchewan.ca/datasets/bedrock-geology) into ArcMap 10 (www.esri.com), identified Campanian and Maastrichtian strata, and assigned an environmental attribute determining whether strata were deposited in a terrestrial, marine, or mixed setting (OSM). These data, along with maximum green vegetation fraction (MGVF) and slope, both proxies for exposure, were imported into R (version 3.5.0). Methods for generating MGVF and slope are provided in OSM. Level plots of total outcrop area, terrestrial, mixed and marine outcrop area, slope, and MGVF were produced using the levelplot() function of the rasterVis() package [55] (electronic supplementary material, figure S2). To investigate the power of each or a combination of these variables to explain the dinosaur occurrence data, we carried out a model-fitting approach using generalized least-squares regression (GLS). Ceratopsian and hadrosaur occurrences from the PBDB were counted in each 1-degree latitudinal bin (latitude is modern latitude, rather than palaeolatitude). Models compared latitudinal changes in ceratopsian and hadrosaur occurrences to changes in four different measures of outcrop area (see OSM), MGVF, and slope. GLS autoregressive models were fitted to combinations of potential explanatory variables. We used a first-order autoregressive model (corARMA) fitted to the data to account for spatial autocorrelation using the function gls() in the R package nlme v. 3.1–150 [56]. GLS reduces the chance of overestimating the statistical significance of regression lines due to serial correlation in the latitudinal series. Data series were ln-transformed prior to analysis to ensure normality and homoskedasticity of residuals. We calculated likelihood-ratio-based pseudo-R2 values using the function r.squaredLR() of the R package MuMIn [57]. Results were compared using Akaike's information criterion for small sample sizes (AICc) and Akaike weights were calculated to identify the best combination of explanatory variables from those tested. AICc was calculated using the function AICc() of the R package qpcR [58], and Akaike weights calculated using the aic.w() function of the R package phytools [59]. To test the impact of the Campanian bimodal sampling distribution on pBC results, we ran a second pBC test where we randomly removed 95% of ceratopsian taxa from the Maastrichtian that occurred between 35 and 50 degrees of latitude. We chose these latitudinal boundaries to enforce a similar bimodal latitudinal diversity gradient on the Maastrichtian data as seen in the Campanian (see Results). The remaining distribution of occurrences was used to re-run pBC analyses (with a µ of 10), and this process was repeated 1000 times for increased accuracy of results. pBC scores were recorded for each run, and the resulting distribution was used to calculate the mean pBC to compare against the original Maastrichtian ceratopsian pBC score and produce a probability density curve to estimate the probability of different values of pBC scores. The observed value of pBC for Campanian ceratopsids was 0.05, while that for Campanian hadrosaurs was 0.11, and for Maastrichtian ceratopsids the observed value was 0.16. pBC was therefore lower for ceratopsids in the Campanian than in the Maastrichtian, and endemism was correspondingly higher, in agreement with previous studies [33,36]. Jack-knifed distributions of ceratopsid pBC for the Campanian and Maastrichtian overlap (figure 1a–c), but their median values are strongly significantly different from each other (Wilcox Test, W = 60235, p = 0.00). Rarefaction curves for ceratopsids for the Campanian and Maastrichtian indicate a much higher pBC in the Maastrichtian than in the Campanian at equivalent levels of sampling, although the confidence intervals do overlap, particularly at lower sampling levels (figure 1d). This demonstrates that the higher pBC of the Maastrichtian is not a consequence of sampling lower numbers of species in that interval in comparison to the Campanian. Higher pBC equates to more cosmopolitan faunas, and thus this result supports lower endemism in Laramidia during the Maastrichtian when compared to the Campanian. Figure 1. (a–c) Jack-knifed distributions of pBC values for Campanian (green) and Maastrichtian (purple) ceratopsids. (d) Rarefaction curves for Campanian (green circles) and Maastrichtian (purple triangles) ceratopsids. Error bars show the 95% confidence intervals of values obtained during rarefaction. (Online version in colour.)
Values for pBC for both Campanian and Maastrichtian ceratopsid data are statistically significantly lower than for datasets in which the geographic areas are randomized (Campanian, p = 0.00; Maastrichtian, p = 0.015; electronic supplementary material, figure S3), and the same is true for the Campanian hadrosaur data (p = 0.00; electronic supplementary material, figure S3). This indicates that endemism was statistically significantly higher than in all randomized datasets across both time intervals and supports previous qualitative hypotheses of distinct northern and southern provinces in Laramidia (e.g. [33–35]) Curves of raw occurrence data with latitude for hadrosaurs and ceratopsids in both the Campanian and the Maastrichtian correlate strongly and statistically significantly with both dinosaur-bearing and tetrapod-bearing collections (figure 2; electronic supplementary material, figure S4; OSM). During the Campanian, sampling and occurrences are focused at two latitudes: 51–49 degrees north, which corresponds with the Dinosaur Park, Oldman and, to a lesser extent, the Foremost formations, and 36–37 degrees north, which corresponds primarily with the Kirtland/Fruitland, Aguja, and Kaiparowits Formations (figure 2a,b; electronic supplementary material, figure S4a,b). These two areas have been sampled orders of magnitude better than the surrounding latitudinal bins [16], although there are tetrapod- and dinosaur-bearing formations across the majority of the Western Interior at this time (figure 2a,b; electronic supplementary material, figure S4a,b). In the Maastrichtian, sampling is more evenly spread across the range of latitudes for which we have hadrosaur and ceratopsid body fossils (figure 2c,d; electronic supplementary material, figure S4c,d; [16]). These data are strongly indicative that the provinciality observed based on raw data in the Campanian could be due to intensive sampling in the Dinosaur Park Formation and Kirtland/Fruitland Formations with a lack of sufficient sampling between, and our observed increase in pBC (= reduced endemism) in the Maastrichtian is due to increased latitudinal coverage of sampling. Figure 2. (a,c) Ceratopsid occurrences and dinosaur-bearing collections with latitude in the (a), Campanian and (c), Maastrichtian. (b,d) Hadrosaurid occurrences and dinosaur-bearing collections with latitude in the (b), Campanian and (d), Maastrichtian. τ = Kendall's tau; ρ = Spearman's rho; DBCs, dinosaur-bearing collections. (Online version in colour.)
The mean pBC score of Maastrichtian ceratopsians subjected to a Campanian-style sampling distribution was 0.0351 with a standard deviation of 0.0476, significantly lower than the original pBC score of 0.16. The probability of a pBC score less than or equal to 0.8 was 0.78 (OSM and electronic supplementary material, figure S5). These results provide a further indication that sampling bias is driving pBC scores of Campanian fauna. A lack of sampling in the area between 49 degrees north and 37 degrees north (the ‘sampling peaks’) in Campanian strata could be caused by a variety of factors. It has long been known that rock outcrop area is strongly correlated with raw diversity (e.g. [60,61]); if there is less outcrop, there are fewer opportunities for palaeontologists to sample the rocks, and fewer fossils found as a consequence. As terrestrial organisms, the vast majority of dinosaur fossils are found in formations that were deposited on land. If Campanian rocks between the sampling peaks are primarily marine, there will be fewer opportunities for dinosaur fossils to be preserved, and thus fewer opportunities for them to be sampled by palaeontologists. Fossils are primarily found where bare rock is exposed at the surface. If less rock is exposed between the sampling peaks than in the areas of the peaks themselves, there will be fewer opportunities for fossils to come to light. GLS analyses recovered the following best models (highest AICc weights) for outcrop and tetrapod occurrence masks (see OSM for additional results): Campanian hadrosaurs, summed outcrop area + MGVF + slope; Campanian ceratopsians, non-marine total outcrop area; Maastrichtian hadrosaurs, null model; Maastrichtian ceratopsians, null model. However, in nearly all cases the correlations are non-significant (OSM) and only the Campanian hadrosaur model results had a strong overall explanatory power (OSM). This indicates that the potential sampling bias with latitude in the Campanian that we have identified cannot be fully explained by any of these variables and other sources of sampling bias that are hard to quantify may additionally be responsible. Several authors have suggested that the apparent faunal provincialism in Laramidia during the Late Cretaceous is an artefact, either because the formations in which dinosaurian taxa have been found are not contemporaneous [17,37,40] or due to uneven sampling of the fossil record [15,16]. Our results show that while the raw data clearly supports faunal endemicity, particularly in the Campanian, this pattern is driven by a lack of sampling outside of two specific latitudinal belts on Laramidia (51–49 degrees north, which corresponds with the Dinosaur Park, Oldman and Foremost formations, and 36–37 degrees north, which corresponds primarily with the Kirtland/Fruitland, Aguja, and Kaiparowits formations). This sampling bias cannot be fully explained by differences in outcrop area across the region, or by differences in slope or vegetation, which are factors that affect rock exposure. There are numerous other factors that can bias sampling, but these are very difficult to quantify. Low sampling between the northern and southern sampling peaks could occur if palaeontologists have yet to prospect the area to the same degree that they have in the north and south. The Late Cretaceous of the Western Interior has been intensively sampled for dinosaur fossils for over 100 years, and it is now probably the best-known Late Cretaceous ecosystem anywhere on Earth [15,16]. It is therefore highly unlikely that large parts of it remain unexplored for dinosaurs, and the fact that dinosaur fossils are known from the area between the sampling peaks during the Maastrichtian suggest the area has been explored. The lack of exploration for fossils is therefore unlikely to be the primary driver of the uneven sampling patterns we have observed. The ‘common cause hypothesis’ (e.g. [62]) suggests that correlations between raw diversity and sampling proxies (e.g. numbers of formations) are driven by a third factor, usually sea level. Although initially formulated for marine environments, the possibility of a sea-level driven common cause on land has also been discussed (e.g. [63]). During sea-level high stands sediment flux to inner shelves and marginal marine areas is high; this results in both high potential for the preservation of fossils due to rapid burial and high diversity due to habitat fragmentation leading to endemism and increased beta diversity. Conversely, sediment bypasses inner shelf environments during low stands, reducing sediment flux and leading to poorer preservation of fossils due to a lower chance of rapid burial, while diversity is lower due to cosmopolitanism. Although the effect of eustatic sea-level changes on the global terrestrial fossil record of vertebrates has been questioned [63], Chiarenza et al. [16] demonstrated that the areas of our northern and southern sampling peaks correlated with high sediment fluxes and low runoff rates during the Campanian. It is therefore possible that reduced sampling between our sampling peaks is because this area was less suitable for fossil preservation in the Campanian. Indeed, Chiarenza et al. [16] suggested that faunal provinicialism in the Campanian was a sampling bias at least partially due to variation in climatic-induced taphonomic suitability between northern and southern regions. Historical collecting practises and/or land ownership might also play a role in the sampling patterns we have observed. If the proportion of outcrop on public land was reduced in the areas outside of our sampling peaks, this might mean palaeontologists have less access to explore there for fossils. Furthermore, if there is a particularly field-active palaeontological institution close to an area of Campanian outcrop, or long-term agreements in place with landowners, this may have allowed prospecting to occur more regularly over a longer period of time in specific areas. A bias may also be introduced by uneven regional entry of data into online databases such as the PBDB. Such a bias could stem from monographs on specific formations or museums whose collections focus on specific areas that also have online databases. Data from these sources are comparatively easy to enter into the PBDB and thus could be contributing to the sampling patterns we observed. It seems highly likely that a combination of available outcrop area, rocks suitable for the fossilization of vertebrate remains, and an interplay between climate, topography, and historical collection and data entry practises is responsible for variations in sampling across the Western Interior, which have resulted in apparent northern and southern faunal provinces on Laramidia. Despite the fact that we find faunal provincialism in the Late Cretaceous Western Interior to mostly likely be due to sampling bias based on currently available data, it is clear that different taxa are found in the northern and southern areas of Laramidia [33–35]. This is especially clear in chasmosaurine ceratopsids, where there is almost no overlap at all between taxa found in the north and those found in the south [34], but see [36,39]. It has been demonstrated that many of these taxa were not contemporaneous [40], which would at least partially explain taxonomic differences. But, in addition, the study area covers 12 degrees of latitude and climate would have varied significantly over that area, even in a greenhouse world where latitudinal temperature gradients were reduced relative to today [39,64]. General circulation models for the Campanian show significant variation in mean annual temperature and rainfall patterns with latitude across Laramidia [16] and recent research has suggested elevated temperature gradients in a climatic transition zone between the northern and southern faunal provinces [38]. Given that there is evidence for both spatial [26] and functional [29] niche partitioning in Laramidia's large herbivores, taxonomic differences between the north and south could be related to climatic preference, and there may well have been a latitudinal biodiversity gradient across the area. Unfortunately, we have demonstrated here that that raw data is currently too influenced by sampling biases for such biodiversity patterns to be reconstructed. We show that data quality of Campanian and Maastrichtian ceratopsids and hadrosaurs, two of the most abundant clades of dinosaurs in the Late Cretaceous of North America, is currently too poor to enable fair tests of endemicity and provincialism. In order to effectively test hypotheses regarding the causative mechanisms of biodiversity distribution, palaeontologists must demonstrate either that the fossil record preserves true biodiversity patterns at high levels of temporal resolution, or that methods exist that can adequately overcome sampling biases. The Western Interior region represents probably the most densely sampled Late Cretaceous terrestrial region worldwide [15,16], but even in this intensively sampled area, it is not currently possible to reconstruct diversity patterns at the regional scale. In order for palaeontologists to make a meaningful contribution to ecological hypotheses about future biodiversity change, we must focus our efforts on smaller scale case studies, where temporal resolution is high, stratigraphic correlation is well established, and where sampling biases are likely to be more homogeneous and can be more easily quantified. A good example of a recent such study is [65]. The results of multiple high-resolution case studies can then be compared globally to establish the rules that governed past biodiversity distributions. Raw data, additional methods, results, and figures can be found in the electronic supplementary material. Further raw data, all code and a copy of the electronic supplementary material are available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.bcc2fqzbz [66]. S.C.R.M. and R.J.B. designed the study. S.C.R.M., R.J.B., and C.D.D. wrote the manuscript. R.I.M. and S.C.R.M. collected data and generated supertrees. S.C.R.M., R.J.B., and C.D.D. ran analyses. All authors provided critical comments on the manuscript. We declare we have no competing interests. R.J.B. was funded by European Union's Horizon 2020 research and innovation programme under grant agreement no. 637483 during the course of this work. Graeme Lloyd (University of Leeds) and David Button (Natural History Museum) wrote the original code to implement pBC. Caleb Brown (Royal Tyrrell Museum of Palaeontology) and an anonymous reviewer provided detailed and thoughtful comments that improved this manuscript. This is Paleobiology Database official publication 404. FootnotesElectronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5459154. © 2021 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 11Flowers are the most successful plant reproductive structures ever to evolve on land [1] and angiosperms (flowering plants) presently are the most abundant and diverse clade of vascular plants [2], currently consisting of over 369 000 described species [1]. This lineage probably originated very early during the Cretaceous, with robust molecular phylogenies placing the origins of the clade at 139.35–136 Ma [3]. This timing is consistent with the earliest documented appearance of angiosperm pollen around 136 Ma [4] and earliest known intact flowers at 125 Ma, from the Early Cretaceous of northeastern China [2,5]. The best documented and earliest known bisexual flower, the ‘Rose Creek Flower’, from the Early Cretaceous Dakota Formation of the United States, examined in this report, is approximately 103 million years (103 Ma) in age [6,7]. Despite the long and sporadic record of fossil flowers and given the abundance and diversity of Dakota Formation flowers, the time is propitious for examination of insect florivory (electronic supplementary material, text S1), which is the consumption of flowers prior to seed coat formation [8], and associated pollination in probably the earliest angiosperm deposit conducive to such an assessment. Such an evaluation could provide a better understanding of the role that mutualisms and antagonisms of angiosperms, effected by their insect pollinators, had on their joint diversification [9]. Here, three hypotheses are posed that we seek to test in our study of Dakota Formation flowers and inflorescences. The first hypothesis is ‘general features of insect florivory and related pollination, as measured by damage patterns on Dakota flowers, are very similar or the same as those made by modern florivores’. The second hypothesis is ‘major taxonomic groups of insect florivores and pollinators from the Dakota Formation are very similar or the same as those of today’. The third, more focused, hypothesis is ‘nectar robbing was present on Dakota flowers’. By addressing these three hypotheses, we place this study in a broader context of Cretaceous angiosperm pollination and provide a glimpse into early angiosperm florivory and associated pollination about 30 million years after the earliest appearance of angiosperms. All compression or impression specimens of flowers and other reproductive material (heretofore termed flower specimens) described in this report were collected from the Dakota Formation of midcontinental United States at eight localities, each bearing a flora of mid-Cretaceous (latest Albian to earliest Cenomanian) age [7], equivalent to approximately 104–97 Ma [10], depending on the locality (electronic supplementary material, figure S1 and table S1). Flower specimens were collected during the 1970s and 1980s by D. Dilcher and colleagues and are stored at the Florida Museum of Natural History, in Gainesville, Florida. The localities are distributed along a younger to older north–northeast to south–southwest transect [11] (electronic supplementary material, figure S1). The localities were deposited along the eastern coast of the North American Mid Cretaceous Seaway that extended from the Gulf of Mexico to the Canadian Arctic [12,13]. Lithological facies within the Dakota Formation are dominated by shales and sandstones of the upper Janssen Clay Member and subjacent claystones of the lower Terra Cotta Clay Member [11,13]. Dakota depositional environments are represented by brackish estuaries, freshwater swamps, low energy channels, floodplain ponds and ox-bow lakes (electronic supplementary material, table S1). These biotic environments consisted of humid forests and woodland [11,13]. The Rose Creek I locality, where most of the flowers were deposited, represented a mangal-like marsh similar to extant communities in southeast Asia [14]. Dakota Formation floras occur in a northeastern to southeastern trend at Courtland, Minnesota; Pleasant Dale, Nebraska (NE); Rose Creek I and II, NE; Braun's Ranch, Kansas (KS); Linnenberger Brothers' Ranch, KS; Acme Brick Quarry, KS; and Hoisington, KS [11] (electronic supplementary material, figure S1). Known vascular plant species diversity throughout the Dakota localities consists of 134 species. An earlier conservative estimate is 150–200 angiosperm species/morphotypes present in the Dakota Formation, as there is less than 25% species overlap between any two localities [6,7,11,15–18] (electronic supplementary material, table S1). Previous occasional descriptions of flowers, inflorescences, infructescences and fruits are known from several studies [7,19,20], yet affiliations of vegetative taxa with reproductive taxa remain largely unknown, although a few mostly wind-pollinated flowers, are associated with vegetative material that have been described [21]. Except for one locality, localities with higher abundance and diversity of foliage also appear to have a higher diversity and abundance of floral morphotypes (electronic supplementary material, figure S2 and table S1). (In this report, we use the term, floral, to refer to a flower, and not a bulk flora involving principally foliage.) The only described and best-known flower from the Dakota Formation is Dakotanthus cordiformis [7]. This flower has five, crescent-shaped, nectariferous pads that occur at the base of the gynoecium, each of which is aligned with a sepal. Dakotanthus, the most abundant morphotype in our dataset, is a member of the Rosidae 1 clade [7] and apparently very similar to a modern taxon with a lobed nectary disc. Other Dakota flower morphotypes show poor development or apparent absence of nectaries or nectary-like structures. However, leaf taxa occurring in the same localities as the unaffiliated Dakota flowers and infructescences have been assigned to extant families within Austrobaileyales, Chloranthales, Canellales, Magnoliales, Laurales and Rosidae 1 [15], which share a common pattern of fluid rewards for pollinating insects [22]. This pattern consists of: (i) staminoidal appendages (sterile stamens) that produce at their base glandular secretions of nectar-like fluids, mucilage, or ‘viscous substances’; (ii) nectariferous glands at the base or tips of fertile stamens; (iii) stigmas that secrete nectar-like substances, usually at their tips; (iv) nectar secreting, parenchymatous tissue on the adaxial surfaces of petals or sepals; and (v) large, substantive glands at the base of stamens that would qualify as true nectaries [22]. From these observations, it is highly likely that Dakota flower morphotypes produced nectar or other secretory, nectar-like fluids that attracted insect florivores and pollinators. The collection of florivory data is analogous to data for foliage or other vegetative organs and follows the same system of evaluating plant–insect associations [23], extensively used in fossil herbivory studies [24]. This system uses the functional-feeding group (FFG)–damage type (DT) system in which the overarching unit of herbivory is the FFG, examples of which are hole feeding, margin feeding, surface feeding and piercing-and-sucking for Dakota flower damage. Each FFG encompasses several or more DTs, which are the basic units of damage for fossil herbivory studies. A DT may be used in three ways. First, a DT may be used in terms of DT richness, referring to the kinds of DTs present; or as DT occurrences, as in the individual instances of damage of on a leaf; or as a formal name, such as DT405, which is a defined, specific mode of margin feeding damage. Details of photodocumentation and statistical methods are given in the electronic supplementary material, text S2. For Dakota plants, previous assessments of herbivory involved almost entirely mining damage on leaves [25–27]. However, Dakota plants, similar to amber deposits [28,29], provide considerable indirect evidence for flower–insect associations in the fossil record. Insect visitors to flowers are of two fundamental groups, florivores and pollinators [30]. Not all florivores are pollinators and not all pollinators are florivores, and the relationships between these two ecological guilds are complex [8]. Florivores typically leave damage on flowers, overwhelmingly on petals [31], often resulting in negative interactions [32]. However, some florivore interactions are neutral or even positive [31,33], as petals occasionally contain nutritive or highly scented tissues designed for consumption by florivores as pollinators [34,35]. Florivory can be a form of predation if plant embryonic tissues are destroyed before the opening of the flower, or if there is the consumption of immature pollen, features that do not appear present in bowl-shaped Dakota flowers, as the damage is overwhelmingly on inner petal surfaces. Consequently, florivores such as Orthoptera (katydids), Hemiptera (aphids, bugs), Thysanoptera (thrips), Coleoptera (beetles) and Hymenoptera (sawflies, wasps, bees) with mandibulate, stylate or similarly modified mouthparts [36], provide good proxy data for the broad spectrum of pollinator interactions on flowers [37] (table 1). However, a substantial component such as most adult Diptera [51] and Lepidoptera are nondamagers, as they do not leave damage on flowers.
The eight localities of the Dakota Formation consist of approximately 7858 total plant specimens, which yielded 645 (8.2%) flower specimens that were assessed for insect damage, some of which were photographically documented (electronic supplementary material, text S2). Flower specimens previously were identified to morphotype by Dilcher and Manchester [10,36], and Xiao, but mostly by the latter. The plant specimens were categorized into 32 flower morphotypes, one of which was Dakotanthus cordiformis (figure 1; electronic supplementary material, figures S2, S3, S6 and appendix 1). The 32 morphotypes consisted of Dakotanthus cordiformis [6,7], 14 flower morphotypes, eight inflorescence–infructescence morphotypes, five reproductive morphotypes, two flower–seed–fruit morphotypes and two Braun Ranch flower morphotypes. Unidentifiable specimens and poor preserved morphotypes, not assigned to one of the 32 morphotypes, were Acme unidentified inflorescence–infructescence, unidentified flower and unidentified stamen, which amounted to one, eight and nine specimens, respectively, attributable to a very limited local sample size or poor preservation. Based on the diversity and abundance of floral morphotypes (electronic supplementary material, figures S2 and S3), our estimate of the log normal fit is 182 species. We also obtained a Fisher's α value of 7.94 (electronic supplementary material, appendix S2). Figure 1. The florivore assemblage on Dakotanthus cordiformis (a–l), displaying petal DTs from all four FFGs of hole feeding, margin feeding, surface feeding and piercing-and-sucking. Specimen UF-12941 at (a) shows DT13 margin feeding cusps on the distal edge of the petal, enlarged in (b); smaller versions of DT12 margin feeding and DT46 punctures enlarged in (c). From the same specimen are several, V-shaped, margin feeding, DT405 notches along the petal edge, enlarged in (d), and further enlarged, including a DT46 ovate puncture, at (e). Specimen UF-3522 with several stamens at (f) shows a series of DT405 edge notches, enlarged in (g). Specimen UF-5612 displays DT12 margin feeding at the bottom of the petal at (h), enlarged at (i). Dakotanthus cordiformis (UF-5773) at (j) are DT02 hole-feeding damage, enlarged in (k), and, together with small dark, ovate punctures of DT46, enlarged in (l). Note well-developed reaction rim surrounding DT01. Scale bars: white, 5 mm; black, 1 mm; empty, 0.5 mm. (Online version in colour.)
Insect damage was present on 109 of the 645 examined flower and related specimens, for a specimen-based florivory rate of 17.2%. This damage was represented by four FFGs, 11 DTs [23] and 207 individual DT occurrences (electronic supplementary material, table S2). Some DTs occurred multiple times on the same specimen. One or more DTs on a specimen was present on nine of the 32 flower morphotypes (28.1%) (electronic supplementary material, table S2). The four FFGs present were hole feeding, margin feeding, surface feeding and piercing-and-sucking (see the electronic supplementary material, tables S2, S3 and text S3 for additional DT occurrence details). The distribution of DTs on plant hosts revealed three levels of host specificity [23]. Borrowing from studies of fossil herbivory as an example [24], host specificity is categorized as specialized damage if three or more occurrences of the same DT are present on the same host morphotype or on a very closely related host; damage is of intermediate specificity if the distribution of three or more occurrences of the same DT are present on more distantly related hosts; and generalized damage if three or more occurrences of the same DT are present on unrelated hosts [23]. For Dakota folivory data, because the phylogenetic relationships among flower morphotypes are unknown, terms expressing host specificity are referenced to the distribution of DTs on the flower species and morphotypes (electronic supplementary material, table S2). The three examples of specialist damage are small hole-feeding DT01 on flower morphotype 4 that hosts 12 of 13 (92.3%) of all occurrences; circular holes between 1 and 5 mm in diameter on Dakotanthus that hosts all eight (100%) of occurrences; and notched margin feeding of DT405 along the petal edges on Dakotanthus, which hosts 66 of 67 (98.5%) of all occurrences. The single example of damage of intermediate specificity is DT12, evidenced by too few DT distributions across three flower morphotypes. Seven examples of generalized damage are present. They are single, random, piercing-and-sucking damage assigned to DT46 on Dakotanthus and flower morphotypes 6, 9 and 10; clustered piercing-and-sucking assigned to DT402 on Dakotanthus and flower morphotypes 1, 5 and 8; and DT13, DT29, DT48, DT138 and DT383 that defaults to generalized specificity, each having only one or two occurrences on Dakotanthus and flower morphotypes 4, 5 and 7. This pattern of host specificity indicates three examples of specialized damage, one of intermediate specificity damage, and seven of generalized damage (electronic supplementary material, table S2). Of the eight localities examined, flower morphotypes from three localities—Rose Creek I, Rose Creek II and Braun's Ranch—showed evidence of florivory. The combined Rose Creek I and II localities exhibited three flower morphotypes with 117 DT occurrences, whereas the Braun's Ranch locality showed a higher diversity of seven flower morphotypes and 90 DT occurrences (electronic supplementary material, table S2). The number of florivorized to total flower morphotypes at each locality (electronic supplementary material, table S2)—three of 10 (30%) at Rose Creek I and II, and seven of 13 (53.8%) at Braun's Ranch—are distinctly significant subsets of the number of available hosts at each locality. At Rose Creek, Dakotanthus overwhelmingly was the dominant flower morphotype present, which displayed a rich spectrum of damage, with three of the four FFGs and seven of the 11 DTs represented (electronic supplementary material, table S2). By comparison, Braun's Ranch showed two florivorized flower morphotypes with a less rich spectrum of FFGs and DTs. Flower morphotype 4 had three of four FFGs and five of 11 DTs present. Similarly, flower morphotype 5 displayed three of four FFGs and six of 11 DTs. These latter flower morphotypes from the Rose Creek and Braun's Ranch localities exhibited a similar distribution of FFGs and DTs. Hole feeding on flowers of the Dakota Formation mostly is single, small and circular perforations of the entire petal thickness that are ovate or circular in shape (electronic supplementary material, text S4). Reaction rims are variably developed and occasionally associated with necroses of adjacent petal tissue. Damage type DT01 consists of holes 1 mm or less in diameter and is associated with flower morphotype 4 (electronic supplementary material, figure S5H–J). DT02 consists of holes between 1 and 5 mm in diameter and occurs on Dakotanthus (figure 1j–k). Dakota hole feeding consists of 21 perforations of DT01 and DT02 that represent 10.1% of all DT occurrences among flower morphotypes. Of all hole-feeding occurrences, 38.1% was present on Dakotanthus cordiformis from the Rose Creek I and II localities and 57.1% on flower morphotype 4 from the Braun's Ranch locality. It was noted that 71.4% of hole feeding was present on the lower half, rather than the upper half of the petals. Margin feeding on Dakota Formation flowers is represented by DTs DT12, DT13 and DT405. DT12 and DT13 consist of cuspate to U-shaped excisions, typically several mm in chord length, occurring along the edges of petals and sepals (electronic supplementary material, text S4). The cut edge, in addition to a bordering rim of dark reaction tissue, occasionally displays micromorphological features such as protruding veinal stringers, necrotic tissue flaps and cuspules within the overall cut edge, analogous to damage on foliage. DT12 occurred along the petal side edges on Dakotanthus (figure 1c,h,i), flower morphotype 4 (electronic supplementary material, figure S5E,F), flower morphotype 5 (electronic supplementary material, figure S4A) and flower morphotype 17 (electronic supplementary material, figure S6F). DT13 was present on the tips of petals of Dakotanthus (figure 1a,b) and flower morphotype 5 (electronic supplementary material, figure S4A). Careful examination of DT13 was required to determine if the damage was present, to eliminate confusion with a retuse or apically embayed margin. DT405 is a newly described DT (electronic supplementary material, text S5) and previously has not been recorded in the fossil record. The only example of a surface feeding FFG on a Dakota flower morphotype was DT DT29 (not illustrated) occurring on flower morphotype 7. DT29 is highly variable in size and shape, featuring polylobate to ovate patches of surface-fed petal tissue with distinct development of a reaction rim resulting from abrasion, scraping or delamination of a surface tissue layer (electronic supplementary material, text S4). Piercing-and-sucking damage of Dakota floral morphotypes is represented by the five DTs of DT46, DT48, DT138, DT383 and DT402. These DTs consist of various patterns of punctures that penetrate or slice into shallow to deep floral tissues (electronic supplementary material, text S4). DT46 and DT48 are single, randomly dispersed punctures less than 1 mm in diameter, present on petals or other flower elements. DT46 is a circular, concave mark with a crater-like rim and occurs on Dakotanthus (figure 1c,e,j,l), flower morphotype 1 (electronic supplementary material, figure S5A,B), and flower morphotype 5 (electronic supplementary material, figure S5C,D,M–O). By contrast, rare DT48 (not illustrated) is an elliptical puncture, with either a cratered rim or a convex central boss. DT138 are linear rows of punctures that occur on flower morphotype 5 (electronic supplementary material, figure S4D,E). DT383 are compact circular to polylobate clusters of punctures that probably accessed deeper tissues (electronic supplementary material, table S5). DT402, a newly described DT (electronic supplementary material, text S5), represents typically elongate, compact clusters of punctures in shallow tissues that occur especially along petal or sepal edges. DT402 occurs on flower morphotype 4 (electronic supplementary material, figure S5F,G and S5K–M,O,P) and flower morphotype 5 (electronic supplementary material, figure S4B,C,F–I). No preferential occurrences of piercing-and-sucking DTs were noted for the five DTs occurring on the upper versus lower halves of the petals. Associations were established between the pattern of Dakota insect damage with the relevant insect mouthpart class borne by an insect that would have produced that damage [36,52,53]. Such relationships, based on modern data [36] (electronic supplementary material, table S3), were made to better constrain the identities of potential florivores and pollinators (electronic supplementary material, text S6). However, Dakota damage caused by adult ectognathate, larval ectognathate and sericterate mouthpart classes typically could not be separated from each other. These three mandibulate (chewing) mouthpart classes account for 37.7% of all DT occurrences. Damage caused by the maxillolabiate (a complex apparatus for nectar-extraction) and rhynchophorate (perforating) mouthpart classes that result in hole-feeding damage cannot be distinguished from each other, but these two mouthpart classes that create holes on leaves account for 10.1% of all DT occurrences. Damage caused by the segmented beak mouthpart class, responsible for puncturing deeper tissues with stylate mouthparts, involve piercing-and-sucking feeding and account for 31.4% of all DT occurrences. Damage attributable to mouthcone mouthparts modified for punch and sucking of shallow, epidermal tissues account for 20.8% of DT occurrences. These seven mouthpart classes, responsible for four major feeding styles, indicate that florivory was dominated by edge feeders and tissue penetrating piercer and suckers, and less so by hole feeders and shallowly penetrating piercer and suckers, reconstructed in figure 2. Figure 2. A reconstruction of the insect pollinator community on Dakotanthus cordiformis [7] based on patterns of florivory. This scene is from the Rose Creek locality of the Early Cretaceous (late Albian) Dakota Formation of Southwestern Nebraska, USA. Painted by Xiaoran Zuo. (Online version in colour.)
For the five most prevalent flower morphotypes, the two metrics expressing the abundance (DT occurrences) and percentage of DTs present are a near-exact match of each other (electronic supplementary material, table S4). This near duplication is shown in (i) the percentage contribution of each morphotype to the total, (ii) the morphotype rank order, and (iii) the cumulative totals. Much less similar is the third metric of the percentage of specimens that are florivorized, which departs from the two other metrics in exhibiting greater differences in the contribution of each morphotype to the total, a different morphotype rank order after the first two most florivorized morphotypes of Dakotanthus and flower morphotype 4, and a higher cumulative total of 97.7 versus the first two of 86.6 and 85.3. A plot of the percentage of florivory of the three dominant FFGs of hole feeding, margin feeding and piercing-and-sucking was made for the major florivorized flower morphotypes of Dakotanthus, flower morphotype 4 and flower morphotype 5 for the three Dakota localities of Rose Creek I, Rose Creek II and Braun’ Ranch (figures 3 and 4). This analysis revealed three patterns. First, the florivore–pollinator communities of flower morphotypes 4 and 5 at Braun's Ranch locality are dominated by the piercing-and-sucking FFG, whereas the hole feeding and margin feeding FFGs played a minor role at both localities. By contrast, the Dakotanthus florivore faunas are dominated by the margin feeding FFG at the Rose Creek I and II localities, whereas the hole feeding and piercing-and-sucking FFGs have minor roles at both localities, although piercing-and-sucking appears to be subdominant at the Rose Creek II locality. Second, whereas margin feeding or piercing-and-sucking may have played a dominant role, depending on locality, hole feeding always had a minor role in the florivory spectrum across all localities. Third, with the exception of hole feeding, the florivore communities at Braun's Ranch versus the Rose Creek I and II localities, were largely feeding inversions of each other, suggesting heterogeneity in the pollinator assemblages by locality. Figure 3. Percentage representation of folivory for the five most insect-damaged flower morphotypes.
Figure 4. Per cent representation of florivore FFGs for flower morphotypes and sites. HF, hole feeding; MF, margin feeding; SF, surface feeding; PS, piercing-and-sucking. (Online version in colour.)
The presence of four FFGs subsuming 11 DTs in Dakota Formation plants (figures 1 and 2; electronic supplementary material, figures S4–S6) provides a very modern cast to the documented florivory (electronic supplementary material, table S3, text S4 and figure S7). In terms of the distinctiveness of the insect-mediated damage, the mouthpart classes responsible for the damage (electronic supplementary material, text S6), the proportional distribution of the damage on the flower morphotypes (electronic supplementary material, table S4), and their near-identical comparison to damage produced by known modern taxa [8] affirms the hypothesis that Dakota florivory is indistinguishable from its modern equivalent. This similarity suggests a similar pollinator community, although evidence for fluid feeding, nondamaging adult taxa is circumstantial. The richness of Dakota florivory and suggested associated pollinators provides, to our knowledge, the first extensive evidence for a community of insect visitors on the earliest, well-documented bowl-shaped flowers (electronic supplementary material, table S5, figure S7 and text S7). The florivore component of Dakota insect pollinators is established based on (i) distinctive DTs, (ii) extant lineages that were present during Dakota time from fossil occurrences, and (iii) evidence of relevant fossil occurrences or presence of closely related clades (phylogenetic bracketing) (electronic supplementary material, texts S6 and S7). The florivore assemblage consisted of a variety of insects with mandibulate and piercing-and-sucking mouthparts that produced recognizable damage patterns on flowers (electronic supplementary material, text S6). This pattern affirms the hypothesis that major taxonomic groups of insect florivores were very similar to their extant counterparts (electronic supplementary material, table S6 and text S8). Probably associated with this community of insect damagers of floral tissue, the indirect evidence indicates presence of lineages of fluid feeding, adult and nondamager taxa that left no trace on Dakota flowers. Based on the Dakota florivory data and modern studies, the core pollinators were heteropterans (especially pentatomorphs), thrips, polyphagan beetles (principally scarab and leaf beetles, and weevils), and bees. A subordinate component of early-diverging moths, sawflies, several major lineages of nematoceran and brachyceran flies, and perhaps parasitoid wasps probably were present (electronic supplementary material, table S6 and text S7). The data on hole feeding across the flower morphotype hosts is intriguing. Hole feeding consisting of 23 perforations of DTs DT01 and DT02 represent 10.1% of all DT occurrences among Dakota flower morphotypes. Of hole-feeding occurrences, 95.2% were present on Dakotanthus from the Rose Creek I locality and on flower morphotype 4 from the Braun's Ranch locality. There is a distinct preference for the lower half (71.4%) rather than the upper half of the petal. All holes were circular or nearly so, and no evidence of slits or tears was observed as holes on the petals. Although these data suggest that nectar robbing was present, nectar robbing would be ineffectual if flowers are open, bowl shaped, and with rewards such as nectar readily available to insect visitors [54]. All modern nectar robbing occurs with tubular or similar flowers that that are highly enclosed and have access through a narrowly throated corolla [54]. The floral morphology of Dakota flowers, exemplified by Dakotanthus, is inconsistent with nectar robbing, and holes are circular and rounded in contrast with modern nectar robbers that construct slits or holes with jagged outlines [55,56]. Consequently, hypothesis 3, that nectar robbing was present on Dakota flowers, is rejected. Given this outcome, a more productive search for the earliest nectar robbing would be among early occurrences of tubular or otherwise enclosed flowers in the younger Late Cretaceous [2]. The early fossil history of angiosperm pollination [2] is illustrated by the bowl and similarly shaped floras in the 21 localities listed in geochronological order that provide data on inferred insect pollinator lineages and their functional-feeding group membership (electronic supplementary material, table S5 and text S7). This list shows that the Dakota insect pollinator fauna probably had taxonomic similarities to that of extant basal angiosperms (electronic supplementary material, table S6; texts S7 and S8). The next, four-million year-younger assemblage of pollinating insects originates from the very geographically, ecologically and taphonomically different locality of Myanmar amber, which shows a distinct pollinator fauna (electronic supplementary material, table S5 and figure S8). Subsequent, Late Cretaceous floras originate from a variety of localities that reveal the expansion of dicot angiosperms and an associated pollinator fauna. This study provides a new approach for the study of pollination in the fossil record. By examining exquisitely preserved insect damage on Dakota flowers, the three hypotheses have been tested that were initially proposed in this study. The first hypothesis—florivore damage patterns on Dakota flowers are similar to those of today—is supported. One proviso to this conclusion is that evidence for pollinating insects on Dakota flowers result from insect immatures and adults capable of leaving detectable damage, principally those with mandibulate and stylate mouthparts or their modifications, but excludes adult fluid feeding insects whose mouthparts are incapable of damaging tissues. Nevertheless, the indirect evidence of phylogenetic bracketing indicates the existence of these nondamaging, fluid feeding and pollinating clades. The second hypothesis—major taxonomic groups Dakota insect florivores are similar to those of today—is supported. Again, a caveat is that major insect taxa lacking the ability to inflict damage on flowers remain undetected. Circumstantial evidence based on fossil occurrences and phylogenetic bracketing indicates that typical nectar-feeding taxa were present. Last, the third hypothesis—nectar robbing was present—is rejected. The reason for rejection is not that a great preponderance of hole feeding was located at petal bases, but rather the open, bowl shape of Dakota flowers would obviate the need for nectar robbing. While addressing these three hypotheses extends understanding of florivory and pollination to the Dakota Formation, important gaps in knowledge of early angiosperm pollination remain. Highly relevant deposits from this time interval (electronic supplementary material, table 5, text S7 and figure S7) should be further explored. The data supporting the analyses of this article consist of (i) table 1 in this article; (ii) more extensive, linked electronic supplementary material that includes raw and related data among the three fossil localities as well as modern data; (iii) appendix 1 that provides detailed descriptions of flower morphotypes/species; and (iv) appendix 2 that provides R code files and Excel data files in a zip file. The data are provided in the electronic supplementary material. L.X.: conceptualization, data curation, formal analysis, investigation, methodology, supervision, validation, visualization, writing—original draft, writing—review and editing; C.L.: conceptualization, funding acquisition, methodology, project administration, supervision, validation, visualization, writing—original draft, writing—review and editing; D.D.: conceptualization, data curation, investigation, resources, supervision, validation, visualization, writing—original draft, writing—review and editing; D.R.: conceptualization, data curation, funding acquisition, project administration, resources, supervision, validation, visualization, writing—original draft, writing—review and editing. All authors gave final approval for publication and agreed to be held accountable for the work performed therein. We declare that we have no competing interests. This research was supported by the National Natural Science Foundation of China (grant nos. 31730087, 32020103006 and 41688103). Fieldwork was supported in part by National Science Foundation of United States grant nos. DEB 75-02268, DEB 75-19849, DEB 77-04846, DEB 10720 and EAR 79-00898 to D.D. We especially thank two reviewers who give us important comments. We thank Gussie Maccracken, Terry Lott, Finnegan Marsh, Steven Manchester, Hongshan Wang, Nareerat Boonchai and Xiaodan Lin who individually provided support in the Labandeira and Manchester laboratories. Sandra Schachat is thanked for producing figures 3 and 4. Gussie Maccracken provided a critique of a preliminary draft of the manuscript and offered needed overall advice. We are especially grateful to Steven Manchester for sponsoring the senior author at the Florida Museum of Natural History, in Gainesville for two months and their dedicated staff for general assistance. We thank Huayan Chen who provided photos of modern florivores and Torsten Dikow who identified the fly in electronic supplementary material, figure S6. We thank Liang Chen who helped us in figure construction. The Fossilworks database is acknowledged. This is contribution 378 to the Evolution of Terrestrial Ecosystems consortium at the National Museum of Natural History, in Washington, D.C. FootnotesElectronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.5448682. © 2021 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 12You have access Proc. R. Soc. B283, 20160477. (Published 25 May 2016). (doi:10.1098/rspb.2016.0477) (1) In the original paper by Start D, Gilbert B., “Host – parasitoid evolution in a metacommunity.”, we stated that we measured the maximum gall diameter of all collected galls. For this measurement, we used two methods that we then averaged. First, we used digital calipers that measure to the 100th of a millimetre. Second, we used a series of pre-drilled holes with known circumferences and calculated the diameter by dividing the circumference by π. The smallest hole that a gall fit through was considered its maximum diameter. In some cases, we destructively sampled galls for larvae before conducting the second measurement and therefore only report caliper data. The estimates in these cases have a lower precision (number of decimal places) than those that are from the averaged measures. (2) The data usage notes available from the Dryad Digital Repository refer to ‘cumulative distance moved by the fly’ and is given at specific time since release. These data should instead be described as ‘total distance from the origin’. Total distance from the origin at a given time point will therefore not always increase monotonically with time, if a fly reversed direction. This is the measure we used to parametrize our metapopulation model. Page 13Proc. R. Soc. B287, 20201671. (Published 21 October 2020). (doi:10.1098/rspb.2020.1671) (1) Under the section in the main text Methods: (d) Effect of intrusions of captive-bred fish on population productivity, the fourth line of this paragraph should read ‘This figure ranged from 0.01 to 0.61.’ (2) In the main text, Figure 1b should appear as such: Figure 1. (a) Overall and cohort-specific comparisons of RRS for captive- and wild-bred Atlantic salmon in the Burrishoole catchment, Ireland. Overall RRS comparison estimated as the weighted geometric mean of the six cohort point estimates. Significance of the overall comparison determined using FCPT, where χ2 = 117.94 with 12 degrees of freedom. Significance of cohort-specific comparisons was determined using one-tailed permutation tests. Horizontal line for emphasis of increase/decrease in reproductive success of captive-bred fish relative to wild-bred fish. Numbers on top of bars represent the number of captive-bred (left number) salmon and wild-bred (right number) salmon used in cohort-specific comparisons. *p < 0.05, **p < 0.01, ***p < 0.001. (b) Productivity of the mixed population as a function of the annual proportion of potentially spawning fish that were captive-bred. The solid line represents the line-of-best fit from a linear model, and shading represents the 95% confidence interval. (Online version in colour.)
(3) In the main text under the section Methods: (d) Effect of intrusions of captive-bred fish on population productivity, the first part of the sentence starting on the second line of the second column should be corrected as ‘While this density measure was only poorly correlated with the annual proportion of captive-bred fish (Pearson's correlation: r= 0.057, t= 0.37, d.f. = 41, p= 0.72)’ (4) In the main text under the section Results, the first sentence of the second paragraph should read ‘The population-level analysis revealed a significant negative relationship between our density-independent population productivity measure and the proportion of captive-bred fish in a spawning cohort (adjusted R2= 0.09, F1,41= 5.15, p-value= 0.0285).’ (5) In the main text under the section Results, the second sentence of the second paragraph should read ‘Population productivity at the mean value of the proportion captive-bred fish across the 43-year period (0.15) was reduced, on average, by 9.52% (back-transformed from the log scale), relative to a hypothetical pure population (proportion captive-bred fish=0)’. Page 14You have accessExpression of concern Proc. R. Soc. B285, 20180384. (Published online 25 July 2018). (doi:10.1098/rspb.2018.0384) Following the publication of ‘Predator macroevolution drives trophic cascades and ecosystem functioning’, by Start D, Proc. R. Soc. B285: 20180384 (published 25 July 2018), it has been brought to the attention of the Editorial team that there are concerns regarding this study. Questions have been raised with regard to missing or incomplete data, taxonomic identification of samples and the relationship between data published in this paper and another paper (Start et al. 2018 Physiology underlies the assembly of ecological communities. PNAS115 (23): 6016–6021). Proceedings B has approached the author and requested clarification on these points. Because the author has not been able to address the questions in a reasonable time, the journal is publishing this editorial Expression of Concern. Notification will be published if the concerns related to this paper are addressed. The journal welcomes the active participation of authors in the correction of any problems identified in the scientific record. Page 15
You have to enable JavaScript in your browser's settings in order to use the eReader. Or try downloading the content offline DOWNLOAD |
The evolution of the angiosperms lead to the emergence of a wide diversity of primates because
Pos Terkait
Copyright © 2024 ketajaman Inc.
