Philippa Kaur

Climate change drove evolution, biogeography, and demographyPhylogenetic results (Fig. 1 and Supplementary Fig. 2) confirm previous findings, recovering Aptenodytes (king and emperor penguins) as the sister clade to all other crown penguins, with brush-tailed (Pygoscelis) penguins in turn sister to two clades uniting the banded (Spheniscus) and little (Eudyptula) penguins and the yellow-eyed (Megadyptes) and crested (Eudyptes) penguins6,7,9. Biogeographical reconstructions (Fig. 1, Supplementary Figs. 3–4 and Supplementary Data 1) support a Zealandian origin for penguins6,7. Stem penguins radiated extensively in Zealandia before dispersing to South America and Antarctica multiple times, following the eastward-flowing direction of the Antarctic Circumpolar Current (ACC) (Fig. 1). Crown penguins most likely arose from descendant lineages in South America, before dispersing back to Zealandia at least three times. Interestingly, at least two such dispersals occurred before the inferred onset of the ACC system, suggesting that early stem penguins were not dependent on currents to disperse over long distances. A second pulse of speciation coincides with the onset of the ACC, though understanding whether this pattern is real or an artifact of fossil sampling requires more collecting from early Eocene localities. We infer an age of ~14 Ma for the origin of crown penguins, which is more recent than the ~24 Ma age recovered in genomic analyses, not including fossil taxa7 (Supplementary Fig. 2b) and coincides with the onset of global cooling during the middle Miocene climate transition4,10 (Supplementary Fig. 3a). This young age suggests that expansion of Antarctic ice sheets and the onset of dispersal vectors such as the Benguela Current11 during the middle to late Miocene facilitated crown penguin dispersal and speciation, as hinted at by fossil evidence12.Incongruences between species trees and gene trees were identified, e.g., alternate topologies occurred at high frequencies ( >10%) for several internal branches (Fig. 1c; Supplementary Fig. 5). These patterns indicate that gene tree discordance may be caused by incomplete lineage sorting (ILS) or introgression events. By quantifying ILS and introgression via branch lengths from over 10,000 gene trees, we found that the rapid speciation within crown penguins was accompanied by >5% ILS content within the ancestors of Spheniscus, Eudyptula, Eudyptes, and several subgroups within Eudyptes (Fig. 2a). Our dated tree provides a temporal framework for this rapid radiation: the four extant Spheniscus taxa are all inferred to have split from one another within the last ~3 Ma, and likewise the nine extant Eudyptes taxa likely split from one another in that same time (Fig. 1b). Many closely related penguin species/lineages are known to hybridize in the wild (see supplementary methods). Consistent with this, multiple analyses suggest that introgression also contributes to species tree—gene tree incongruence (Supplementary Figs. 6–9 and Supplementary Data 2; also see Supplementary Methods for further details). This could explain the most notable conflict in previous phylogenetic results, which showed inconsistency over whether Aptenodytes alone7 or Aptenodytes and Pygoscelis together4,5 represent the sister clade to all other extant penguins. Introgression was detected between the ancestor of Aptenodytes and the ancestor of other extant penguins, and is inferred to have occurred when the range of these ancestors overlapped in South America (Fig. 2a and Supplementary Data 2). Introgression ( >9%) was also detected between Eudyptula novaehollandiae and Eudyptula minor, and several introgression events were especially pervasive in Eudyptes (Fig. 2a and Supplementary Fig. 6).Fig. 2: Incomplete lineage sorting, introgression events, and demographic history among penguins.a Model of incomplete lineage sorting (ILS) and introgression events estimated from QuIBL and hybrid pairwise sequentially Markovian coalescent (hPSMC) results. hPSMC was only run for 20 species pairs (see b). Numbers on branches represent the proportion (%) of ILS (orange branches) or introgression (blue lines, blue dashed lines, and blue dotted lines) detected by QuIBL. Proportions 50 km; see Supplementary Data 117), while taxa that decreased towards the end of the LGP (e.g., S. humboldti, S. demersus, M. a. antipodes and likely M. a. richdalei) tend to be residential, and forage inshore; see Supplementary Data 1. Taxa that disperse farther may have overcome local impacts of global climate cooling during the LGP (e.g., changes in sea-ice extent, prey abundance and terrestrial glaciation, however see18) largely by relocating to lower latitudes (e.g.,14), whereas locally-restricted taxa may have been more prone to sudden population collapses.Penguins have the slowest evolutionary rates among birdsThe integrated evolutionary speed hypothesis (IESH) proposes that temperature, water availability, population size, and spatial heterogeneity influence evolutionary rate19. Life history traits also impact the evolutionary rate, but such relationships remain incompletely understood in birds20. Penguins are long-lived, large-bodied, and produce few offspring, thus providing an ideal case study in how life history may impact evolutionary rate. We tested the IESH using three proxies for evolutionary rate: substitution rate, P and K2P distances between lineages and their ancestors (Supplementary Fig. 12 and Supplementary Data 3). We found that penguins and their sister group (Procellariiformes) had the lowest evolutionary rates of the 17 avian orders sampled by21 (Fig. 3a, Supplementary Fig. 13, and Supplementary Data 3). Because other aquatic orders also show slow rates (e.g., the aquatic Anseriformes show a significantly slower rate than their terrestrial sister group Galliformes), we hypothesize that the rate in penguins represents the culmination of a gradual slowdown associated with increasingly aquatic ecology. Intriguingly, we detected a trend toward decreasing rate over the first ~10 Ma of crown penguin evolution, followed by a marked uptick ~2 Ma, which suggests the onset of glacial-interglacial cycles contributed to a recent increase in evolutionary rates in penguins (Fig. 3b).Fig. 3: Evolutionary rates in birds.a Evolutionary rate in avian orders based on a ~19 Mbp alignment of highly conserved genome regions. Sphenisciformes and Procellariiformes have the lowest evolutionary rate among modern bird orders (One-sided Wilcoxon Rank sum test, P values  0.1)). Numbers at the tips represent the sample size in each group. Numbers at nodes represent the divergence times (Ma) between each order and its sister taxon and red dots within the boxplots indicate average values. We did not attempt to estimate the evolutionary rates for orders containing less than three sampled species (gray font; Musophagiformes, Mesitornithiformes, and Struthioniformes). Boxplots show the median with hinges at the 25th and 75th percentile and whiskers extending 1.5 times the interquartile range. Some bird images were downloaded from phylopic.org and were licensed under the Creative Commons (CC0) 1.0 Universal Public Domain Dedication. b Evolutionary rates inferred for extant penguin lineages at internal nodes from the maximum clade credibility tree, calculated using a 500 Mbp genome alignment. Gray shadows represent the 95% credible intervals. c–e Correlations between c, body mass and generation time (P value  More

FAO. The State of World Fisheries and Aquaculture 2020., https://doi.org/10.4060/ca9229en (FAO, 2020).Watson, J. W. & Kerstetter, D. W. Pelagic Longline Fishing Gear: A Brief History and Review of Research Efforts to Improve Selectivity. Mar. Technol. Soc. J. 40, 6–11 (2006).Article 

Google Scholar 
Hare, S. R. et al. The western an our d central Pacific tuna fishery: 2019 overview and status of stocks. (SPC, 2020).Wang, S.-P. Stock assessment of blue marlin in the Indian Ocean using Stock Synthesis. (IOTC, 2019).Ohshimo, S. et al. Horizontal distribution and habitat of Pacific bluefin tuna, Thunnus orientalis, larvae in the waters around Japan. Bull. Mar. Sci. 93, 769–787 (2017).ADS 
Article 

Google Scholar 
Margulies, D., Scholey, V. P., Wexler, J. B. & Stein, M. S. Chapter 5 – Research on the Reproductive Biology and Early Life History of Yellowfin Tuna Thunnus albacares in Panama. In Advances in Tuna Aquaculture: From Hatchery to Market (eds. Benetti, D. D., Partridge, G. J. & Buentello, A.) 77–114, https://doi.org/10.1016/B978-0-12-411459-3.00004-7 (Academic Press, 2016).Madigan, D. J. et al. Intrinsic tracers reveal recent foraging ecology of giant Pacific bluefin tuna at their primary spawning grounds. Mar. Ecol. Prog. Ser. 553, 253–266 (2016).ADS 
CAS 
Article 

Google Scholar 
Ward, T. M., Staunton-Smith, J., Hoyle, S. & Halliday, I. A. Spawning patterns of four species of predominantly temperate pelagic fishes in the sub-tropical waters of southern Queensland. Estuar. Coast. Shelf Sci. 56, 1125–1140 (2003).ADS 
Article 

Google Scholar 
Kiyofuji, H. & Ochi, D. Proposal of alternative spatial structure for skipjack stock assessment in the WCPO. (WCPFC, 2016).Reglero, P., Tittensor, D., Álvarez-Berastegui, D., Aparicio-González, A. & Worm, B. Worldwide distributions of tuna larvae: revisiting hypotheses on environmental requirements for spawning habitats. Mar. Ecol. Prog. Ser. 501, 207–224 (2014).ADS 
Article 

Google Scholar 
Schilling, H. T. et al. Multiple spawning events promote increased larval dispersal of a predatory fish in a western boundary current. Fish. Oceanogr. 29, 309–323 (2020).Article 

Google Scholar 
Richardson, D. E., Hare, J. A., Overholtz, W. J. & Johnson, D. L. Development of long-term larval indices for Atlantic herring (Clupea harengus) on the northeast US continental shelf. ICES J. Mar. Sci. 67, 617–627 (2010).Article 

Google Scholar 
Muhling, B. A. et al. Overlap between Atlantic bluefin tuna spawning grounds and observed Deepwater Horizon surface oil in the northern Gulf of Mexico. Mar. Pollut. Bull. 64, 679–687 (2012).CAS 
Article 

Google Scholar 
Brown, S. K., Buja, K. R., Jury, S. H., Monaco, M. E. & Banner, A. Habitat Suitability Index Models for Eight Fish and Invertebrate Species in Casco and Sheepscot Bays, Maine. North Am. J. Fish. Manag. 20, 408–435 (2000).Article 

Google Scholar 
Hernández, C. M. et al. Evidence and patterns of tuna spawning inside a large no-take Marine Protected Area. Sci. Rep. 9, 1–11 (2019).Article 

Google Scholar 
Janßen, H. et al. Integration of fisheries into marine spatial planning: Quo vadis? Estuar. Coast. Shelf Sci. 201, 105–113 (2018).ADS 
Article 

Google Scholar 
Richardson, A. J. et al. Residency and reproductive status of yellowfin tuna in a proposed large-scale pelagic marine protected area. Aquat. Conserv. Mar. Freshw. Ecosyst. 28, 1308–1316 (2018).Article 

Google Scholar 
CBD. Protected areas and other effective area-based conservation measures (Decision 14/8). https://www.cbd.int/doc/decisions/cop-14/cop-14-dec-08-en.pdf (2018).Gurney, G. G. et al. Biodiversity needs every tool in the box: use OECMs. Nature 595, 646–649 (2021).ADS 
CAS 
Article 

Google Scholar 
Hiltz, E., Fuller, S. & Mitchell, J. Disko Fan Conservation Area: a Canadian case study. Parks 24, 17–30 (2018).Article 

Google Scholar 
IUCN-WCPA Task Force on OECMs. Recognising and reporting other effective area-based conservation measures. https://portals.iucn.org/library/sites/library/files/documents/PATRS-003-En.pdf (2019).Berkeley, S. A., Hixon, M. A., Larson, R. J. & Love, M. S. Fisheries Sustainability via Protection of Age Structure and Spatial Distribution of Fish Populations. Fisheries 29, 23–32 (2004).Article 

Google Scholar 
Hall, S. In A Fishery manager’s guidebook – Second Edition (eds. Garcia, S. M. & Cochrane, K. L.) 196–219 Ch.8 – Area and time restrictions. (Wiley-Blackwell, 2009).Jonas, H. D., Barbuto, V., Jonas, H. C., Kothari, A. & Nelson, F. New steps of change: looking beyond protected areas to consider other effective area-based conservation measures. Parks 20, 111–128 (2014).Article 

Google Scholar 
Dunn, D. C., Maxwell, S. M., Boustany, A. M. & Halpin, P. N. Dynamic ocean management increases the efficiency and efficacy of fisheries management. Proc. Natl. Acad. Sci. 113, 668–673 (2016).ADS 
CAS 
Article 

Google Scholar 
Le Quesne, W. J. F. & Codling, E. A. Managing mobile species with MPAs: the effects of mobility, larval dispersal, and fishing mortality on closure size. ICES J. Mar. Sci. 66, 122–131 (2009).Article 

Google Scholar 
Richardson, D. E. et al. Discovery of a spawning ground reveals diverse migration strategies in Atlantic bluefin tuna (Thunnus thynnus). Proc. Natl. Acad. Sci. 113, 3299–3304 (2016).ADS 
CAS 
Article 

Google Scholar 
Dueri, S., Bopp, L. & Maury, O. Projecting the impacts of climate change on skipjack tuna abundance and spatial distribution. Glob. Change Biol. 20, 742–753 (2014).ADS 
Article 

Google Scholar 
Asch, R. G. Climate change and decadal shifts in the phenology of larval fishes in the California Current ecosystem. Proc. Natl. Acad. Sci. 112, E4065–E4074 (2015).ADS 
CAS 
Article 

Google Scholar 
Mountain, D. G. & Murawski, S. A. Variation in the distribution of fish stocks on the northeast continental shelf in relation to their environment, 1980–1989. ICES mar. Sci. Symp. 195, 424–432.Muhling, B. A. et al. Potential impact of climate change on the Intra-Americas Sea: Part 2. Implications for Atlantic bluefin tuna and skipjack tuna adult and larval habitats. J. Mar. Syst. 148, 1–13 (2015).Article 

Google Scholar 
Asch, R. G., Stock, C. A. & Sarmiento, J. L. Climate change impacts on mismatches between phytoplankton blooms and fish spawning phenology. Glob. Change Biol. 25, 2544–2559 (2019).ADS 
Article 

Google Scholar 
Pankhurst, N. W. & Munday, P. L. Effects of climate change on fish reproduction and early life history stages. Mar. Freshw. Res. 62, 1015–1026 (2011).CAS 
Article 

Google Scholar 
Lehodey, P., Senina, I., Nicol, S. & Hampton, J. Modelling the impact of climate change on South Pacific albacore tuna. Deep Sea Res. Part II Top. Stud. Oceanogr. 113, 246–259 (2015).ADS 
CAS 
Article 

Google Scholar 
Muhling, B. A. et al. Collection of Larval Bluefin Tuna (Thunnus Thynnus) Outside Documented Western Atlantic Spawning Grounds. Bull. Mar. Sci. 87, 687–694 (2011).Article 

Google Scholar 
Nishikawa, Y., Honma, M., Ueyanagi, S. & Kikawa, S. Average Distribution of Larvae of Oceanic Species of Scombroid Fishes, 1956–1981. (Far Seas Fisheries Research Laboratory, 1985).Maunder, M. N. & Punt, A. E. Standardizing catch and effort data: a review of recent approaches. Fish. Res. 70, 141–159 (2004).Article 

Google Scholar 
Collette, B. B. In Annotated Checklist of Fishes Vol. 19 Family Scombridae Rafinesque 1815 – mackerels, tunas, and bonitos. (California Academy of Sciences, 2003).Fricke, R., Eschmeyer, W. N. & Van der Laan, R. Eschmeyer’s catalog of fishes; Genera, species, references. Electronic version. (California Academy of Sciences, 2022).R Core Team. R: A language and environment for statistical computing. (R Foundation for Statistical Computing, 2022).Buenafe, KCV. tinbuenafe/DigitizingNishikawa: Digitizing Nishikawa v3.0, Zenodo, https://doi.org/10.5281/zenodo.6592148 (2022).Hijmans, R. J. et al. terra: Spatial Data Analysis. (2022).Pebesma, E. Simple Features for R: Standardized Support for Spatial Vector Data. R J. 10, 439–446 (2018).Article 

Google Scholar 
Blondel, E. & Billet, N. RFigisGeo: A R package to handle utilities for geospatial processing. (2022).FAO. FAO Major Fishing Areas. https://www.fao.org/fishery/en/area/search (2022).Richards, W. J. & Potthoff, T. Analysis of the Taxonomic Characters of Young Scombrid Fishes, Genus Thunnus. In The Early life history of Fish (ed. Blaxter, J. H. S.) 623–648, https://doi.org/10.1007/978-3-642-65852-5_50 (Springer, 1974).Luthy, S. A., Cowen, R. K., Serafy, J. E. & McDowell, J. R. Toward identification of larval sailfish (Istiophorus platypterus), white marlin (Tetrapturus albidus), and blue marlin (Makaira nigricans) in the western North Atlantic Ocean. Fish. Bull. 103 (2004). More

Díaz, S. et al. The global spectrum of plant form and function. Nature 529, 167–171 (2016).ADS 
Article 

Google Scholar 
Céréghino, R. et al. Constraints on the functional trait space of aquatic invertebrates in bromeliads. Funct. Ecol. 32, 2435–2447 (2018).Article 

Google Scholar 
Winemiller, K. O., Fitzgerald, D. B., Bowler, L. & Pianka, E. R. Functional traits, convergent evolution, and periodic tables of niches. Ecol. Lett. 18, 737–751 (2015).Article 

Google Scholar 
Díaz, S. et al. Functional traits, the phylogeny of function, and ecosystem service vulnerability. Ecol. Evol. 3, 2958–2975 (2013).Article 

Google Scholar 
Leimar, O. Evolutionary change and Darwinian demons. Selection 2, 65–72 (2001).Article 

Google Scholar 
Cummins, K. W. & Klug, M. J. Feeding ecology of stream invertebrates. Annu. Rev. Ecol. Syst. 10, 147–172 (1979).Article 

Google Scholar 
Pianka, E. R., Vitt, L. J., Pelegrin, N., Fitzgerald, D. B. & Winemiller, K. O. Toward a periodic table of niches, or exploring the lizard niche hypervolume. Am. Nat. 190, 601–616 (2017).Article 

Google Scholar 
Rosenberg, D. M. & Resh, V. H. Freshwater Biomonitoring and Benthic Macroinvertebrates (Chapman and Hall, 1993).
Google Scholar 
Allan, J. D. & Castillo, M. M. Stream Ecology. Structure and Function of Running Waters 2nd edn. (Springer, 2007).
Google Scholar 
Wallace, J. B. & Webster, J. R. The role of macroinvertebrates in stream ecosystem function. Annu. Rev. Entomol. 41, 115–139 (1996).CAS 
Article 

Google Scholar 
Southwood, T. R. E. Habitat, the templet for ecological strategies?. J. Anim. Ecol. 46, 336–365 (1977).Article 

Google Scholar 
Townsend, C. R. & Hildrew, A. G. Species traits in relation to a habitat templet for river systems. Freshw. Biol. 31, 265–275 (1994).Article 

Google Scholar 
Violle, C. et al. Let the concept of trait be functional!. Oikos 116, 882–892 (2007).Article 

Google Scholar 
Schmera, D., Heino, J., Podani, J., Erős, T. & Dolédec, S. Functional diversity: A review of methodology and current knowledge in freshwater macroinvertebrate research. Hydrobiologia 787, 27–44 (2017).Article 

Google Scholar 
Diehl, S. Fish predation and benthic community structure: The role of omnivory and habitat complexity. Ecology 73, 1646–1661 (1992).Article 

Google Scholar 
Brucet, S. et al. Effects of temperature, salinity and fish in structuring the macroinvertebrate community in shallow lakes: Implications for the effects of climate change. PLoS ONE 7, e30877 (2012).ADS 
CAS 
Article 

Google Scholar 
Usseglio-Polatera, P., Bournaud, M., Richoux, P. & Tachet, H. Biological and ecological traits of benthic freshwater macroinvertebrates: Relationships and definition of groups with similar traits. Freshw. Biol. 43, 175–205 (2000).Article 

Google Scholar 
Poff, N. L. et al. Functional trait nichees of North American lotic insects: Trait-based ecological applications in light of phylogenetic relationships. J. North Am. Soc. 25, 730–755 (2006).Article 

Google Scholar 
Bonada, N., Dolédec, S. & Statzner, B. Taxonomic and biological trait differences of stream macroinvertebrate communities between Mediterranean and temperate regions: Implications for future climatic scenarios. Glob. Change Biol. 13, 1658–1671 (2007).ADS 
Article 

Google Scholar 
Stazner, B., Bonada, N. & Dolédec, S. Biological attributes discriminating invasive from native European stream macroinvertebrates. Biol. Invasions 10, 517–530 (2008).Article 

Google Scholar 
Schmidt-Kloiber, A. & Hering, D. www.freshwaterecology.info—An online tool that unifies, standardises and codifies more than 20,000 European freshwater organisms and their ecological preferences. Ecol. Indic. 53, 271–282 (2015).Article 

Google Scholar 
Verberk, W. C. E. P., Siepel, H. & Esselink, H. Life-history strategies in freshwater macroinvertebrates. Freshw. Biol. 53, 1722–1738 (2008).Article 

Google Scholar 
Dolédec, S., Statzner, B. & Frainay, V. Accurate description of functional community structure: Identifying stream invertebrates to species-level?. Bull. North Am. Benthol. Soc. 15, 154–155 (1998).
Google Scholar 
Podani, J., Kalapos, T., Barta, B. & Schmera, D. Principal component analysis of incomplete data—A simple solution to an old problem. Ecol. Inform. 61, 101235 (2021).Article 

Google Scholar 
Podani, J. Introduction to the Exploration of Multivariate Biological Data (Backhuys Publishers, 2000).MATH 

Google Scholar 
Tachet, H., Richoux, P., Bournaud, M. & Usseglio-Polatera, P. Invertébrés d’eau douce—systématique, biologie, écologie 600 (CNRS Editions, 2010).
Google Scholar 
Chevenet, F., Dolédec, S. & Chessel, D. A fuzzy coding approach for the analysis of long-term ecological data. Freshw. Biol. 31, 295–309 (1994).Article 

Google Scholar 
Schmidt-Kloiber A. & Hering, D. www.freshwaterecology.info—The Taxa and Autecology Database for Freshwater Organisms, Version 7.0. (Accessed on 12.09.2019) (2019).Schmera, D., Podani, J., Heino, J., Erős, T. & Poff, N. L. A proposed unified terminology of species traits in stream ecology. Freshw. Sci. 34, 823–830 (2015).Article 

Google Scholar 
Schmera, D., Podani, J., Erős, T. & Heino, J. Combining taxon-by-trait and taxon-by-site matrices for analysing trait patterns of macroinvertebrate communities: A rejoinder to Monaghan & Soares (2014). Freshw. Biol. 59, 1551–1557 (2014).Article 

Google Scholar 
Bonada, N. et al. Do Mediterranean genera not included in Tachet et al. 2002 have Mediterranean characteristics?. Limnetica 30, 129–142 (2011).Article 

Google Scholar 
de Jong, Y. et al. Fauna Europaea—All European animal species on the web. Biodivers. Data J. 2, e4034 (2014).Article 

Google Scholar 
de Bello, F., Botta-Dukát, Z., Leps, J. & Fibich, P. Towards a more balanced combination of multiple traits when computing functional differences between species. Methods Ecol. Evol. 12, 443–448 (2021).Article 

Google Scholar 
Legendre, P. & Legendre, L. Numerical Ecology 2 English. (Elsevier, 1998).MATH 

Google Scholar 
Cornwell, W. K., Schwilk, D. W. & Ackerly, D. A trait-based test for habitat filtering: Convex hull volume. Ecology 87, 1465–1471 (2006).Article 

Google Scholar 
Anderson, M. J. A new method for non-parametric multivariate analysis of variance. Austral. Ecol. 26, 32–46 (2001).
Google Scholar 
Anderson, M. J. Distance-based tests for homogeneity of multivariate dispersions. Biometrics 62, 245–253 (2006).MathSciNet 
Article 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2020).
Google Scholar 
Kindt, R. & Coe, R. Tree Diversity Analysis. A Manual and Software for Common Statistical Methods for Ecological and Biodiversity Studies (World Agroforestry Centre (ICRAF), 2005).
Google Scholar 
Habel, K., Grasman, R., Gramacy, R. G., Mozharovskyi, P. & Sterratt, D. C. Geometry: Mesh Generation and Surface Tessellation. R package version 0.4.5. https://CRAN.R-project.org/package=geometry (2019).Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., Minchin, P. R., O’Hara, R. G., Simpson, G. L., Solymos, P., Stevens, M. H. H., Szoecs, E. & Wagner, H. vegan: Community Ecology Package. R package version 2.5-6. https://CRAN.R-project.org/package=vegan (2019). More

Temporal differentiation of wild populations of T. hispidus and climatic parametersComparison of morphometric features of T. hispidus shells collected in different years in two geographic regions, i.e., Wrocław and Lubawka, showed significant differences depending on the year of collection. The largest number of differences was revealed in shells from Wrocław (Figs. 1 and 2A; Additional file 2: Table S1). Out of 210 comparisons (15 pairs of collection years × 14 features), 84 were statistically significant (Additional file 2: Table S2), e.g., shell diameter (D) was significantly different in 11 cases, shell height (H) and shell width (W) in 10 cases, body whorl height (bwH), the number of whorls (whl), umbilicus major (U) and minor (u) diameters in 9 cases and aperture height/width ratio (h/w) in 7 cases. Nine features obtained more than 10% difference between shells in at least one comparison of mean values, e.g., U 24%, u 19%, H 16% and D 15% (Additional file 2: Table S2). Umbilicus major (U) and minor (u) diameters showed the largest average percentage difference, i.e., 12% and 10%, respectively, in comparisons of all years.Figure 1Shells of Trochulus hispidus collected in different years in Wrocław.Full size imageFigure 2Changes in: mean values of selected morphometric features of shells collected in various years in Wrocław (A) as well as the mean temperature (B) and the relative humidity (C) recorded in four seasons in Wrocław in eight-year period. Abbreviations: D—shell diameter (in mm), H—shell height (in mm), h/w—aperture height/width ratio, whl—number of whorls. The summary statistics for A is included in Table S1 and original data in Table S10 in Additional file 2.Full size imageFor snails from Lubawka, out of 84 comparisons (6 pairs of collection years × 14 features) only 8 were statistically significant (Additional file 2: Table S3). The shells differed significantly in their aperture height (h) and width (w) in 3 comparisons. The h feature showed the percentage difference up to 9% (Additional file 2: Table S3) and the largest average difference was 4.5%.Besides the phenotypic variation, climatic parameters also showed high fluctuations in the studied period (Fig. 2B,C, Additional file 2: Table S4). The maximum difference reported between temperature parameters in some years prior to sample collection in Wrocław was up to 3.7 °C for the maximum winter temperature, while the maximum difference in the relative humidity was up to 11% for autumn. The maximum temperature difference in Jelenia Góra close to Lubawka was up to 3.5 °C for the minimum winter temperature, while the relative humidity differed at most by up to 8% in summer.Differences in shell morphometry under various climatic conditionsThe distinction between shells collected in individual years and changes in climatic parameters along the same period suggest that these differences can be associated with the climate. Therefore, we calculated the average value of a given climatic parameter for each season and studied region and next divided the collected shell data into two groups according to this value. The first group included the shells that developed in conditions above this average and the second below this average (Additional file 2: Table S5). The differences between these groups were statistically significant for 15 out of 16 considered climatic parameters for at least two shell features (Fig. 3). Similarly, each of 14 features significantly separated the groups based on at least two climatic conditions. The results demonstrated that the mean winter temperature substantially influenced nine morphometric shell features, whereas eight characters were changed due to the maximum winter temperature as well as the mean and minimum temperatures in spring, summer and autumn. Umbilicus major (U) and minor (u) diameters as well as umbilicus relative diameter (U/D) were significantly different in 14 pairs of groups characterized by various climatic parameters. In 11 pairs, the height/width ratio (H/W) was significantly different and shell height (H) in 10 pairs.Figure 3Mean percentage differences in morphometric features between shells that were grown in different conditions. The shells were divided into two groups according to the average value of a given climatic parameter for each season and studied region. The first group included the shells that developed in conditions above this average and the second below this average. Positive values indicate that the given feature was greater in the first group, whereas negative values indicate that this feature was greater in the second group. Dendrograms cluster the features and the parameters according to their similarity in the percentage differences. Values marked in bold indicate statistically significant differences between the compared groups of shells. Values at the dendrogram nodes indicate significance assessed according to approximately unbiased test (au) and bootstrap resampling (bp).Full size imageThe umbilicus diameters (u and U) as well as umbilicus relative diameter (U/D) clustered together in the dendrogram based on the mean percentage difference, which indicates that they similarly responded to climatic conditions (Fig. 3). The features u and U revealed the strongest average increase of all features, from 4.1 to 10.5% in shells developed in higher temperatures in all seasons. The largest percentage difference exceeding 10% was recorded for groups separated according to the mean summer and autumn temperatures as well as the maximum summer and minimum autumn temperatures. The U/D ratio was also significantly greater with the mean percentage difference of 2.8–7.6% in shells grown under high temperatures for all seasons and almost all temperature types. On the other hand, the u and U diameters as well as the U/D ratio were on average by 3.7–6.0% significantly smaller in shells developed under higher humidity in summer and winter.The height/width shell ratio (H/W) was grouped with H and bwH features in the dendrogram and was on average by up to 3.6% significantly smaller in shells grown under higher temperatures in all seasons for almost all types of parameters. The maximum winter temperature caused a significant increase, on average by ca. 3%, in shell height (H) and body whorl height (bwH), whereas higher temperatures in other seasons led to their decrease by up to 3.4% (Fig. 3).The shells that were grown in autumn with a relatively high maximum temperature were characterized by ca. 3% significantly smaller aperture height (h) and aperture height/width ratio (h/w), which were clustered together in the dendrogram (Fig. 3).Other four features, shell diameter (D), number of whorls (whl) as well as shell (W) and aperture width (w), formed an additional cluster in the dendrogram (Fig. 3). All of them were on average significantly greater in shells collected one year after winter that was characterized by relatively higher mean and maximum temperatures. The percentage difference was greater, with 3.6–3.9% for W and D.In the dendrogram, the climatic parameters were clustered in several groups indicating their similar influence on the morphometric features of shells (Fig. 3). There are separate clusters for temperature and humidity parameters with the exception of the autumn maximum temperature and autumn humidity, which are grouped together. The other temperature parameters for warmer seasons are separated from those for winter, which indicates that they differently influenced the shell morphometry.Correlations between morphometric shell features and climatic parametersThe influence of climatic conditions on the shells collected in individual years was also assessed using Spearman’s correlation coefficient between the morphometric features and climatic parameters (Fig. 4). Of 224 potential relationships 113 were statistically significant. The spring mean temperature was significantly correlated with 10 morphometric features. Summer humidity and six temperature parameters, i.e., the minimum temperatures as well as the spring and winter maximum temperatures, significantly correlated with eight shell features. Minor umbilicus diameter (u) and umbilicus relative diameter (U/D) were significantly correlated with almost all climatic parameters, i.e., 15, umbilicus major diameter (U) and height/width ratio (H/W) with 13 and the ratio of umbilicus minor to its major diameter (u/U) with 11.Figure 4Spearman’s correlation coefficients between morphometric features of shells with climatic parameters under which the snails were grown. Dendrograms cluster the features and the parameters according to their similarity in the coefficients. Values marked in bold are statistically significant. Values at the dendrogram nodes indicate significance assessed according to approximately unbiased test (au) and bootstrap resampling (bp).Full size imageAs in the case of percentage difference, we can also recognize groups of morphometric features that were similarly correlated with climatic parameters (Fig. 4). Features U/D, U and u were significantly positively correlated with all or almost all temperature parameters for four seasons with the coefficients up to 0.34, 0.30 and 0.36, respectively. On the other hand, the significant correlation coefficients between these features and the humidity in spring, summer and winter were negative and reached − 0.34.Another group of features included shell height/width ratio (H/W), shell height (H) and body whorl height (bwH) (Fig. 4). All of them showed significant negative correlations with all temperature parameters for spring and summer as well as the minimum autumn temperature, and H/W also with the mean and maximum autumn temperatures as well as the mean and minimum winter temperatures. The correlation coefficients reached − 0.28, − 0.27 and − 0.28, respectively. These three features significantly correlated with summer and spring humidity, at up to 0.23.The number of whorls (whl), shell width (W), shell diameter (D), demonstrated a similar correlation with climatic parameters (Fig. 4). They showed the largest and significant correlation coefficients with winter temperatures: up to 0.24, 0.22 and 0.22, respectively. The ratio of umbilicus minor to its major diameter (u/U) showed significant positive correlation up to 0.22 with temperature of warmer seasons.The climatic parameters were grouped into several clusters indicating their similar relationships with morphometric features (Fig. 4). Humidity parameters of warmer seasons formed a separate cluster and temperature parameters were grouped according to seasons. The winter parameters were connected with autumn humidity and separated from temperatures for warmer seasons.Modelling relationships between morphometric shell features and climatic parametersThe joint influence of many climatic parameters on morphometry of shells collected in individual years was studied using a linear mixed-effects (LME) model after exclusion of correlated parameters and a linear ridge regression (LRR) model including all climatic parameters. The latter allows for the inclusion of correlated variables. We separately investigated the seasonal maximum, mean and minimum temperature parameters in combination with seasonal humidity parameters (Additional file 2: Table S6) because they are obviously correlated.Umbilicus minor (u) and major (U) diameters as well as umbilicus relative diameter (U/D) turned out best explained by the climatic parameters (with R2  > 0.15) in two models (Additional file 2: Table S6). Moreover, u, U and U/D were described in LME models by the largest number of significant climatic parameters, i.e., 15. The features u and U had also the largest number of significant parameters in LRR models, i.e., 18 out of 24 possibilities. The largest average values of temperature coefficients for the LRR models were 0.66 for D, 0.58 for W, 0.32 for H, 0.26 for U and 0.22 for u. Thus, all the above-mentioned features were under the strongest influence of the climatic conditions.In the case of LRR models, the coefficients at the winter mean temperature were most often selected as significant, in 12 out of 14 possibilities (Additional file 2: Table S6). The humidity coefficients for autumn were significant in 30 cases of 42 possibilities. The highest average absolute values of coefficients in climatic variables were those for the summer (0.63), spring (0.31) and autumn (0.24) minimum temperatures as well as the summer mean temperature (0.31). Thus, the temperatures of warmer seasons were more important for developing shell morphology. Seasonal humidity coefficients showed similar values compared to each other.Comparison of shell morphometry of T. hispidus and T. sericeus kept under various conditionsIn order to verify the influence of different climatic parameters on Trochulus shell morphometry in selected conditions, we compared shells from three groups of T. hispidus, which represented several subsequent generations: (1) parental snails collected in the wild in Wrocław-Jarnołtów, (2) their offspring bred in the laboratory for two generations and (3) offspring of the second laboratory-bred generation transplanted again into a garden in Wrocław (Fig. 5A–C). The comparison of the group 2 and 1 was to verify if laboratory conditions with controlled temperature and humidity can influence the shell morphometry within only one generation, whereas including the group 3 in the comparison, we wanted to check if snails raised in wild garden conditions can recover the original phenotype. Furthermore, we transplanted into the same garden conditions T. sericeus, which was collected in the wild in Muszkowice (Fig. 5D,E). In this case, we verified if two originally different ecophenotypes T. hispidus and T. sericeus, develop the same shell morphometry under the same conditions.Figure 5Shells of two Trochulus ecophenotypes: parental T. hispidus from wild habitat in Wrocław (A), the first generation of T. hispidus raised in laboratory (B); T. hispidus reared in garden in Wrocław (C); T. sericeus from wild habitat in Muszkowice (D); T. sericeus reared in garden in Wrocław (E).Full size imageConditions in which these snails developed were different. According to WorldClim, the wild environment of T. hispidus in Wrocław was generally warmer than that of T. sericeus in Muszkowice (Additional file 2: Table S7). The largest difference was 1.4 °C for the maximum summer temperature. Relative humidity was lower in Wrocław by up to 2% for warmer seasons but was higher in winter by 1.6%. The difference between the wild and garden localities in Wrocław was much smaller and did not exceed 0.41 °C. The garden conditions were less humid, by up to 2%. However, data from WorldClim are generalized over a longer period and wider regions, so may not well reflect local conditions in the studied places. Actually, the Wrocław site was an open habitat covered with a nettle community like a garden patch, while the Muszkowice site was overgrown by a beech forest, which most likely maintained a higher humidity and a more stable temperature.Laboratory temperatures were substantially different from those in the field, especially for winter (by 18–19.7 °C) as well as for spring and autumn (by 8.2–12 °C). Laboratory humidity was by up to 4.5% lower compared to winter and 5.9–9.9% higher than in spring and summer.A discriminant function analysis (DFA) for the defined groups of snails provided their interesting grouping and separation (Fig. 6). The analysis identified three significant discriminant functions (p  More

Angilletta, M. J., Niewiarowski, P. H. & Navas, C. A. The evolution of thermal physiology in ectotherms. J. Therm. Biol 27, 249–268. https://doi.org/10.1016/S0306-4565(01)00094-8 (2002).Article 

Google Scholar 
Ebersole, J. L., Liss, W. J. & Frissell, C. A. Cold water patches in warm streams: physicochemical characteristics and the influence of shading. JAWRA J. Am. Water Resour. Assoc. 39, 355–368. https://doi.org/10.1111/j.1752-1688.2003.tb04390.x (2003).ADS 
Article 

Google Scholar 
Comte, L. & Grenouillet, G. Do stream fish track climate change? Assessing distribution shifts in recent decades. Ecography 36, 1236–1246. https://doi.org/10.1111/j.1600-0587.2013.00282.x (2013).Article 

Google Scholar 
Kurylyk, B. L., MacQuarrie, K. T. B., Linnansaari, T., Cunjak, R. A. & Curry, R. A. Preserving, augmenting, and creating cold-water thermal refugia in rivers: Concepts derived from research on the Miramichi River, New Brunswick (Canada). Ecohydrology 8, 1095–1108. https://doi.org/10.1002/eco.1566 (2015).Article 

Google Scholar 
Ebersole, J. L., Quiñones, R. M., Clements, S. & Letcher, B. H. Managing climate refugia for freshwater fishes under an expanding human footprint. Front. Ecol. Environ. 18, 271–280. https://doi.org/10.1002/fee.2206 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Caissie, D. The thermal regime of rivers: a review. Freshw. Biol. 51, 1389–1406. https://doi.org/10.1111/j.1365-2427.2006.01597.x (2006).Article 

Google Scholar 
Dick, J. J., Tetzlaff, D. & Soulsby, C. Landscape influence on small-scale water temperature variations in a moorland catchment. Hydrol. Process. 29, 3098–3111. https://doi.org/10.1002/hyp.10423 (2015).ADS 
Article 

Google Scholar 
Fullerton, A. H. et al. Rethinking the longitudinal stream temperature paradigm: region-wide comparison of thermal infrared imagery reveals unexpected complexity of river temperatures. Hydrol. Process. 29, 4719–4737. https://doi.org/10.1002/hyp.10506 (2015).ADS 
Article 

Google Scholar 
Fullerton, A. H. et al. Longitudinal thermal heterogeneity in rivers and refugia for coldwater species: Effects of scale and climate change. Aquatic Sci. 80, 3. https://doi.org/10.1007/s00027-017-0557-9 (2018).Article 

Google Scholar 
Segura, C., Caldwell, P., Sun, G., McNulty, S. & Zhang, Y. A model to predict stream water temperature across the conterminous USA. Hydrol. Process. 29, 2178–2195. https://doi.org/10.1002/hyp.10357 (2015).ADS 
Article 

Google Scholar 
Jonkers, A. R. T. & Sharkey, K. J. The differential warming response of Britain’s rivers (1982–2011). PLOS One 11, e0166247. https://doi.org/10.1371/journal.pone.0166247 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Jackson, F. L., Hannah, D. M., Fryer, R. J., Millar, C. P. & Malcolm, I. A. Development of spatial regression models for predicting summer river temperatures from landscape characteristics: Implications for land and fisheries management. Hydrol. Process. 31, 1225–1238. https://doi.org/10.1002/hyp.11087 (2017).ADS 
Article 

Google Scholar 
Maheu, A., Poff, N. L. & St-Hilaire, A. A classification of stream water temperature regimes in the conterminous USA. River Res. Appl. 32, 896–906. https://doi.org/10.1002/rra.2906 (2016).Article 

Google Scholar 
Steel, E. A., Sowder, C. & Peterson, E. E. Spatial and temporal variation of water temperature regimes on the Snoqualmie River network. J. Am. Water Resour. Assoc. 52, 769–787. https://doi.org/10.1111/1752-1688.12423 (2016).Article 

Google Scholar 
Kearney, M. R., Matzelle, A. & Helmuth, B. Biomechanics meets the ecological niche: The importance of temporal data resolution. J. Exp. Biol. 215, 922–933. https://doi.org/10.1242/jeb.059634 (2012).Article 
PubMed 

Google Scholar 
Burgmer, T., Hillebrand, H. & Pfenninger, M. Effects of climate-driven temperature changes on the diversity of freshwater macroinvertebrates. Oecologia 151, 93–103. https://doi.org/10.1007/s00442-006-0542-9 (2007).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Isaak, D. J., Young, M. K., Nagel, D. E., Horan, D. L. & Groce, M. C. The cold-water climate shield: Delineating refugia for preserving salmonid fishes through the 21st century. Glob. Change Biol. 21, 2540–2553. https://doi.org/10.1111/gcb.12879 (2015).ADS 
Article 

Google Scholar 
Steel, E. A., Beechie, T. J., Torgersen, C. E. & Fullerton, A. H. Envisioning, quantifying, and managing thermal regimes on river networks. Bioscience 67, 506–522. https://doi.org/10.1093/biosci/bix047 (2017).Article 

Google Scholar 
Budescu, D. V. Dominance analysis: A new approach to the problem of relative importance of predictors in multiple regression. Psychol. Bull. 114, 542–551. https://doi.org/10.1037/0033-2909.114.3.542 (1993).Article 

Google Scholar 
Singhal, B. B. S. & Gupta, R. P. Applied Hydrogeology of Fractured Rocks. 2 edn, 408 (Springer, 2010).Shimizu, T. Relation between scanty runoff from mountainous watershed and geology, slope and vegetation (in Japanese with English summary). Bull. Forestry Forest Prod. Res. Inst. 310, 109–128 (1980).
Google Scholar 
Iwasaki, K., Nagasaka, Y. & Nagasaka, A. Geological effects on the scaling relationships of groundwater contributions in Forested Watersheds. Water Resour. Res. 57, e2021WR029641. https://doi.org/10.1029/2021WR029641 (2021).ADS 
Article 

Google Scholar 
Ishiyama, N. et al. The role of geology in creating stream climate-change refugia along climate gradients. bioRxiv, 2022.2005.2002.490355, https://doi.org/10.1101/2022.05.02.490355 (2022).Kanno, Y., Vokoun, J. C. & Letcher, B. H. Paired stream-air temperature measurements reveal fine-scale thermal heterogeneity within headwater brook trout stream networks. River Res. Appl. 30, 745–755. https://doi.org/10.1002/rra.2677 (2014).Article 

Google Scholar 
Snyder, C. D., Hitt, N. P. & Young, J. A. Accounting for groundwater in stream fish thermal habitat responses to climate change. Ecol. Appl. 25, 1397–1419. https://doi.org/10.1890/14-1354.1 (2015).Article 
PubMed 

Google Scholar 
Carslaw, D. C. & Ropkins, K. Openair—an R package for air quality data analysis. Environ. Model. Softw. 27–28, 52–61. https://doi.org/10.1016/j.envsoft.2011.09.008 (2012).Article 

Google Scholar 
Pinheiro, J. C. & Bates, D. M. Mixed-Effects Models in S and S-PLUS. (Springer, 2000).Gelman, A. & Hill, J. Data Analysis Using Regression and Multilevel/Hierarchical Models. (Cambridge University Press, 2006).Harrison, X. A. et al. A brief introduction to mixed effects modelling and multi-model inference in ecology. PeerJ 6, e4794. https://doi.org/10.7717/peerj.4794 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Clarke, P. When can group level clustering be ignored? Multilevel models versus single-level models with sparse data. J. Epidemiol. Commun. Health 62, 752. https://doi.org/10.1136/jech.2007.060798 (2008).CAS 
Article 

Google Scholar 
Theall, K. P. et al. Impact of small group size on neighbourhood influences in multilevel models. J. Epidemiol. Commun. Health 65, 688–695. https://doi.org/10.1136/jech.2009.097956 (2011).Article 

Google Scholar 
Nakagawa, S. & Schielzeth, H. A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods Ecol. Evol. 4, 133–142. https://doi.org/10.1111/j.2041-210x.2012.00261.x (2013).Article 

Google Scholar 
Nakagawa, S., Johnson, P. C. D. & Schielzeth, H. The coefficient of determination R2 and intra-class correlation coefficient from generalized linear mixed-effects models revisited and expanded. J. Royal Soc. Interface 14, 20170213. https://doi.org/10.1098/rsif.2017.0213 (2017).Article 

Google Scholar 
Lüdecke, D., Ben-Shachar, M. S., Patil, I., Waggoner, P. & Makowski, D. performance: An R package for assessment, comparison and testing of statistical models. J. Open Source Softw. 6, 3139. https://doi.org/10.21105/joss.03139 (2021).ADS 
Article 

Google Scholar 
Hair, J. F., Black, W. C., Babin, B. J. & Anderson, R. E. Multivariate Data Analysis: A Global Perspective. 7 edn, (Prentice Hall, 2009).Azen, R. & Budescu, D. V. The dominance analysis approach for comparing predictors in multiple regression. Psychol. Methods 8, 129–148. https://doi.org/10.1037/1082-989x.8.2.129 (2003).Article 
PubMed 

Google Scholar 
Grömping, U. Estimators of relative importance in linear regression based on variance decomposition. Am. Stat. 61, 139–147. https://doi.org/10.1198/000313007X188252 (2007).MathSciNet 
Article 

Google Scholar 
Luo, W. & Azen, R. Determining predictor importance in hierarchical linear models using dominance analysis. J. Educ. Behav. Stat. 38, 3–31. https://doi.org/10.3102/1076998612458319 (2013).Article 

Google Scholar 
R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. (2020).Erickson, T. R. & Stefan, H. G. Linear air/water temperature correlations for streams during open water periods. J. Hydrol. Eng. 5, 317–321. https://doi.org/10.1061/(ASCE)1084-0699(2000)5:3(317) (2000).Article 

Google Scholar 
Webb, B. W., Clack, P. D. & Walling, D. E. Water–air temperature relationships in a Devon river system and the role of flow. Hydrol. Process. 17, 3069–3084. https://doi.org/10.1002/hyp.1280 (2003).ADS 
Article 

Google Scholar 
Gu, Z., Gu, L., Eils, R., Schlesner, M. & Brors, B. Circlize implements and enhances circular visualization in R. Bioinformatics.
30, 2811–2812. https://doi.org/10.1093/bioinformatics/btu393 (2014).Sugimoto, S., Nakamura, F. & Ito, A. Heat budget and statistical analysis of the relationship between stream temperature and riparian forest in the Toikanbetsu River Basin, Northern Japan. J. For. Res. 2, 103–107. https://doi.org/10.1007/BF02348477 (1997).Article 

Google Scholar 
Dugdale, S. J., Malcolm, I. A., Kantola, K. & Hannah, D. M. Stream temperature under contrasting riparian forest cover: Understanding thermal dynamics and heat exchange processes. Sci. Total Environ. 610–611, 1375–1389. https://doi.org/10.1016/j.scitotenv.2017.08.198 (2018).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Timm, A., Ouellet, V. & Daniels, M. Riparian land cover, water temperature variability, and thermal stress for aquatic species in urban streams. Water 13, 2732. https://doi.org/10.3390/w13192732 (2021).Article 

Google Scholar 
Mitchell, S. A simple model for estimating mean monthly stream temperatures after riparian canopy removal. Environ. Manage. 24, 77–83. https://doi.org/10.1007/s002679900216 (1999).CAS 
Article 
PubMed 

Google Scholar 
Horne, J. P. & Hubbart, J. A. A spatially distributed investigation of stream water temperature in a contemporary mixed-land-use watershed. Water 12, 1756. https://doi.org/10.3390/w12061756 (2020).Article 

Google Scholar 
Graham, C. B., Barnard, H. R., Kavanagh, K. L. & McNamara, J. P. Catchment scale controls the temporal connection of transpiration and diel fluctuations in streamflow. Hydrol. Process. 27, 2541–2556. https://doi.org/10.1002/hyp.9334 (2013).ADS 
Article 

Google Scholar 
Sun, H., Kasahara, T., Otsuki, K., Saito, T. & Onda, Y. Spatio-temporal streamflow generation in a small, steep headwater catchment in Western Japan. Hydrol. Sci. J. 62, 818–829. https://doi.org/10.1080/02626667.2016.1266635 (2017).Article 

Google Scholar 
Sophocleous, M. Interactions between groundwater and surface water: The state of the science. Hydrogeol. J. 10, 52–67. https://doi.org/10.1007/s10040-001-0170-8 (2002).ADS 
CAS 
Article 

Google Scholar 
Arnott, S., Hilton, J. & Webb, B. W. The impact of geological control on flow accretion in lowland permeable catchments. Hydrol. Res. 40, 533–543. https://doi.org/10.2166/nh.2009.017 (2009).Article 

Google Scholar 
Calvache, M. L., Duque, C., Fontalva, J. M. G. & Crespo, F. Processes affecting groundwater temperature patterns in a coastal aquifer. Int. J. Environ. Sci. Technol. 8, 223–236. https://doi.org/10.1007/BF03326211 (2011).Article 

Google Scholar 
Nejadhashemi, A. P., Wardynski, B. J. & Munoz, J. D. Evaluating the impacts of land use changes on hydrologic responses in the agricultural regions of Michigan and Wisconsin. Hydrol. Earth Syst. Sci. 2011, 3421–3468, https://doi.org/10.5194/hessd-8-3421-2011 (2011).Macedo, M. N. et al. Land-use-driven stream warming in southeastern Amazonia. Philos. Trans. R Soc. Lond. B Biol. Sci. 368, 20120153–20120153. https://doi.org/10.1098/rstb.2012.0153 (2013).Article 
PubMed 
PubMed Central 

Google Scholar 
Carlson, K. M. et al. Influence of watershed-climate interactions on stream temperature, sediment yield, and metabolism along a land use intensity gradient in Indonesian Borneo. J. Geophys. Res. Biogeosci. 119, 1110–1128. https://doi.org/10.1002/2013JG002516 (2014).Article 

Google Scholar 
Wang, Y. I., He, B. I. N. & Takase, K. Effects of temporal resolution on hydrological model parameters and its impact on prediction of river discharge. Hydrol. Sci. J. 54, 886–898. https://doi.org/10.1623/hysj.54.5.886 (2009).Article 

Google Scholar 
Levin, S. A. The problem of pattern and scale in ecology: The Robert H MacArthur award lecture. Ecology 73, 1943–1967. https://doi.org/10.2307/1941447 (1992).Article 

Google Scholar 
García Molinos, J. & Donohue, I. Downscaling the non-stationary effect of climate forcing on local-scale dynamics: The importance of environmental filters. Clim. Change 124, 333–346. https://doi.org/10.1007/s10584-014-1077-4 (2014).ADS 
Article 

Google Scholar 
Newman, E. A., Kennedy, M. C., Falk, D. A. & McKenzie, D. Scaling and complexity in landscape ecology. Front. Ecol. Evolution https://doi.org/10.3389/fevo.2019.00293 (2019).Article 

Google Scholar 
Atkinson, S. E., Woods, R. A. & Sivapalan, M. Climate and landscape controls on water balance model complexity over changing timescales. Water Resour. Res. 38, 50-51–50-17, https://doi.org/10.1029/2002WR001487 (2002).Engel, M. et al. Controls on spatial and temporal variability in streamflow and hydrochemistry in a glacierized catchment. Hydrol. Earth Syst. Sci. 23, 2041–2063. https://doi.org/10.5194/hess-23-2041-2019 (2019).ADS 
Article 

Google Scholar 
Karlsen, R. H. et al. Landscape controls on spatiotemporal discharge variability in a boreal catchment. Water Resour. Res. 52, 6541–6556. https://doi.org/10.1002/2016WR019186 (2016).ADS 
Article 

Google Scholar 
Parmesan, C. & Yohe, G. A globally coherent fingerprint of climate change impacts across natural systems. Nature 421, 37–42. https://doi.org/10.1038/nature01286 (2003).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Weiskopf, S. R. et al. Climate change effects on biodiversity, ecosystems, ecosystem services, and natural resource management in the United States. Sci. Total Environ. 733, 137782. https://doi.org/10.1016/j.scitotenv.2020.137782 (2020).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Radchuk, V. et al. Adaptive responses of animals to climate change are most likely insufficient. Nat. Commun. 10, 3109. https://doi.org/10.1038/s41467-019-10924-4 (2019).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Kingsford, R. T. Conservation management of rivers and wetlands under climate change—a synthesis. Mar. Freshw. Res. 62, 217–222. https://doi.org/10.1071/MF11029 (2011).CAS 
Article 

Google Scholar  More

Hoegh-Guldberg, O. et al. Coral reefs under rapid climate change and ocean acidification. Science 318, 1737–1742. https://doi.org/10.1126/science.1152509 (2007).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Hughes, T. P. et al. Climate change, human impacts, and the resilience of coral reefs. Science 301, 929–933. https://doi.org/10.1126/science.1085046 (2003).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Hughes, T. P. et al. Global warming and recurrent mass bleaching of corals. Nature 543, 373–377. https://doi.org/10.1038/nature21707 (2017).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Bahr, K. D., Jokiel, P. L. & Toonen, R. J. The unnatural history of Kāne‘ohe Bay: coral reef resilience in the face of centuries of anthropogenic impacts. PeerJ 3, e950. https://doi.org/10.7717/peerj.950 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Montilla, L. M., Ascanio, A., Verde, A. & Croquer, A. Systematic review and meta-analysis of 50 years of coral disease research visualized through the scope of network theory. PeerJ 7, e7041. https://doi.org/10.7717/peerj.7041 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
Williams, S. D., Walter, C. S. & Muller, E. M. Fine Scale temporal and spatial dynamics of the stony coral tissue loss disease outbreak within the lower Florida keys. Front. Mar. Sci. https://doi.org/10.3389/fmars.2021.631776 (2021).Article 

Google Scholar 
Harrison, P. L. in Coral Reefs: An Ecosystem in Transition (eds Zvy Dubinsky & Noga Stambler) 59–85 (Springer Netherlands, 2011).Richmond, R. H. & Hunter, C. L. Reproduction and recruitment of corals: comparisons among the Caribbean, the Tropical Pacific, and the Red Sea. Mar. Ecol. Prog. Ser. 60, 185–203 (1990).ADS 
Article 

Google Scholar 
Humphrey, C., Weber, M., Lott, C., Cooper, T. & Fabricius, K. Effects of suspended sediments, dissolved inorganic nutrients and salinity on fertilisation and embryo development in the coral Acropora millepora (Ehrenberg, 1834). Coral Reefs 27, 837–850. https://doi.org/10.1007/s00338-008-0408-1 (2008).ADS 
Article 

Google Scholar 
Fabricius, K. E. Effects of terrestrial runoff on the ecology of corals and coral reefs: review and synthesis. Mar. Pollut. Bull. 50, 125–146. https://doi.org/10.1016/j.marpolbul.2004.11.028 (2005).CAS 
Article 
PubMed 

Google Scholar 
Ayalon, I. et al. Coral gametogenesis collapse under artificial light pollution. Curr. Biol. 31, 413-419.e413. https://doi.org/10.1016/j.cub.2020.10.039 (2021).CAS 
Article 
PubMed 

Google Scholar 
Szmant, A. M. & Gassman, N. J. The effects of prolonged “bleaching” on the tissue biomass and reproduction of the reef coral Montastrea annularis. Coral Reefs 8, 217–224. https://doi.org/10.1007/BF00265014 (1990).ADS 
Article 

Google Scholar 
Fisch, J., Drury, C., Towle, E. K., Winter, R. N. & Miller, M. W. Physiological and reproductive repercussions of consecutive summer bleaching events of the threatened Caribbean coral Orbicella faveolata. Coral Reefs 38, 863–876. https://doi.org/10.1007/s00338-019-01817-5 (2019).ADS 
Article 

Google Scholar 
Hughes, T. P. et al. Global warming impairs stock–recruitment dynamics of corals. Nature 568, 387–390. https://doi.org/10.1038/s41586-019-1081-y (2019).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Barfield, S., Aglyamova, G. V. & Matz, M. V. Evolutionary origins of germline segregation in Metazoa: evidence for a germ stem cell lineage in the coral Orbicella faveolata (Cnidaria, Anthozoa). Proc. Biol. Sci. https://doi.org/10.1098/rspb.2015.2128 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Highsmith, R. C. Reproduction by fragmentation in corals. Mar. Ecol. Prog. Ser. 7, 207–226 (1982).ADS 
Article 

Google Scholar 
Baums, I. B. A restoration genetics guide for coral reef conservation. Mol. Ecol. 17, 2796–2811. https://doi.org/10.1111/j.1365-294X.2008.03787.x (2008).Article 
PubMed 

Google Scholar 
Babcock, R. C. et al. Synchronous spawnings of 105 scleractinian coral species on the Great Barrier Reef. Mar. Biol. 90, 379–394. https://doi.org/10.1007/BF00428562 (1986).Article 

Google Scholar 
Baird, A. H., Guest, J. R. & Willis, B. L. Systematic and biogeographical patterns in the reproductive biology of scleractinian corals. Annu. Rev. Ecol. Evol. Syst. 40, 551–571. https://doi.org/10.1146/annurev.ecolsys.110308.120220 (2009).Article 

Google Scholar 
Levitan, D. R. & Petersen, C. Sperm limitation in the sea. Trends Ecol. Evol. 10, 228–231. https://doi.org/10.1016/S0169-5347(00)89071-0 (1995).CAS 
Article 
PubMed 

Google Scholar 
Yund, P. O. How severe is sperm limitation in natural populations of marine free-spawners?. Trends Ecol. Evol. 15, 10–13. https://doi.org/10.1016/S0169-5347(99)01744-9 (2000).CAS 
Article 
PubMed 

Google Scholar 
Bouwmeester, J. et al. Multi-species spawning synchrony within scleractinian coral assemblages in the Red Sea. Coral Reefs 34, 65–77. https://doi.org/10.1007/s00338-014-1214-6 (2015).ADS 
Article 

Google Scholar 
Guest, J. R., Baird, A. H., Goh, B. P. L. & Chou, L. M. Reproductive seasonality in an equatorial assemblage of scleractinian corals. Coral Reefs 24, 112–116. https://doi.org/10.1007/s00338-004-0433-7 (2005).Article 

Google Scholar 
Chelliah, A. et al. First record of multi-species synchronous coral spawning from Malaysia. PeerJ 3, e777. https://doi.org/10.7717/peerj.777 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Hunter, C. L. in Proceedings of the 6th International Coral Reef Symposium Vol. 2, 727–732 (1988).Jokiel, P. L., Ito, R. Y. & Liu, P. M. Night irradiance and synchronization of lunar release of planula larvae in the reef coral Pocillopora damicornis. Mar. Biol. 88, 167–174. https://doi.org/10.1007/BF00397164 (1985).Article 

Google Scholar 
Willis, B. L., Babcock, R. C., Harrison, P. L., Oliver, J. K. & Wallace, C. C. in Proceedings of the 5th International Coral Reef Congress Vol. 4, 343–348 (1985).Brady, A. K., Hilton, J. D. & Vize, P. D. Coral spawn timing is a direct response to solar light cycles and is not an entrained circadian response. Coral Reefs 28, 677–680. https://doi.org/10.1007/s00338-009-0498-4 (2009).ADS 
Article 

Google Scholar 
Mendes, J. M. & Woodley, J. D. Timing of reproduction in Montastraea annularis: relationship to environmental variables. Mar. Ecol. Prog. Ser. 227, 241–251. https://doi.org/10.3354/meps227241 (2002).ADS 
Article 

Google Scholar 
van Woesik, R. Calm before the spawn: global coral spawning patterns are explained by regional wind fields. Proc. Biol. Sci. 277, 715–722. https://doi.org/10.1098/rspb.2009.1524 (2010).Article 
PubMed 

Google Scholar 
Twan, W.-H. et al. Hormones and reproduction in scleractinian corals. Comp. Biochem. Physiol. A Mol. Integr. Physiol. 144, 247–253. https://doi.org/10.1016/j.cbpa.2006.01.011 (2006).CAS 
Article 
PubMed 

Google Scholar 
Tan, E. S., Izumi, R., Takeuchi, Y., Isomura, N. & Takemura, A. Molecular approaches underlying the oogenic cycle of the scleractinian coral, Acropora tenuis. Sci. Rep. 10, 9914. https://doi.org/10.1038/s41598-020-66020-x (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Baird, A. H. & Marshall, P. A. Mortality, growth and reproduction in scleractinian corals following bleaching on the Great Barrier Reef. Mar. Ecol. Prog. Ser. 237, 133–141. https://doi.org/10.3354/meps237133 (2002).ADS 
Article 

Google Scholar 
Levitan, D. R., Boudreau, W., Jara, J. & Knowlton, N. Long-term reduced spawning in Orbicella coral species due to temperature stress. Mar. Ecol. Prog. Ser. 515, 1–10. https://doi.org/10.2307/24894795 (2014).ADS 
Article 

Google Scholar 
Hagedorn, M. et al. Potential bleaching effects on coral reproduction. Reprod. Fertil. Dev. https://doi.org/10.1071/rd15526 (2016).Article 

Google Scholar 
Ward, S., Harrison, P. L. & Hoegh-Guldberg, O. in Proceedings of the Ninth International Coral Reef Symposium Vol. 2, 1123–1128 (2002).Lager, C. V. A., Hagedorn, M., Rodgers, K. S. & Jokiel, P. L. The impact of short-term exposure to near shore stressors on the early life stages of the reef building coral Montipora capitata. PeerJ 8, e9415. https://doi.org/10.7717/peerj.9415 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Vermeij, M. J. A., Fogarty, N. D. & Miller, M. W. Pelagic conditions affect larval behavior, survival, and settlement patterns in the Caribbean coral Montastraea faveolata. Mar. Ecol. Prog. Ser. 310, 119–128. https://doi.org/10.3354/meps310119 (2006).ADS 
Article 

Google Scholar 
Torres, J. L., Armstrong, R. A. & Weil, E. Enhanced ultraviolet radiation can terminate sexual reproduction in the broadcasting coral species Acropora cervicornis (Lamarck). J. Exp. Mar. Biol. Ecol. 358, 39–45. https://doi.org/10.1016/j.jembe.2008.01.022 (2008).Article 

Google Scholar 
Wellington, G. M. & Fitt, W. K. Influence of UV radiation on the survival of larvae from broadcast-spawning reef corals. Mar. Biol. 143, 1185–1192. https://doi.org/10.1007/s00227-003-1150-4 (2003).CAS 
Article 

Google Scholar 
Paxton, C. W., Baria, M. V. B., Weis, V. M. & Harii, S. Effect of elevated temperature on fecundity and reproductive timing in the coral Acropora digitifera. Zygote 24, 511–516. https://doi.org/10.1017/S0967199415000477 (2016).Article 
PubMed 

Google Scholar 
Nozawa, Y. Annual variation in the timing of coral spawning in a high-latitude environment: influence of temperature. Biol. Bull. 222, 192–202. https://doi.org/10.1086/BBLv222n3p192 (2012).Article 
PubMed 

Google Scholar 
Bouwmeester, J. et al. Latitudinal variation in monthly-scale reproductive synchrony among Acropora coral assemblages in the Indo-Pacific. Coral Reefs https://doi.org/10.1007/s00338-021-02129-3 (2021).Article 

Google Scholar 
Mendes, J. M. & Woodley, J. D. Effect of the 1995–1996 bleaching event on polyp tissue depth, growth, reproduction and skeletal band formation in Montastraea annularis. Mar. Ecol. Prog. Ser. 235, 93–102. https://doi.org/10.3354/meps235093 (2002).ADS 
Article 

Google Scholar 
Johnston, E. C., Counsell, C. W. W., Sale, T. L., Burgess, S. C. & Toonen, R. J. The legacy of stress: Coral bleaching impacts reproduction years later. Funct. Ecol. 00, 1–11. https://doi.org/10.1111/1365-2435.13653 (2020).Article 

Google Scholar 
Holbrook, S. J. et al. Recruitment drives spatial variation in recovery rates of resilient coral reefs. Sci. Rep. 8, 7338. https://doi.org/10.1038/s41598-018-25414-8 (2018).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
de la Cruz, D. W. & Harrison, P. L. Enhanced larval supply and recruitment can replenish reef corals on degraded reefs. Sci. Rep. 7, 13985. https://doi.org/10.1038/s41598-017-14546-y (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Villanueva, R. D., Baria, M. V. B. & de la Cruz, D. W. Growth and survivorship of juvenile corals outplanted to degraded reef areas in Bolinao-Anda Reef Complex, Philippines. Mar. Biol. Res. 8, 877–884. https://doi.org/10.1080/17451000.2012.682582 (2012).Article 

Google Scholar 
Chamberland, V. F. et al. Restoration of critically endangered elkhorn coral (Acropora palmata) populations using larvae reared from wild-caught gametes. Global Ecol. Conserv. 4, 526–537. https://doi.org/10.1016/j.gecco.2015.10.005 (2015).Article 

Google Scholar 
Hunter, C. L. & Evans, C. W. Coral reefs in Kaneohe Bay, Hawaii: two centuries of western influence and two decades of data. Bull. Mar. Sci. 57, 501–515 (1995).
Google Scholar 
Rodgers, K. S., Jokiel, P. L., Brown, E. K., Hau, S. & Sparks, R. Over a decade of change in spatial and temporal dynamics of Hawaiian coral reef communities. Pac. Sci. 69, 1–13. https://doi.org/10.2984/69.1.1 (2015).Article 

Google Scholar 
Maragos, J. E. et al. 2000–2002 rapid ecological assessment of corals (Anthozoa) on shallow reefs of the Northwestern Hawaiian Islands. Part 1: species and distribution. Pac. Sci. 58, 211–230. https://doi.org/10.1353/psc.2004.0020 (2004).Article 

Google Scholar 
Richards Donà, A. Investigation into the functional role of chromoproteins in the physiology and ecology of the Hawaiian stony coral Montipora flabellata in Kāne‘ohe Bay, O‘ahu, University of Hawaiʻi at Mānoa, (2019).Padilla-Gamiño, J. L. & Gates, R. D. Spawning dynamics in the Hawaiian reef-building coral Montipora capitata. Mar. Ecol. Prog. Ser. 449, 145–160. https://doi.org/10.3354/meps09530 (2012).ADS 
Article 

Google Scholar 
Krupp, D. A. Sexual reproduction and early development of the solitary coral Fungia scutaria (Anthozoa: Scleractinia). Coral Reefs 2, 159–164. https://doi.org/10.1007/BF00336722 (1983).ADS 
Article 

Google Scholar 
Kolinski, S. P. & Cox, E. F. An update on modes and timing of gamete and planula release in Hawaiian scleractinian corals with implications for conservation and management. Pac. Sci. 57, 17–27. https://doi.org/10.1353/psc.2003.0005 (2003).Article 

Google Scholar 
Heyward, A. J. Sexual reproduction in five species of the coral Montipora. Coral Reef Popul. Biol. Hawaii Inst. Mar. Biol. Tech. Rep. 37, 170–178 (1985).
Google Scholar 
Marquis, R. J. Phenological variation in the neotropical understory shrub Piper arielanum: causes and consequences. Ecology 69, 1552–1565. https://doi.org/10.2307/1941653 (1988).Article 

Google Scholar 
Padilla-Gamiño, J. L. et al. Sedimentation and the reproductive biology of the Hawaiian reef-building coral Montipora capitata. Biol. Bull. 226, 8–18. https://doi.org/10.1086/BBLv226n1p8 (2014).Article 
PubMed 

Google Scholar 
Humason, G. L. Animal Tissue Techniques. 661 (W. H. Freeman & Co, 1979).Abramoff, M. D., Magalhães, P. J. & Ram, S. J. Image processing with ImageJ. Biophotonics Int. 11, 36–42 (2004).
Google Scholar 
Szmant-Froelich, A., Reutter, M. & Riggs, L. Sexual reproduction of Favia fragum (Esper): lunar patterns of gametogenesis, embryogenesis and planulation in Puerto Rico. Bull. Mar. Sci. 37, 880–892 (1985).
Google Scholar 
Vargas-Ángel, B., Colley, S. B., Hoke, S. M. & Thomas, J. D. The reproductive seasonality and gametogenic cycle of Acropora cervicornis off Broward County, Florida, USA. Coral Reefs 25, 110–122. https://doi.org/10.1007/s00338-005-0070-9 (2006).ADS 
Article 

Google Scholar 
Baird, A. H., Blakeway, D. R., Hurley, T. J. & Stoddart, J. A. Seasonality of coral reproduction in the Dampier Archipelago, northern Western Australia. Mar. Biol. 158, 275–285. https://doi.org/10.1007/s00227-010-1557-7 (2011).Article 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, Vienna, Austria, https://www.R-project.org, 2019).An {R} Companion to Applied Regression. Third Edition (Sage, Thousand Oaks (CA), 2019).Mangiafico, S. rcompanion: Functions to Support Extension Education Program Evaluation. R package version 2.3.7. https://CRAN.R-project.org/package=rcompanion. (2019).Bahr, K. D., Jokiel, P. L. & Rodgers, K. S. The 2014 coral bleaching and freshwater flood events in Kāneʻohe Bay, Hawaiʻi. PeerJ 3, e1136. https://doi.org/10.7717/peerj.1136 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Bahr, K. D., Rodgers, K. S. & Jokiel, P. L. Impact of three bleaching events on the reef resiliency of Kāne‘ohe Bay, Hawai‘i. Front. Mar. Sci. 4, 398. https://doi.org/10.3389/fmars.2017.00398 (2017).Article 

Google Scholar 
Bachtiar, I. Reproduction of three scleractinian corals (Acropora cytherea, A. nobilis, Hydnophora rigida) in easter Lombok Strait, Indonesia. Indones. J. Mar. Sci. 6, 18–27 (2001).
Google Scholar 
Baird, A. H., Marshall, P. A. & Wolstenholme, J. in Proceedings of the 9th International Coral Reef Symposium Vol. 1, 385–389 (2002).Mangubhai, S. & Harrison, P. L. Asynchronous coral spawning patterns on equatorial reefs in Kenya. Mar. Ecol. Prog. Ser. 360, 85–96. https://doi.org/10.3354/meps07385 (2008).ADS 
Article 

Google Scholar 
Prasetia, R., Sinniger, F. & Harii, S. Gametogenesis and fecundity of Acropora tenella (Brook 1892) in a mesophotic coral ecosystem in Okinawa, Japan. Coral Reefs 35, 53–62. https://doi.org/10.1007/s00338-015-1348-1 (2016).ADS 
Article 

Google Scholar 
Parker, G. A. Why are there so many tiny sperm? Sperm competition and the maintenance of two sexes. J. Theor. Biol. 96, 281–294. https://doi.org/10.1016/0022-5193(82)90225-9 (1982).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Hayward, A. & Gillooly, J. F. The cost of sex: quantifying energetic investment in gamete production by males and females. PLoS ONE 6, e16557. https://doi.org/10.1371/journal.pone.0016557 (2011).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Fenner, D. P. Corals of Hawai’i. A field guide to the hard, black, and soft corals of Hawai’i and the northwest Hawaiian Islands, including Midway (Mutual Publishing Company, 2005).
Google Scholar 
Darling, E. S., Alvarez-Filip, L., Oliver, T. A., McClanahan, T. R. & Côté, I. M. Evaluating life-history strategies of reef corals from species traits. Ecol. Lett. 15, 1378–1386. https://doi.org/10.1111/j.1461-0248.2012.01861.x (2012).Article 
PubMed 

Google Scholar 
Okubo, N., Motokawa, T. & Omori, M. When fragmented coral spawn? effect of size and timing on survivorship and fecundity of fragmentation in Acropora formosa. Mar. Biol. 151, 353–363. https://doi.org/10.1007/s00227-006-0490-2 (2006).Article 

Google Scholar 
Szmant-Froelich, A., Yevich, P. & Pilson, M. E. Q. Gametogenesis and early development of the temperate coral Astrangia danae (Anthozoa: Scleractinia). Biol. Bull. 158, 257–269. https://doi.org/10.2307/1540935 (1980).Article 

Google Scholar 
Kojis, B. L. Sexual reproduction in Acropora (Isopora) (Coelenterata: Scleractinia). Mar. Biol. 91, 311–318. https://doi.org/10.1007/BF00428624 (1986).Article 

Google Scholar 
Neves, E. & Pires, D. Sexual reproduction of Brazilian coral Mussismilia hispida (Verrill, 1902). Coral Reefs 21, 161–168. https://doi.org/10.1007/s00338-002-0217-x (2002).Article 

Google Scholar 
Pennington, J. T. The ecology of fertilization of Echinoid eggs: the consequences of sperm dilution, adult aggregation, and synchronous spawning. Biol. Bull. 169, 417–430. https://doi.org/10.2307/1541492 (1985).Article 
PubMed 

Google Scholar 
Oliver, J. & Babcock, R. C. Aspects of the fertilization ecology of broadcast spawning corals: sperm dilution effects and in situ measurements of fertilization. Biol. Bull. 183, 409–417. https://doi.org/10.2307/1542017 (1992).CAS 
Article 
PubMed 

Google Scholar 
Lasker, H. R. et al. In situ rates of fertilization among broadcast spawning Gorgonian corals. Biol. Bull. 190, 45–55. https://doi.org/10.2307/1542674 (1996).CAS 
Article 
PubMed 

Google Scholar 
Gorbunov, M. Y. & Falkowski, P. G. Photoreceptors in the cnidarian hosts allow symbiotic corals to sense blue moonlight. Limnol. Oceanogr. 47, 309–315. https://doi.org/10.4319/lo.2002.47.1.0309 (2002).ADS 
Article 

Google Scholar 
van Woesik, R., Lacharmoise, F. & Köksal, S. Annual cycles of solar insolation predict spawning times of Caribbean corals. Ecol. Lett. 9, 390–398. https://doi.org/10.1111/j.1461-0248.2006.00886.x (2006).Article 
PubMed 

Google Scholar 
Wolstenholme, J. K. Temporal reproductive isolation and gametic compatibility are evolutionary mechanisms in the Acropora humilis species group (Cnidaria; Scleractinia). Mar. Biol. 144, 567–582. https://doi.org/10.1007/s00227-003-1209-2 (2004).Article 

Google Scholar 
Colley, S. B., Feingold, J. S., Peña, J. & Glynn, P. W. in Proceedings of the 9th International Coral Reef Symposium Vol. 1, 23–27 (2000).Chaves-Fonnegra, A., Maldonado, M., Blackwelder, P. & Lopez, J. V. Asynchronous reproduction and multi-spawning in the coral-excavating sponge Cliona delitrix. J. Mar. Biol. Assoc. U.K. 96, 515–528. https://doi.org/10.1017/S0025315415000636 (2016).Article 

Google Scholar 
Chamberland, V. F., Snowden, S., Marhaver, K. L., Petersen, D. & Vermeij, M. J. A. The reproductive biology and early life ecology of a common Caribbean brain coral, Diploria labyrinthiformis (Scleractinia: Faviinae). Coral Reefs 36, 83–94. https://doi.org/10.1007/s00338-016-1504-2 (2017).ADS 
Article 

Google Scholar 
Sherman, C. D. H. Mating system variation in the hermaphroditic brooding coral, Seriatopora hystrix. Heredity 100, 296–303. https://doi.org/10.1038/sj.hdy.6801076 (2008).CAS 
Article 
PubMed 

Google Scholar 
Yeoh, S.-R. & Dai, C.-F. The production of sexual and asexual larvae within single broods of the scleractinian coral, Pocillopora damicornis. Mar. Biol. 157, 351–359. https://doi.org/10.1007/s00227-009-1322-y (2010).Article 

Google Scholar 
Henley, E. M. et al. Reproductive plasticity of Hawaiian Montipora corals following thermal stress. Sci. Rep. 11, 12525. https://doi.org/10.1038/s41598-021-91030-8 (2021).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Coma, R. & Lasker, H. R. Effects of spatial distribution and reproductive biology on in situ fertilization rates of a broadcast-spawning invertebrate. Biol. Bull. 193, 20–29. https://doi.org/10.2307/1542733 (1997).CAS 
Article 
PubMed 

Google Scholar 
Westneat, M. W. & Resing, J. M. Predation on coral spawn by planktivorous fish. Coral Reefs 7, 89–92. https://doi.org/10.1007/BF00301646 (1988).ADS 
Article 

Google Scholar 
Fitzhugh, G. R., Shertzer, K. W., Kellison, G. T. & Wyanski, D. M. Review of size- and age-dependence in batch spawning: implications for stock assessment of fish species exhibiting indeterminate fecundity. Fish. Bull. 110, 413–425 (2012).
Google Scholar 
Alvarado, E. M., García, R. & Acosta, A. Sexual reproduction of the reef-building coral Diploria labyrinthiformis (Scleractinia:Faviidae), in the Colombian Caribbean. Rev. Biol. Trop. 52, 859–868 (2004).PubMed 

Google Scholar 
Maragos, J. E. A Study of the Ecology of Hawaiian Reef Corals, University of Hawaiʻi at Mānoa, (1972).Jokiel, P. L. & Brown, E. Global warming, regional trends and inshore environmental conditions influence coral bleaching in Hawaii. Glob. Chang. Biol. 10, 1627–1641. https://doi.org/10.1111/j.1365-2486.2004.00836.x (2004).ADS 
Article 

Google Scholar 
Courchamp, F., Clutton-Brock, T. & Grenfell, B. Inverse density dependence and the Allee effect. Trends Ecol. Evol. 14, 405–410. https://doi.org/10.1016/S0169-5347(99)01683-3 (1999).CAS 
Article 
PubMed 

Google Scholar 
Berec, L., Angulo, E. & Courchamp, F. Multiple Allee effects and population management. Trends Ecol. Evol. 22, 185–191. https://doi.org/10.1016/j.tree.2006.12.002 (2007).Article 
PubMed 

Google Scholar 
Hughes, T. P. et al. Spatial and temporal patterns of mass bleaching of corals in the Anthropocene. Science 359, 80. https://doi.org/10.1126/science.aan8048 (2018).ADS 
CAS 
Article 
PubMed 

Google Scholar  More

For a schematic overview of the experimental design, see Fig. 2.Experimental evolutionProtocolThe stock population (Stock population below) was allowed to expand for one generation and from this we established eight replicate experimental evolution populations, four selected for fighter morphs (F-lines) and four selected for scrambler morphs (S-lines). Each population was founded by 1,000 recently eclosed adults (500 random females and 500 random males of the desired morph). The classification of the morphs was based on visual inspection using a stereoscopic microscope and was unambiguous due to the discontinuous distribution of the phenotypes (Classifying male morphs below). Adults were allowed to interact freely for 6 days, all surviving adults (with previously laid eggs discarded) were transferred to a new container for 24 h of egg laying, after which adults were removed. The resulting offspring were allowed to mature over 13 days and 1,000 individuals from the newly eclosed adults selected for founding the following generation, again 500 random females and 500 random males of the desired morph, with this protocol repeated every generation (Extended Data Fig. 1). The isolation of nymphs to use virgins was unfeasible with our experimental design and population sizes. However, the period of 6 days after selecting the founders of the next generation and collecting eggs for the next generation was probably enough to displace most sperm stored by females mated with any unselected males due to the high number of remating that will be occurring over this duration (females on average remate after 80 min, ref. 88) and last male sperm precedence89. The timing of generation was chosen to reflect maturation rates from our stock population to avoid indirect selection on this trait. Moreover, a previous study90 showed there was no difference between male morphs in maturation rates and that over similar lengths of time to the protocol here the fertility of both morphs remains similar. Therefore, our protocol was not likely to impose strong differential selection on morph life histories.Tracking morph proportionWe assayed the proportion of male morph in each population every 6–7 generations, by isolating 200 larvae (ten per vial) from the container, allowing maturation within vials and recording the morph of all males that eclosed (mean n = 86 per population, per generation, range 71–109). Our selection protocol was highly effective in driving an increase in the frequency of the desired male morph to >90% after 20 generations in both treatments, with this effect considerably faster within F-lines indicted by a significant two-way interaction between proportion of the desired morph and generation (χ2 = 39.9, d.f. = 6, P 90% is probably a consequence of a longer interaction period (3 versus 6 days) in which the stored sperm of males before selection was able to be displaced and/or because selection was acting more efficiently in our larger populations. The difference between rate of changes in morph proportion between F- and S-lines in the current study, and also found by Plesnar-Bielak et al.39, may be associated with the genetic architecture of morph expression. Alternatively, selection could be less effective in scrambler lines if they are less efficient than fighters in displacing sperm of previous females’ partners, but this is unlikely as R. robini male morphs have previously been demonstrated to not differ in their sperm competiveness89.Stock populationWe established a stock population by mixing three laboratory populations that were collected from three sites in Poland (Krakόw, collected in 1998 and 2008, Kwiejce, collected in 2017 and Mosina, collected in 2017; Extended Data Fig. 1), where the line derived for material used in creating the reference genome (below) was also established from the same collections at Mosina in 2017. All populations were maintained in cultures with several hundred individuals per generation before mixing and establishment of the stock population. The mixing of distinct populations increased the genetic variance in the stock population, which otherwise would probably have been limited due to founder events and the limited population size of each of the contributing populations73, thus decreasing our power to detect the effects of SSTs on genetic variation. The newly mixed stock population was maintained with several hundred individuals per generation for roughly 12 generations before the onset of this experiment. This time period is probably enough to break linkage disequilibria that could have arisen due to mixing (for unlinked loci, linkage disequilibrium should decay by half each generation91).One generation before establishing experimental evolution populations the proportion of male morphs was determined from 176 random males, indicating a roughly equal morph ratio (95 fighters, 81 scramblers) of the stock population (Extended Data Fig. 3).General housing and husbandryThe stock population and experimental evolution populations were maintained in plastic containers (approximate length, 9.5 cm; width, 7 cm; height, 4.5 cm), filled with roughly 1 cm of plaster-of-Paris. The same containers were used when sampling mites for sequencing for the reference genome or resequencing from experimental evolution populations, but either replaced the plaster-of-Paris with 5% agarose gel or added a thin layer of 5% agarose gel above the plaster-of-Paris, respectively. The agarose gel was used to reduce the number of contaminates within our samples and on the basis of preliminary extractions that indicated that small pieces of plaster-of-Paris may reduce the quality of DNA during extractions. Individuals, pairs and small groups of ten mites were housed in glass vials (approximate height, 2 cm; diameter, 0.8 cm) and large groups of 60 or 150 mites in plastic containers (approximate height, 1.5 cm; diameter, 2 cm diameter or height, 1.5 cm; diameter, 3.5 cm diameter, respectively) all with an approximate 1 cm base of plaster-of-Paris. All plaster-of-Paris bases were completely soaked in water before mites were transferred into them. All mites were reared at a constant 23 °C, at high humidity ( >90%) and were provided an excess of powdered yeast ad libitum.Classifying male morphsTo illustrate the discontinuous distribution of the weapon and to demonstrate that this classification based on visual inspection is non-subjective, we performed phenotypic measurements from male mites from a population collected near Krakόw, Poland, that had previously been fixed onto microscope slides for a separate study66. The measurements taken were idiosoma (body without mouthparts) length and width of third proximal segment of the third right leg (genu). Measurements were preformed using Lecia DM5500B microscope and Lecia Application Suite v.4.6.1. We then performed an analysis to, first, determine whether the allometric relationship between idiosoma length and width of third pair of legs is best described as discontinuous and, second, to verify that classification by simple visual inspection matches the same classification from allometric analysis. One researcher performed all the measurements and classified each male as a fighter (n = 50) or scrambler (n = 50), a separate researcher was then given the measurements but not the classification of the male morph.Broadly, guidelines for the analysis of non-linear allometries92 were followed. The log–log scatterplots of idiosoma length against leg width were visualized, which showed there was clear evidence for non-linear scaling relationships. Next histograms of idiosoma length, leg width and relative leg width (leg width/idiosoma length) were visualized (Extended Data Fig. 2a–c). Where a normal distribution of idiosoma length, and a binomial distribution in leg width and relative leg width are further indications of a discontinuous relationship. On the basis of the lowest point between the two peaks of the density plot of relative leg width (Extended Data Fig. 2c) males were classified as scramblers (relative leg width 0.125). Replotting the log–log scatterplot of idiosoma length and leg width, and using the classification of morph described above clearly demonstrates the discontinuous allometric relationship of idiosoma length and leg width in R. robini (Extended Data Fig. 2d). Moreover, on the basis of the Akaike information criterion (AIC), the discontinuous model where males were assigned a morph (AIC = 646.5) clearly has a substantially better fit than a simple linear and quadratic models (AIC = 918.5 and 920.2, respectively). Further models were omitted from comparison (for example, breakpoint or sigmoidal) due to the clear discontinuous allometry observed. Finally, all 100 males were assigned the same morph by visual inspection and blind allometric analysis, demonstrating that the former is effective and accurate in classifying male morph.Phenotypic assaysFecundity assays were performed using experimental evolution females at F20 and F32. Eggs laid by females between days 4–8 were counted, encompassing the window of time of most evolutionary relevance for female fitness during maintenance of selection lines (that is, egg laying period in selection lines was between days 6–7) and also likely to capture variation in lifetime fecundity that remains largely consistent throughout the first 3 weeks of life93. Nymphs were individually isolated to gain virgin females, which on maturation females from each experimental evolution population (n = 30) were paired with a male from the stock population (15 with fighters and 15 with scramblers). Pairs were transferred to a new vial on day 4, with the pair being removed from the second vial after a further 4 days and all eggs in the second vial counted. If the male had died in the first vial, they were replaced with a stock male of the same morph. Any female deaths in the first or second vials were recorded.Longevity assays were also performed at F20 and F32. At F20, females used in fecundity assays, including the stock male they were paired with (replaced if dead), were transferred to a new vial at day 8. After this point, vials were then checked every 2 days for female deaths and pairs were moved to new vials every 4 days. Males were replaced with stock males of the same morph if found dead. Similarly, at F20, on maturation males from experimental evolution populations (n = 30) were paired with stock females, vials were checked every 2 days and changed every 4 days, with females being replaced if dead. At F32, only female longevity was determined and was performed in groups; 30 experimental evolution females and 30 stock males (15 of each morph) were placed in plastic containers, two per experimental population. This logistically easier estimate of longevity was done due to local restrictions during the SARS-CoV-2 pandemic and the imposed limitations on people working closely together. Groups were checked for dead females every other day and all remaining live mites transferred to a new container every 4 days. When mites were transferred to a new container the sex and morph ratio were balanced to that of the remaining females, by either removing or adding males of the desired morph from the stock population.To determine whether the survival of mites differed between F- and S-lines when competition between males was allowed, at F45 we created small colonies from each population and survival of males and females recorded over 6 days, the same period as used between selecting founders of the next generation and subsequent egg laying period. Colonies were at a 50:50 sex ratio, established with 150 newly eclosed mites placed into small plastic containers. This was approximately the same density after selection of the next generations founders during the maintenance of experimental evolution populations (150 mites in roughly 9.5 cm2 = 16 mites per 1 cm2; 1,000 mites in roughly 67 cm2 = 15 mites per 1 cm2). After 3 days, all colonies were checked and any dead mites identified by sex. After another 3 days, again dead mites were recorded and all surviving mites sexed and counted.Additionally, at F45 we performed further fecundity assays to obtain estimates of inbreeding depression within experimental evolution populations. To establish family groups, larvae were isolated and on maturation F0 males and females (n = 16) from within the same experimental evolution population were paired together. Pairs were allowed to produce eggs for 48 h, after which adults were removed from vials. After hatching from each pair, 12 F1 larvae were isolated into new vials. On their maturation, these F1 mites were either paired with a full sibling, that is, from the same family, or with an individual from a different family but from the same experimental evolution population. When possible, we made two inbred and two outbred pairs with same family lines used. Again, pairs were allowed to produce eggs for 48 h before their removal for the vial. After a further 5 days, vials were checked for larvae, if larvae were present in the first vial six were individually isolated and the second vial discarded, if no larvae were present in the first vial the second vial was checked for larvae and, if present, they were isolated. This protocol therefore produced inbred and outbred individuals from within the same experimental evolution population. Which, as above, on maturation F2 inbred and outbred females were paired with stock males (fighter males only) and number of eggs laid between days 4 and 8 counted. Only a single female from each unique inbred or outbred family was used. Either due to pairs failing to produce offspring or there being no F2 females, samples sizes were not exactly equal. In total, 59 outbred and 55 inbred females from F-lines, and 56 outbred and 54 inbred females from S-lines were paired with stock males.Phenotypic assay statistical analysesAll phenotypic analysis was conducted using R statistical software94 (v.3.5.2) and data were visualized using ggplot2 (ref. 95).Analysis of male morph proportion was performed using a generalized linear mixed model with binomial error structure, fitted using lme4 (ref. 96). Where the proportion of desired morph was compared in model with morph selection and generation (as a factor) including their two-way interaction as explanatory variables, and population included as a random effect.All fecundity data were analysed using generalized linear mixed models with Poisson error structures, fitted using lme4. Due to the differences in stock population males used between F45 and earlier generations, and slightly different rearing conditions between females in the fecundity assays from generations F20 and F32, they were analysed separately from data collected in F45. However, we noted that the fecundity of females in Fig. 5a was comparable to the outbred females in Fig. 5b. Explanatory variables fitted to fecundity data from F20 and F32 were, morph selection treatment, generation, including their two-way interaction term, and stock male morph. The explanatory variables fitted to fecundity data from inbreeding depression data were, morph selection treatment and status of female (that is, inbred or outbred), including their two-way interaction term. In both analyses, we included population as a random effect and an observation level random effect to account for overdispersion, we omitted fitting random slopes due to issues with increasing the complexity of random effects close to reaching a singular fit. Females that died before the end of the fecundity assay and those that laid zero eggs were removed from analysis. This excluded five females from F20 (three F-line and two S-line), 20 from F32 (13 F-line and seven S-line) and 16 from F45 (three inbred and three outbred F-line, and nine inbred and one outbred S-line).Longevities of females at F20 and F32, and males at F20, were analysed separately using mixed effects Cox models, fitted using coxme97. In all analyses, we included a random effect of population, with morph selection treatment as an explanatory variable and extra variable of male morph included in female longevity analysis at F20. Survival of mites over 6 days at F45 was analysed using a GLM with counts of dead and surviving mites fitted with a quasibinomial error structure, the model included morph selection treatment and sex, including their interaction term, as explanatory variables. If individuals were lost due to handling error (that is, killed or escaped) they were right-censored during analysis.Genome assemblySample originA line of R. robini originated from a wild-collected population from the Mosina region (Wielkopolska, Poland). In October 2017, onions were collected from the field and approximately 200 individuals of R. robini were identified under dissecting microscope. The line used for DNA isolation in the genome sequencing project was developed from full sib × sib mating for 14 generations (to maximize homozygosity) following and continuing the protocol described in ref. 67.DNA extractionFor DNA extraction we used only mite eggs, that were laid by 500 females, collected in a container (see above for a description) Females were kept in this container for 3 days. After that time, they were removed, and eggs were filtered using fine sieves and washed for 1 min in 0.3% sodium hypochlorite solution and in Milli-Q water for 2 × 2 min to remove any potential foreign DNA contamination. These eggs were collected in 1.5 ml Eppendorf tube and after short centrifugation, the remains of the water (supernatant) removed with a pipette. The sample was immediately transferred to ice and prepared for DNA extraction. DNA was extracted using Bionano Prep Animal Tissue DNA Kit for HMW DNA isolation according to the manufacturer’s instructions. Briefly, eggs were smashed with a sterile pellet pestle on ice in 500 μl homogenization buffer; the sample was fixed with 500 μl cold ethanol and incubated 60 min on ice, after that time the sample was centrifuged at 1500g for 5 min at 4 °C and the supernatant was discarded. Next, after resuspension in a homogenization buffer pellet, this was cast in four agarose plugs as described in the original protocol. Agarose plugs were incubated with Proteinase K and Lysis buffer solution for 2 h with intermittent mixing. After that time, the digestion solution was replaced with a freshly made one and incubated overnight with intermittent mixing. According to the original protocol, after RNase A digestion and plug washing, DNA was recovered by incubation of the plugs in TE buffer, followed by plug melting and addition of agarase. Recovered DNA was dialysed and homogenized on a membrane for 45 min at room temperature and transferred to a clean tube with a wide bore tip.SequencingSequencing was done using Oxford Nanopore Technologies (ONT, MinION). Isolated DNA purified using AMPure XP beads and resuspended in H2O before library preparation. Two separate libraries were prepared using ligation sequencing kit, SQK-LSK109 and Rapid Sequencing Kit SQK-RAD004, respectively, according to the manufacturer’s protocols and were sequenced on a FLO-MIN106 R9.4.1 SpotON flow cell on a MinION Mk 1B sequencer (ONT). The total yields from sequencing were 484,700 reads (2,417,068,187 nt) with a read-N50 of 10,044 nt (ranging from 216,403 to 100). Base calling of the raw reads was done using Guppy (v.3.3) resulting in a total sum of the reads 7,979,616,172, equivalent to 26× coverage aiming for a genome of 300 megabases (Mb). The reads N50/N90 were estimated at 7,958/1,719.Assembling reference genomeReads aligning with the Mitochondrion genome were identified using BLASTN and filtered from the raw reads before assembling the genome. The remaining ONT reads were assembled using the Flye software (v.2.6), with –min-overlap 3,000 to increase stringency at the initial overlay step, and default parameters including five rounds of polishing through consensus, contigs were additionally polished two times with Medaka (v.0.11.2). Illumina paired-end RNA dataset is assembled using CLC Assembler (CLC Assembly Cell). Both RNA assemblies and paired-end 10X genomic dataset (unpublished data) were mapped onto the contigs using minimap2 (v.2.16) and BWA mapper (v.0.7.17), respectively, and the assembly was further polished using PILON (v.1.20) to error correct potential low-quality regions. The resulting assembly yielded a genome of 307 Mb, assembled into 1,533 contigs ranging from 10,840,357 to 100 basepairs (bp) and an assembly-N50 of 1.670 Mb. Moreover, the BUSCO completeness analysis using the Arachnida (odb10) reference set confirmed our assembly represents the complete genome C:94.8%(S:89.1%,D:5.7%),F:0.9%,M:4.3%,n:2934 (=arachnida_odb10), only missing 126 genes from the whole reference set. Knowing that BUSCO only gives a rough estimation, we remain confident that this assembly represents well the bulb mite genome.Flow-cytometryWhole individual R. robini were homogenized in 500 μl of ice-cold LB01 detergent buffer along with the head of a male Drosophila melanogaster (1 C = 0.18 pg) as an internal standard. The homogenized tissue was filtered through a 30-μm nylon filter. Then 12 μl of propidium iodide with 2 μl of RNase was added, and stained for 1 h on ice in the dark. All samples were run on an FC500 flow cytometer (Beckman-Coulter) using a 488-nm blue laser, providing output as single-parameter histograms showing relative fluorescence between the standard nuclei and the R. robini nuclei. Six replicate samples were run to account for variation in fluorescence outputs. The genome size of R. robini was estimated at 0.30 pg, or about 293 Mb, and consistent with estimation of size from the genome assembly described above.Mitochondrial genomeONT reads aligned with R. robini mitochondrion genome were de novo assembled with Flye (v.2.6) assembler and polished with Racon. Mitochondrion genome is assembled in one single contig with a size of 15,335 bases.Gene predictionOn the polished final genome, protein coding genes have been predicted. For this, AUGUSTUS was used including hints coming from R. robini RNA-sequencing (RNA-seq) (samples SRR3934324, SRR3934325, SRR3934326, SRR3934327, SRR3934328, SRR3934329, SRR3934330, SRR3934331, SRR3934332, SRR3934333, SRR3934335, SRR3934337, SRR3934338 and SRR3934339 from the PRJNA330592 BioProject deposited at the National Center for Biotechnology Information (NCBI) Short Read Archive) and proteins coming from highly curated Tetranychus urticae (v.2020-03-20) as well as proteins from the previous version of the unpublished, Illumina-sequenced R. robini genome (https://public-docs.crg.eu/rguigo/Data/fcamara/bulbmite.v4a/). The PE RNA-seq reads were mapped on the genome using HISAT2 (-k 1 —no-unal) and further processed with Regtools to extract junction hints and filtered for junctions with a minimum coverage of 10. All the RNA-seq reads were also assembled with CLC Assembly Cell (v.5.2.0) software, setting the word size for the Bruijn graph at 50 and maximum bubble size at 31. The reads were assembled into 689,563 contigs (ranging from 10,675 to 180 bp), which were later mapped on the genome with GenomeTheader to generate complementary DNA hints. Protein hints were generated by using with Exonerate (v.2.2) with Protein2Genome model. To reduce the amount of overprediction due to repeated elements (transposable elements, simple sequence repeats) we de novo predicted high abundant repeats using RepeatModeler. The accompanying parameter file for extrinsic data for AUGUSTUS was adapted to include these hints as well as the softmasking of the genomic sequence. The resulting gene predictions from AUGUSTUS were further curated with EvidenceModeler using the same extrinsic data. The BUSCO analysis confirmed that our gene prediction indeed captured the expected genes well (C:94.6%(S:86.3%,D:8.3%),F:0.4%,M:5.0%,n:2934 (=arachnida_odb10)). The final predicted gene set was subsequently processed to be uploaded into ORCAE (https://bioinformatics.psb.ugent.be/orcae/overview/Rhrob)98.ResequencingGenomic sampling and mappingFor genomic analyses we sampled material from each of the morph selection lines (n = 8) at F1, F12 and F29. Following the experimental evolution protocols, after the first 24 h of egg laying all adults were transferred to a new container (described above) for a second 24 h to lay eggs and from these second dishes genomic material was sampled. On maturation, adults were transferred to and kept for 3 days in containers. Adults were then randomly selected and placed into Digestion Solution for MagJET gDNA Kit (F1&12) or ATL buffer (F29) before freezing at −20 °C. From each population two samples were collected consisting of 100 individuals (1 × 100 females and 1 × 100 males of random morph), the two samples separated by sex were used as technical replicates. The tissue from the 100 individuals within each sample was homogenized and DNA was extracted by Proteinase K digestion (24 h) followed by standard procedures using MagJET Genomic DNA Kit (ThermoScientific, F1&12) or DNeasy Blood and Tissue (Qiagen, F29). DNA concentration was controlled with the Qubit double-stranded DNA HS Assay Kit and DNA quality was assessed on agarose gels. The library preparation was performed using NEBNext Ultra II FS DNA Library Prep kit for Illumina.Whole-genome resequencing was carried out by National Genomics Infrastructure (Uppsala, Sweden) using the Illumina Nova-Seq 6000 platform with S4 flow cell to produce 2 × 150 bp reads (average 160.7 × 106: range 130.7 × 106 − 189.9 × 106). Adaptors were trimmed from reads using Trimmomatic99 software (v.0.39) and unpaired reads discarded. Fastq files were mapped to the assembled genome with bwa mem100 (v.0.7.17-r1188) using default settings. Sam files were converted to bam files, sorted, duplicates marked and ambiguously mapped reads removed using samtools101 (v.1.9). On average, 90% (range, 86–93%) of the reads from each sample were mapped successfully, of which an average of 17% (range, 15–19%) were marked as duplicates. This left us with an average of 117.7 × 106 pair end reads per sample, ranging between 99.6 × 106 and 145.9 × 106 (Supplementary Table 1).Genomic analysisFile preparation and filteringPreparation of files used in genomic analysis was done as follows: bam files were converted to a pile-up file using samtools, following which indels and surrounding windows (5 bp either side) were filtered, using identify-genomic-indel-regions.pl and filter-pileup-by-gtf.pl in PoPoolation102 (v.1.2.2) to avoid false SNPs, with the resulting filtered pile-up file converted to a sync file using mpileup2sync.pl in PoPoolation2 (ref. 103) (v.1.201). Using custom python scripts, the distribution of coverage from each sample (single sex) was determined by recording the coverage of positions every 10 kb across the genome from the sync file to give information on expected coverage (Supplementary Fig. 1). On the basis of this, we filtered the sync and pile-up files to contain only regions within a range of informative coverage, where the mean coverage of all samples at every position was between 50% of the expected coverage and 200% of the expected coverage (56×, range 23−112×). The pile-up and sync files containing individual male and female samples (48 in total) were then merged by sex to give files containing allele frequencies from 24 samples (eight populations across three generations), each consisting of allele frequencies of 200 individuals (100 males and 100 females, above) and used in all subsequent analysis (unless stated otherwise). Similarly, we drew coverage of a position every 10 kb from each sample in the sex-merged sync file to determine a distribution from which we decided to subsample to (Supplementary Fig. 1). We putatively identified X-linked contigs (below) and excluded them autosomal analysis. A similar, but, separate analysis on genes and SNPs from X-linked contigs was performed by using different parameters (below).Estimating nucleotide diversityUsing PoPoolation we determined various estimates of genetic diversity per sample (that is, 24 sex-merged samples). The pile-up file from each sample was subsampled using subsample-pileup.pl to a coverage of 63× (max coverage, 252×) to standardize estimations of genetic diversity across the genome, between populations and across generations. First, nucleotide diversity (Tajima’s Pi, π) and number of segregating sites (Watterson’s theta, ϴ) were estimated within genes. We performed analysis of exons using Syn-nonsyn-at-position.pl, in which genetic diversity of synonymous and non-synonymous positions were determined. Further analysis of overall genetic diversity within exons and introns were performed using Variance-at-position.pl, Tajima’s D (D) also estimated in the former. We used a minimum count of three (equal to a minor allele frequency of roughly 5%) for a SNP to be called, and a phred score >30 and a pool size of 400. Further analysis using 10 kb sliding windows (step size 10 kb) across the genome were performed using Variance-sliding.pl, and also included estimation of D. Estimates of D require the minimum count to be 2, but otherwise all the same parameters were used.We filtered genes to be included in our analysis (and all subsequent analysis) on the basis of a number of criteria. On the basis of extensive RNA-seq data from both males and females (Plesnar-Bielak, unpublished data with NCBI accession number PRJNA796800), we only included genes in our analyses that were expressed at a mean level of fragments per kilobase of transcript per million mapped reads >1 across 72 samples originated from both sexes and both morphs rearing in three different temperatures (18, 23 and 28 °C). A further filtering step was performed to remove genes with inconsistent mapping between samples, only genes with >60% exons mapped to (calculated from positions used to calculate parameters in the Syn-nonsyn-at-position.pl π outputs), with 63−252× coverage, in all 24 samples were included in the analysis. The final dataset contained 13,389 autosomal genes and subsequently used to filter other datasets to retain this set of genes only (see Supplementary Table 8 for a list of genes). Similarly, windows were discarded from outputs if 60% of genes being mapped to in all 24 samples) and reducing the final X-linked dataset to contain fewer than 200 genes. We therefore opted to reduce the target coverage further to 40×, in an attempt to retain more genes. This slight reduction of target coverage increased the number of genes in the final dataset substantially to 587 genes. We therefore opted to use a minimum coverage of 40× in all analysis of X-linked SNPs, genes and windows.Diverging SNPsTo determine divergent SNPs between F- and S-lines, we extracted the allele frequencies of all samples from the sex combined sync file. Samples from F29 were then used to filter the entire dataset to only contain SNPs on the basis of a number of criteria. First, positions within all samples were required to have a coverage >63× and 5% (that is, the average of all samples but not necessarily above >5% in all samples). Thus, our dataset contained only positions with the target coverage in all F29 samples and in which polymorphisms were unlikely to be a consequence of sequencing errors. After this filtering we were left with roughly 6 million SNPs used in further analysis. We performed a GLM, at each position by comparing the count of the major allele against counts of minor alleles at F29, to determine consistent allele frequency changes between treatments70. If any population had minor or major allele count of 0, +1 was added to minor and major alleles from all samples. To correct for multiple testing, we converted P values to q values using the qvalues R package (v.2.14.1)104 and applied a FDR with a q 900,000), GLMs were performed (identical to above) on the simulated major and minor allele counts. Using a FDR with a q  More

Abiotic parametersThe temperature regime experienced by the embryos was purposefully natural and therefore varied between the three air exposure treatments. The subtidal treatment, where embryos were continuously submerged in water, remained around 9.5 °C for the duration of the experiment, while the intertidal treatments experienced dips in temperature during outside air exposure down to 2.5 °C and 0.8 °C, for the low and high intertidal respectively (Fig. 3A). Accumulated thermal units (ATU; days × temperature post collection until hatch) for each air exposure treatment were 79.6, 75.4 and 65.3 for subtidal, low intertidal and high intertidal, respectively. Despite differences in thermal regime, peak hatch was on the same day (March 14, 2021) for all air exposure and CO2 treatments, estimated at 11 dpf.Figure 3Temperature and pH experienced by the herring embryos and larvae throughout the experiment. Hourly measurements of (A) air/water temperature experienced by herring embryos in each of the tidal treatments (subtidal: continuous immersion in 9.5 °C water; low intertidal: 2 × 2 h air exposure; high intertidal: 5 + 9 h air exposure) and (B) pH levels in the tanks for each of the CO2 treatments (greens = control, 400 µatm CO2, yellows = medium, 1600 µatm CO2, reds = high, 3000 µatm CO2); dots are pH levels measured in the individual jars during larval incubation.Full size imageThe pH levels in the tanks were measured hourly and were stable over the course of the embryonic incubation period, with no overlap between treatments, although there was some overlap between individual jars. Control treatment was consistently around a pH of 8, the medium treatment had a pH of 7.4 and the high CO2 treatment had a mean pH of 7.1 (Table 1). After hatch, when the larvae were transferred to the jars, circulation and gas exchange between jars and tank were not as high and CO2 accumulated in the jars over time, leading to pH levels deviating from tank pH levels (Fig. 3B). Although oxygen levels remained high (7–9 mg/L), the pH dropped from a mean 8–7.6 in the control on two occasions, and was brought back up with a partial water exchange from the incubation tank water. The pH in the medium and high CO2 treatments were not as affected (Fig. 3B), however, final water chemistry measurements after completion of the experiment (2 days post water exchange) revealed much higher CO2 levels in all treatments (Table 1: day 15).Table 1 Mean water parameters for each treatment (mean of 3 tanks ± S.D.) at the beginning (day 1, 2021-03-06) and end (day 6, 2021-03-12) of embryonic incubation and mean parameters in the jars (N = 9) at the end of larval incubation (day 15, 2021-03-19); Temperature (T), salinity, pCO2, total CO2 (TCO2) measured at distinct sampling intervals with the BoL; total alkalinity (TA) and pH (on the total scale) calculated with CO2SYS.Full size tableEffect of air exposure and CO2 treatment during embryonic developmentNeither embryonic survival nor growth were significantly affected by treatment in our experiment. Percent daily embryonic mortality was low and not significantly affected by CO2 treatment or air exposure (CO2: p = 0.088, F2 = 2.45; Tide: p = 0.11, F2 = 2.19; CO2*Tide: p = 0.18, F2 = 1.59) . Egg diameter at 6 dpf was also not significantly affected by treatment (CO2: p = 0.38, X2 (2, N = 30) = 1.92; Tide: p = 0.83, X2 (2, N = 30) = 0.33; CO2*Tide: p = 0.08, X2 (2, N = 30) = 8.25). Metabolic rate, as indicated by embryonic heart rate, was significantly affected by air exposure at 6 dpf (p  *; 0.1  >).Full size image More