Ecology

Study site and set-upThe studied hemiboreal riparian forest is a 40-year old Filipendula type gray alder (Alnus incana (L.) Moench) forest stand grown on a former agricultural land. It is situated in the Agali Village (58°17′N; 27°17′E) in eastern Estonia within the Lake Peipsi Lowland50 (Supplementary Figs. 12 and 13).The area is characterized by a flat relief with an average elevation of 32 m a.s.l., formed from the bottom of former periglacial lake systems, it is slightly inclined (1%) towards a tributary of the Kalli River. The soil is Gleyic Luvisol. The thickness of the humus layer was 15–20 cm. The content of total carbon (TC), total nitrogen (TN), nitrate (NO3−–N), ammonia NH4+–N, Ca, and Mg per dry matter in 10 cm topsoil was 3.8 and 0.33%, and 2.42, 2.89, 1487 and 283 mg kg−1, respectively, which was correspondingly 6.3, 8.3, 4.4, 3.6, 2.3, and 2.0 times more than those in 20 cm deep zone (Supplementary Table 3).The long-term average annual precipitation of the region is 650 mm, and the average temperature is 17.0 °C in July and –6.7 °C in January. The duration of the growing season is typically 175–180 days from mid-April to October51.The mean height of the forest stand is 17.5 m, stand density 1520 trees per ha, the mean stem diameter at breast height 15.6 cm, basal area 30.5 m2 ha−1, the growing stock 245 m3 ha−1, and the current annual increment of stems 12.0 m3 ha−1 year−1 (based on Uri et al.52 and Becker et al.53). In the forest floor, the following herbs dominate: Filipendula ulmaria (L.) Maxim., Aegopodium podagraria L., Cirsium oleraceum (L.) Scop., Geum rivale L., Crepis paludosa (L.) Moench, shrubs (Rubus idaeus L., Frangula alnus L., Daphne mezereum L.), and young trees (A. incana, Prunus padus (L.)) dominate. In the moss-layer Climacium dendroides (Hedw.) F. Weber & D. Mohr, Plagiomnium spp and Rhytidiadelphus triquetrus (Hedw.) Warnst are overwhelming.Environmental characteristics of hot momentsBased on high emissions of N2O, dynamics of SWC, and near-ground air temperature, we identified four hot moments and related them to soil and environmental variables (see numbers in Fig. 1): wet (1), dry (2) with drought onset (2a), freeze-thaw (3), and dry-minor (4). The main criterion for the hot moments was a rapid increase in N2O emissions of any source.Anomalies from the mean of each hot moment period illustrate the pattern of fluxes during the hot moments (Supplementary Fig. 2). At the end of the freeze-thaw period, the rising SWC is driven by snowmelt became a leading determinant (Supplementary Fig. 2). During the wet period, the rise in soil emissions was accompanied by a remarkable increase in the EC-based ecosystem fluxes. However, all the other hot moments were isolated to soil surfaces.Soil flux measurementsSoil fluxes were measured using 12 automatic dynamic chambers located at 1–2 m distance from each studied tree and installed in June 2017 (Supplementary Fig. 11, see also54). The chambers were made from polymethyl methacrylate (Plexiglas) covered with non-transparent plastic film. Each soil chamber (volume of 0.032 m³) covered a 0.16 m² soil surface. To avoid stratification of gas inside the chamber, air with a constant flow rate of 1.8 L min−1 was circulated within a closed loop between the chamber and gas analyzer unit during the measurements by a diaphragm pump. The air sample was taken from the top of the chamber headspace and pumped back by distributing it to each side of the chamber. For the measurements, the soil chambers were closed automatically for 9 min each. The flushing time of the whole system with ambient air between measurement periods was 1 min. Thus, there were ~12 measurements per chamber per day, making a total of 144 flux measurements per day. A Picarro G2508 (Picarro Inc., Santa Clara, CA, USA) gas analyzer using cavity ring-down spectroscopy (CRDS) technology was used to monitor N2O gas concentrations in the frequency of ~1.17 measurements per second. The chambers were connected to the gas analyzer using a multiplexer allowing a sequent practically continuous measurement.To account for initial stabilization after chamber closing and flushing time, we used 5 min out of the total 9 min closing time (~350 concentration measurements) to estimate slope change of N2O concentration, which was the basis for soil flux calculations.After the quality check, 105,830 flux values (98.7% of total possible) of soil N2O fluxes could be used during the whole study period.Stem flux measurementsThe tree stem fluxes were measured manually with frequency 1–2 times per week from September 2017 until December 2018. Twelve representative mature gray alder trees were selected for stem flux measurements and equipped with static closed tree stem chamber systems for stem flux measurements20. Soil fluxes were investigated close to each selected tree. The tree chambers were installed in June 2017 in the following order: at the bottom part of the tree stem (~10 cm above the soil) and at 80 and 170 cm above the ground. The rectangular shape stem chambers were made of transparent plastic containers, including removable airtight lids (Lock & Lock Co Ltd, Seoul, Republic of Korea). For the chamber, preparation see Schindler et al.54. Two chambers per profile were set randomly across 180° and interconnected with tubes into one system (total volume of 0.00119 m³) covering 0.0108 m² of stem surface. A pump (model 1410VD, 12 V; Thomas GmbH, Fürstenfeldbruck, Germany) was used to homogenize the gas concentration prior to sampling. Chamber systems remained open between each sampling campaign. During 60 measurement campaigns, four gas samples (each 25 ml) were collected from each chamber system via septum in a 60 min interval: 0/60/120/180 min sequence (sampling time between 12:00 and 16:00) and stored in pre-evacuated (0.3 bar) 12 ml coated gas-tight vials (LabCo International, Ceregidion, UK). The gas samples were analyzed in the laboratory at the University of Tartu within a week using gas chromatography (GC-2014; Shimadzu, Kyoto, Japan) equipped with an electron capture detector for detection of N2O and a flame ionization detector for CH4. The gas samples were injected automatically using Loftfield autosampler (Loftfield Analytics, Göttingen, Germany). For gas-chromatographical settings see Soosaar et al.55.Soil and stem flux calculationFluxes were quantified on a linear approach according to the change of CH4 and N2O concentrations in the chamber headspace over time, using the equation according to Livingston and Hutchison56.Stem fluxes were quantified on a linear approach according to the change of N2O concentrations in the chamber headspace over time. A data quality control (QC) was applied based on R2 values of linear fit for CO2 measurements. When the R2 value for CO2 efflux was above 0.9, the conditions inside the chamber were applicable, and the calculations for N2O gases were also accepted in spite of their R2 values.To compare the contribution of soil and stems, the stem fluxes were upscaled to hectares of ground area based on average stem diameter, tree height, stem surface area, tree density, and stand basal area estimated for each period. A cylindric shape of the tree stem was assumed. To estimate average stem emissions per tree, fitted regression curves for different periods were made between the stem emissions and height of the measurements as previously done by Schindler et al.54.EC instrumentationEC system was installed on a 21 m height scaffolding tower. Fast 3-D sonic anemometer Gill HS-50 (Gill Instruments Ltd., Lymington, Hampshire, UK) was used to obtain three wind components. CO2 fluxes were measured using the Li-Cor 7200 analyser (Li-Cor Inc., Lincoln, NE, USA). Air was sampled synchronously with the 30 m teflon inlet tube and analyzed by a quantum cascade laser absorption spectrometer (QCLAS) (Aerodyne Research Inc., Billerica, MA, USA) for N2O concentrations. The Aerodyne QCLAS was installed in the heated and ventilated cottage near the tower base. A high-capacity free scroll vacuum pump (Agilent, Santa Clara, CA, USA) guaranteed an airflow rate of 15 L min−1 between the tower and gas analyzer during the measurements. Air was filtered for dust and condense water. All measurements were done at 10 Hz and the gas-analyzer reported concentrations per dry air (dry mixing ratios).Eddy-covariance flux calculation and data QCThe fluxes of N2O were calculated using the EddyPro software (v.6.0-7.0, Li-Cor) as a covariance of the gas mixing ratio with the vertical wind component over 30-min periods. Despiking of the raw data was performed following Mauder et al.57. Anemometer tilt was corrected with the double-axis rotation. Linear detrending was chosen over block averaging to minimize the influence of possible fluctuations of a gas analyzer. Time lags were detected using covariance maximization in a given time window (5 ± 2 s was chosen based on the tube length and flow rate). While Webb-Pearman-Leuning (WPL) correction58 is typically performed for the closed-path systems, we did not apply it as water correction was already performed by the Aerodyne and the software reported dry mixing ratios. Both low and high-frequency spectral corrections were applied using fully analytic corrections59,60.Calculated fluxes were filtered out in case they were coming from the half-hour averaging periods with at least one of the following criteria: more than 1000 spikes, half-hourly averaged mixing ratio out of range (300–350 ppb), QC flags higher than 761.The footprint area was estimated using Kljun et al.62 implemented in TOVI software (Li-Cor Inc.). A footprint allocation tool was implemented to flag the non-forested areas within the 90% cumulative footprint and fluxes appointed to these areas were removed from the further analysis.Storage fluxes were estimated using concentration measurements from the eddy system (Eq. (1)), assuming the uniform change within the air column under the tower during every 30 min period63,64:$${mathrm{S}} = {Delta}{mathrm{C}}/{Delta}{mathrm{t}} ast {mathrm{z}}_{mathrm{m}},$$
(1)
where S is storage, ΔC is change in the dry mixing ratio of N2O, Δt is time period (30 min), zm is measurement height (21 m).In the absence of a better estimate or profile measurements, these estimates were used to correct for storage change. Total flux values that were higher than eight times the standard deviation were additionally filtered out (following Wang et al.36). Overall, the QC procedures resulted in 61% data coverage.While friction velocity (u*) threshold is used to filter eddy fluxes of CO265, visual inspection of the friction velocity influence on N2O fluxes demonstrated no effect. Thus, we decided not to apply it, taking into account that the 1–9 QC flag system already marks the times when the turbulence is not sufficient.To obtain the continuous time-series and to enable the comparison to chamber estimates over hourly time scales, gap-filling of N2O fluxes was performed using marginal distribution sampling method implemented in ReddyProcWeb online tool (https://www.bgc-jena.mpg.de/bgi/index.php/Services/REddyProcWeb) (described in detail in Wutzler et al.66).MATLAB (ver. 2018a-b, Mathworks Inc., Natick, MA, USA) was used for all the eddy fluxes data analysis.Ancillary measurementsAir temperature, relative and absolute humidity were measured within the canopy at 10 m height using the HC2A-S3—Standard Meteo Probe/RS24T (Rotronic AG, Bassersdorf, Switzerland) and Campbell CR100 data logger (Campbell Scientific Inc., Logan, UT, USA). The potential amount of dissolved N2O in the atmospheric water was calculated based on the absolute humidity and the maximum solubility of N2O in water67. DPD was calculated from air temperature and estimated dew point temperature to characterize the chance of fog formation within the canopy. The solar radiation data were obtained from the SMEAR Estonia station located at 2 km from the study site68 using the Delta-T-SPN-1 sunshine pyranometer (Delta-T Devices Ltd., Cambridge, UK). The cloudiness ratio was calculated as the ratio of diffuse solar radiation to total solar radiation.Near-ground air temperature, soil temperature (Campbell Scientific Inc.), and SWC sensors (ML3 ThetaProbe, Delta-T Devices, Burwell, Cambridge, UK) were installed directly on the ground and 0–10 cm soil depth close to the studied tree spots. During six campaigns from August to November 2017, composite topsoil samples were taken with a soil corer from a depth of 0–10 cm for physical and chemical analysis using standard methods69.Statistical analysisR version 4.0.2 (R Development Core Team, 2020) was used to examine, analyze and visualize the data. The significance level (alpha) considered for all the tests was 0.05. The “akima” package version 0.6–2.1 was used to create interpolated contour plots representing a three-dimensional surface70 by plotting soil temperature and SWC against soil N2O emissions as the independent variable. Linear regression models were fitted and Spearman’s rank correlation coefficients were shown for change of SWC and soil N2O flux in period drought onset and air temperature and soil N2O flux in period freeze-thaw. Spearman’s rank correlation coefficients were also shown characterizing the relationship between the monthly average number of days with a high chance of sunshine and fog formation and the difference between the N2O flux from soil and ecosystem. Regarding all measurements of soil temperature, SWC, and soil N2O flux, relationships were better represented by nonlinear than linear models. In addition, the Bragg equation with four parameters71 was used for describing the relationship between SWC and soil N2O flux in the dry period. A workflow for the nonlinear regression analysis was used72 and regression models were fitted in R using functions lm, nls, or loess. More

Our study addressed the interaction between the SASL and the operation of the Chinook salmon fishery in the Toltén River, supported by a non-native and invasive species found populating many rivers and adjacent sea of Chile6. Our results demonstrate that the interaction between SASLs and the small-scale fishing communities of Caleta La Barra vary depending on different factors, both operational, such as the number of boats, and environmental, such as moon luminosity during the fishing operations. We were able to determine common patterns among fishers that allow us to establish different social profiles that shape their relationship with the SASL and also their relationship with the natural environment in which they live.Despite the frequent presence of SASL during fishing operations ( > 90%), we assessed that only 35% of them constituted interactions. This frequency of interaction is lower than recorded by other studies in this species18,20,21, but higher than in others17. These differences indicate that there is a high level of variability in the frequency of interactions between the SASL and the artisanal fisheries, both at spatial20,21 and temporal17,18 scales. According to de María et al.40, these differences can be explained by different factors such as the season in which the study was carried out, the productivity of each area, the fishing gear used and the captured species, among others. As the Chinook salmon is only fished during the austral summer months (January and February), it was not possible to analyze potential temporal variations in the frequency of interactions in this study.A positive relationship was observed between the number of SASL and the number of boats operating. A similar relationship was reported by Goetz et al.19, who observed groups of SASLs following boats during fishing operations. This can be explained by a number of factors. First, Szteren and Paéz15 reported that SASL in Uruguay are able to recognize the sound emitted by the boats during a fishing operation, and that is an indicator of prey availability (“dinner bell” effect). Second, it is possible that some individuals have learned to feed during fishing operations, specializing in a type of prey or in a particular feeding strategy associated with fishing operations41. Likewise, due to the general decrease in resource availability due to overexploitation in Chile42, SASLs could be attracted to fishing activities, associating a higher number of boats with greater availability and easy-to-capture prey19,43.We also found a positive relationship between the number of SASLs and the moon’s luminosity. Nights with high lunar luminosity (full moon) are associated with the largest tidal ranges of the month. These tides are related to a high productivity due to the formation of tidal fronts characterized by the abrupt difference in temperature, oxygen and fluorescence, commonly observed during spring tides, increasing the concentration of zooplankton and therefore attracting more predators44. This could result in greater resource availability. Largest tidal ranges may be also linked to influx of returning Chinook salmon in higher numbers than the rest of the month (authors’ unpublished results). However, it cannot be ruled out that the better visibility provided by the increased luminosity would have allowed the observer to count a greater number of SASLs around the boat16,19.Although non-significant, we found that the interactions between SASLs and the fishing operation occurred more often when fishing nets were deeper and the total catch was higher. A relationship between the occurrence of interactions and depth could be related to the time the net is underwater. At greater depths the time of hauling increases, which in turn increase the opportunity for the SASL to interact with the boat. This trend is opposite to what was found by other authors in other species45 who observed more interaction close to the surface as a learned strategy to reduce energetic demands. However, in our study the maximum depth of a fishing net was 10 m which is much lower than the depth SASLs usually dive during their foraging trips46. With respect to the relationship between the occurrence of interactions and total catch, different studies have demonstrated that fishing operations and seals and sea lions coincide in areas where the resource is more abundant15,47, and thus the number of SASL raiding the nets increases when there are more fish caught in the net and easily accessible.We found no relationship between CPUE and the interaction of SASL during fishing operations, i.e. no variations in the standardized catch of Chinook salmon per haul were recorded regardless of interactions with SASL. Similar results have been identified by other authors in studies of gillnet and purse-seine fisheries in other areas, both in Chile and elsewhere for this same species18,19,26,40. This lack of relationship could be explained by a number of factors. Firstly, the number of fish consumed by SASLs in the fishing gear is not high enough to generate differences in the CPUE at the fishery scale15,18,20. Secondly, and as mentioned before, SASLs and fishing operations coincide in areas where the resource is more abundant. Therefore, and even if an interaction with the SASL was reported, the CPUE was maintained or was higher in areas where SASLs are present, in comparison to areas where this predator is absent15. Finally, the volume of fish catches is affected by additional factors besides the presence of SASLs, such as environmental conditions, the presence of other predators, and the abundance of resources in each area, and therefore do not exclusively depend on interactions with the SASL40.South American sea lions are frequently blamed for causing significant impacts to the economy of local fishers26. However, at the scale of the whole fishery, damage to catch from SASLs interactions was recorded in only four of 22 fishing events and in 2,5% of the total catch, suggesting that damage is smaller than perceived by fishers. It is important to note, however, that the assessment of damage by SASLs in this study is conservative, since we do not consider fish that could be wholly removed from the net, which can increase the total biomass lost by SASLs. Competition for fishing resources between SASLs and artisanal fishers has been largely documented in different areas and situations and will likely exacerbate negative perception about this mammal and an obstacle to fishing operations17,24. In a recent study, Oliveira et al.26 demonstrated that the actual economic losses caused by sea lions to the local fishery in Brazil are much smaller (0.8–3% of the productivity of the monitored boats for the analyzed year) compared to the large damage perceived by the fishers. This perception extends to other cases of human-big marine carnivore interactions and competition during fishing activities48,49. Therefore, it is possible that the observed negative perception about SASLs is independent of a particular event, but it rather relates to the continuous presence of SASLs in the area, and also to the general perception that the SASLs can eat hundreds of fish over the course of one interaction26.However, it is important to note that economic losses could be relevant on a boat-to-boat basis, and this clearly contributes to exacerbate a negative perception towards SASLs. In the case of the Chinook salmon fishery from Toltén River, we quickly appraised that SASL interactions may result in significant losses on a boat-to-boat basis. Ex-vessel price of Chinook salmon has been on average CLP$3500 per kg (USD$5; currency rate at 27 April 2021), suggesting that a 10–15 kg salmon may be worth CLP$35,000–53,000 (USD$50–76). Boats with a damaged catch by SASL lost up to 11% of their revenue during the duration of the survey. Considering fishers’ sole reliance on a short fishing season8 and the overexploitation of alternative, native fishery resources42, fishers’ perceptions are likely to be negative towards SASL given these economic losses.The negative perception of SASLs by fishers needs to be understood within the context in which nature has been historically and socially altered. Since Chinook salmon was successfully introduced around 25 years ago50, both human fishing practices and the behavior of SASLs have changed in important ways. Fishing has moved from a year-round communitarian activity performed in the open sea, focused on abundant small native marine species, to a mostly familiar activity, performed in the estuary during a certain period of the year and focused on a big and profitable catch. The SASLs have also learned to feed on this new species and changed their predatory practices, moving from the sea to the river, and predominantly predating on salmon in the short period in which they arrive in the estuary. There is a mutual coproduction31,32,33 process in which these three parties—salmon, fisher, and SASLs—have each modified their behaviors and condition of reproduction and existence. The salmon colonized a new habitat, fishers learned to fish this new attractive species, and the SASL also learned to prey on this new species, colonizing a new space in the estuary, following fishing boats. This coproduction has presented new opportunities for the human population, as they have been able to exploit a new economic resource that has revitalized the economy of the town, but it has also refueled a long-standing conflict between humans and SASLs. Furthermore, these changes in fishing activity may have increased the perceived interactions with SASLs and the damage they cause. Moving fishing operations from the sea to the smaller and calmer estuarine waters makes the sea lion more visible. Also, a fishing based on a big (~ 15–20 kg each) and profitable Chinook salmon species, makes any sea lion attack more damaging than former fishing practices in which SASLs took one or two small fish from big net hauls.In order to address the negative perception of SASLs amongst artisanal fishers, we framed the discourse about SASLs within a broader discourse with nature using the concept of valuation language from Martinez-Alier51. Using Likert scales we observed that, contrary to what was expected and seen in other contexts, fishers do not demonstrate a view in line with the “environmentalism of the poor”. However, they did possess a view that values sustainability, though not in a conservationist context that would be expected from a community that rely on natural resources for their livelihood, but it is not purely instrumental either. This leads us to rethink the Martinez-Alier51 categories in order to describe the combination that was found among fishers in this study: a mix between elements of utilitarian rationality, viewing fish solely as a resource, and elements of a deep caring for the river and sea life, linking local economics with environmental wellbeing. However, this view did not include some of the more substantive and political aspects of an “environmentalism of the poor” perception, such as the right for all natural beings to exist. This discourse was better described as part of the “popular sustainable development” view, as the sustainable development paradigm is based more on an economic, rather than ecological, rationality that transforms ‘nature’ into ‘environment’ and resources52. However, in this case, caring for the river and sea, and economic wellbeing were seen as mutually necessary in a discourse that is accompanied by the respect of natural forces as being important to the life and wellbeing of the community.We also found that most fishers, from all socio-demographic characteristics, hold this “popular sustainable development” view. This was followed by instrumental rationality (4.0), held mostly by older fishers with lower levels of formal education. These results are consistent with those found in other studies24,53,54, in which older fishers have the most negative attitude and perception of SASLs because of the several negative encounters with these species24. This commonality of the popular sustainable development view, accompanied by a more conservationist view held by younger fishers with more years of formal education, opens up the possibility for a better relationship among Chinook salmon, fishers, and SASLs, and a better coexistence. We hypothesize that it is possible to make local effective governance changes to improve this coexistence55, for example, by changing the fishing practices to modify the SASL behavior. Currently, fishers process the fish and dispose of the waste in the river where they fish. In their own words this practice “domesticate” SASLs, and the community is making agreements to move this practice away from the shore. It also seems that SASLs follow the boats as they have learned that boats mean availability of prey. This requires improved governance systems to redistribute fishing places within the river, as fishers already partially and informally do it. Also, it is recommendable the implementation of a long-term education program to the fishing community that include, among other issues, the critical role of this marine top predator on the trophic webs and consequently the negative impacts of their removal18,24. Actions like these could help to move from a paradigm of “defending the fishing from the SASLs” to a better local understanding of the relationship between Chinook salmon, SASLs, and human behavior. In recent years, due to the important economic opportunity that Chinook salmon has meant for the community, young fishers with more years of formal education, and a more sustainable and conservationist view of fishing, have returned to La Barra to occupy leadership roles in the community. We cannot be sure how this leadership will evolve as they grow older and more experienced, but the formal education and experience backgrounds they hold separate them from their parents and may anticipate a different trajectory. We speculate that this is a step in the right direction at resolving SASL-fisheries conflicts. More

1.Reddin, C. J., Kocsis, Á. T., Aberhan, M. & Kiessling, W. Victims of ancient hyperthermal events herald the fates of marine clades and traits under global warming. Glob. Chang. Biol. 27, 868–878 (2021).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
2.Reddin, C. J., Kocsis, Á. T. & Kiessling, W. Marine invertebrate migrations trace climate change over 450 million years. Glob. Ecol. Biogeogr. 27, 704–713 (2018).Article 

Google Scholar 
3.Kordas, R. L., Harley, C. D. G. & O’Connor, M. I. Community ecology in a warming world: The influence of temperature on interspecific interactions in marine systems. J. Exp. Mar. Bio. Ecol. 400, 218–226 (2011).Article 

Google Scholar 
4.Poloczanska, E. S. et al. Responses of marine organisms to climate change across oceans. Front. Mar. Sci. 3, 1–21 (2016).Article 

Google Scholar 
5.Hanken, J. & Wake, D. B. Miniaturization of body size: Organismal consequences and evolutionary significance. Annu. Rev. Ecol. Syst. 24, 501–519 (1993).Article 

Google Scholar 
6.Sheridan, J. A. & Bickford, D. Shrinking body size as an ecological response to climate change. Nat. Clim. Chang. 1, 401–406 (2011).ADS 
Article 

Google Scholar 
7.Forster, J., Hirst, A. G. & Atkinson, D. Warming-induced reductions in body size are greater in aquatic than terrestrial species. Proc. Natl. Acad. Sci. U. S. A. 109, 19310–19314 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
8.Ohlberger, J. Climate warming and ectotherm body size – from individual physiology to community ecology. Funct. Ecol. 27, 991–1001 (2013).Article 

Google Scholar 
9.Garilli, V. et al. Physiological advantages of dwarfing in surviving extinctions in high-CO 2 oceans. Nat. Clim. Chang. 5, 678–682 (2015).ADS 
CAS 
Article 

Google Scholar 
10.Daufresne, M., Lengfellner, K. & Sommer, U. Global warming benefits the small in aquatic ecosystems. Proc. Natl. Acad. Sci. U. S. A. https://doi.org/10.1073/pnas.0902080106 (2009).Article 
PubMed 
PubMed Central 

Google Scholar 
11.Verberk, W. C. E. P. et al. Shrinking body sizes in response to warming: explanations for the temperature–size rule with special emphasis on the role of oxygen. Biol. Rev. 96, 247–268 (2021).PubMed 
Article 
PubMed Central 

Google Scholar 
12.Angilletta, M. J., Steury, T. D. & Sears, M. W. Temperature, growth rate, and body size in ectotherms: Fitting pieces of a life-history puzzle. Integr. Comp. Biol. 44, 498–509 (2004).PubMed 
Article 
PubMed Central 

Google Scholar 
13.Hoving, H. J. T. et al. Extreme plasticity in life-history strategy allows a migratory predator (jumbo squid) to cope with a changing climate. Glob. Chang. Biol. 19, 2089–2103 (2013).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
14.Angilletta, M. J. & Dunham, A. E. The temperature-size rule in ectotherms: simple evolutionary explanations may not be general. Am. Nat. 162, 332–342 (2003).PubMed 
Article 
PubMed Central 

Google Scholar 
15.Vinarski, M. V. On the applicability of Bergmann’s rule to ectotherms: the state of the art. Biol. Bull. Rev. 4, 232–242 (2014).Article 

Google Scholar 
16.Atkinson, D. Temperature and organism size: a biological law for organisms?. Adv. Ecol. Res. 25, 1–58 (1994).Article 

Google Scholar 
17.Atkinson, D. Effects of temperature on the size of aquatic ectotherms: Exceptions to the general rule. J. Therm. Biol. 20, 61–74 (1995).Article 

Google Scholar 
18.Forster, J. & Hirst, A. G. The temperature-size rule emerges from ontogenetic differences between growth and development rates. Funct. Ecol. 26, 483–492 (2012).Article 

Google Scholar 
19.Gillooly, J. F., Brown, J. H., West, G. B., Savage, V. M. & Charnov, E. L. Effects of size and temperature on metabolic rate. Science 293, 2248–2251 (2001).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
20.Schulte, P. M. The effects of temperature on aerobic metabolism: Towards a mechanistic understanding of the responses of ectotherms to a changing environment. J. Exp. Biol. 218, 1856–1866 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
21.Gillooly, J. F., Charnov, E. L., West, G. B., Savage, V. M. & Brown, J. H. Effects of size and temperature on developmental time. Nature 417, 70–73 (2002).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
22.Riemer, K., Anderson-Teixeira, K. J., Smith, F. A., Harris, D. J. & Ernest, S. K. M. Body size shifts influence effects of increasing temperatures on ectotherm metabolism. Glob. Ecol. Biogeogr. 27, 958–967 (2018).Article 

Google Scholar 
23.Rosa, R. et al. Ocean warming enhances malformations, premature hatching, metabolic suppression and oxidative stress in the early life stages of a keystone squid. PLoS ONE 7, e38282 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
24.Pecl, G. T. & Jackson, G. D. The potential impacts of climate change on inshore squid: Biology, ecology and fisheries. Rev. Fish Biol. Fish. 18, 373–385 (2008).Article 

Google Scholar 
25.Twitchett, R. J. The Lilliput effect in the aftermath of the end-Permian extinction event. Palaeogeogr. Palaeoclimatol. Palaeoecol. 252, 132–144 (2007).Article 

Google Scholar 
26.Harries, P. J. & Knorr, P. O. What does the ‘Lilliput Effect’ mean?. Palaeogeogr. Palaeoclimatol. Palaeoecol. 284, 4–10 (2009).Article 

Google Scholar 
27.Metcalfe, B., Twitchett, R. J. & Price-Lloyd, N. Changes in size and growth rate of ‘Lilliput’ animals in the earliest Triassic. Palaeogeogr. Palaeoclimatol. Palaeoecol. 308, 171–180 (2011).Article 

Google Scholar 
28.Chu, D. et al. Lilliput effect in freshwater ostracods during the Permian-Triassic extinction. Palaeogeogr. Palaeoclimatol. Palaeoecol. 435, 38–52 (2015).Article 

Google Scholar 
29.Urbanek, A. Biotic crises in the history of upper silurian graptoloids: a palaeobiological model. Hist. Biol. https://doi.org/10.1080/10292389309380442 (1993).Article 

Google Scholar 
30.Urlichs, M. Stunting in invertebrates from the type area of the Cassian Formation (Early Carnian) of the dolomites (Italy). GeoAlp 8, 164–169 (2011).
Google Scholar 
31.Morten, S. D. & Twitchett, R. J. Fluctuations in the body size of marine invertebrates through the Pliensbachian-Toarcian extinction event. Palaeogeogr. Palaeoclimatol. Palaeoecol. https://doi.org/10.1016/j.palaeo.2009.08.023 (2009).Article 

Google Scholar 
32.Piazza, V., Ullmann, C. V. & Aberhan, M. Temperature-related body size change of marine benthic macroinvertebrates across the Early Toarcian Anoxic Event. Sci. Rep. 10, 1–13 (2020).Article 
CAS 

Google Scholar 
33.Calosi, P., Putnam, H. M., Twitchett, R. J. & Vermandele, F. Marine metazoan modern mass extinction: improving predictions by integrating fossil, modern, and physiological data. Ann. Rev. Mar. Sci. 11, 369–390 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
34.Gerber, S. Comparing the differential filling of morphospace and allometric space through time: the morphological and developmental dynamics of Early Jurassic ammonoids. Paleobiology 37, 369–382 (2011).Article 

Google Scholar 
35.Pálfy, J. & Smith, P. L. Synchrony between Early Jurassic extinction, oceanic anoxic event, and the Karoo-Ferrar flood basalt volcanism. Geology 28, 747–750 (2000).ADS 
Article 

Google Scholar 
36.Caruthers, A. H., Smith, P. L. & Gröcke, D. R. The Pliensbachian-Toarcian (Early Jurassic) extinction, a global multi-phased event. Palaeogeogr. Palaeoclimatol. Palaeoecol. 386, 104–118 (2013).Article 

Google Scholar 
37.Wignall, P. B. Large igneous provinces and mass extinctions. Earth Sci. Rev. 53, 1–33 (2001).ADS 
CAS 
Article 

Google Scholar 
38.Percival, L. M. E. et al. Globally enhanced mercury deposition during the end-Pliensbachian extinction and Toarcian OAE: A link to the Karoo-Ferrar Large Igneous Province. Earth Planet. Sci. Lett. 428, 267–280 (2015).ADS 
CAS 
Article 

Google Scholar 
39.Foster, G. L., Hull, P., Lunt, D. J. & Zachos, J. C. Placing our current ‘hyperthermal’ in the context of rapid climate change in our geological past. Phil. Trans. R. Soc. A 376, 20170086 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
40.Suan, G. et al. Secular environmental precursors to Early Toarcian (Jurassic) extreme climate changes. Earth Planet. Sci. Lett. 290, 448–458 (2010).ADS 
CAS 
Article 

Google Scholar 
41.Fantasia, A. et al. Global versus local processes during the Pliensbachian-Toarcian transition at the Peniche GSSP, Portugal: A multi-proxy record. Earth Sci. Rev. 198, 102932 (2019).CAS 
Article 

Google Scholar 
42.Müller, T. et al. Ocean acidification during the early Toarcian extinction event: Evidence from Boron isotopes in brachiopods. Geology 48, 1184–1188 (2020).ADS 
Article 

Google Scholar 
43.Suan, G., Mattioli, E., Pittet, B., Mailliot, S. & Lécuyer, C. Evidence for major environmental perturbation prior to and during the Toarcian (Early Jurassic) oceanic anoxic event from the Lusitanian Basin Portugal. Paleoceanography 23, A1202 (2008).ADS 
Article 

Google Scholar 
44.Dera, G. et al. High-resolution dynamics of early Jurassic marine extinctions: The case of Pliensbachian-Toarcian ammonites (Cephalopoda). J. Geol. Soc. London 167, 21–33 (2010).CAS 
Article 

Google Scholar 
45.Dera, G. et al. Climatic ups and downs in a disturbed Jurassic world. Geology 39, 215–218 (2011).ADS 
CAS 
Article 

Google Scholar 
46.Miguez-Salas, O., Rodríguez-Tovar, F. J. & Duarte, L. V. Selective incidence of the Toarcian oceanic anoxic event on macroinvertebrate marine communities: a case from the Lusitanian basin Portugal. Lethaia 50, 548–560 (2017).Article 

Google Scholar 
47.Correia, V. F., Riding, J. B., Duarte, L. V., Fernandes, P. & Pereira, Z. The palynological response to the Toarcian Oceanic Anoxic Event (Early Jurassic) at Peniche, Lusitanian Basin, western Portugal. Mar. Micropaleontol. 137, 46–63. https://doi.org/10.1016/j.marmicro.2017.10.004 (2017).ADS 
Article 

Google Scholar 
48.Rita, P., Nätscher, P., Duarte, L. V., Weis, R. & De Baets, K. Mechanisms and drivers of belemnite body-size dynamics across the Pliensbachian-Toarcian crisis. R. Soc. Open Sci. 6, 190494 (2019).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
49.Caswell, B. A. & Coe, A. L. The impact of anoxia on pelagic macrofauna during the Toarcian Oceanic Anoxic Event (Early Jurassic). Proc. Geol. Assoc. 125(4), 383–391. https://doi.org/10.1016/j.pgeola.2014.06.001 (2014).Article 

Google Scholar 
50.Ullmann, C. V., Thibault, N., Ruhl, M., Hesselbo, S. P. & Korte, C. Effect of a Jurassic oceanic anoxic event on belemnite ecology and evolution. Proc. Natl. Acad. Sci. U. S. A. 111, 10073–10076 (2014).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
51.Dera, G., Toumoulin, A. & de Baets, K. Diversity and morphological evolution of Jurassic belemnites from South Germany. Palaeogeogr. Palaeoclimatol. Palaeoecol. 457, 80–97 (2016).Article 

Google Scholar 
52.Neige, P., Weis, R. & Fara, E. Ups and downs of belemnite diversity in the Early Jurassic of Western Tethys. Palaeontology 64, 263–283 (2021).Article 

Google Scholar 
53.Rita, P., De Baets, K. & Schlott, M. Rostrum size differences between Toarcian belemnite battlefields. Foss. Rec. 21, 171–182 (2018).Article 

Google Scholar 
54.Rita, P. et al. Biogeographic patterns of belemnite body size responses to episodes of environmental crisis. PeerJ Prepr. (2019).55.Hoffmann, R. & Stevens, K. The palaeobiology of belemnites – foundation for the interpretation of rostrum geochemistry. Biol. Rev. 95, 94–123 (2020).Article 

Google Scholar 
56.Adams, D. C. & Otárola-Castillo, E. Geomorph: An r package for the collection and analysis of geometric morphometric shape data. Methods Ecol. Evol. 4, 393–399 (2013).Article 

Google Scholar 
57.Schlegelmilch, R. Die Belemniten des süddeutschen Jura. Die Belemniten des süddeutschen Jura https://doi.org/10.1007/978-3-8274-3083-0 (1998).Article 

Google Scholar 
58.McArthur, J. M. et al. Sr-isotope stratigraphy (87Sr/86Sr) of the lowermost Toarcian of Peniche, Portugal, and its relation to ammonite zonations. Newsletters Stratigr. 53, 297–312 (2020).Article 

Google Scholar 
59.Hesselbo, S. P., Jenkyns, H. C., Duarte, L. V. & Oliveira, L. C. V. Carbon-isotope record of the Early Jurassic (Toarcian) Oceanic Anoxic Event from fossil wood and marine carbonate (Lusitanian Basin, Portugal). Earth Planet. Sci. Lett. 253, 455–470 (2007).ADS 
CAS 
Article 

Google Scholar 
60.Klug, C., Schweigert, G., Fuchs, D., Kruta, I. & Tischlinger, H. Adaptations to squid-style high-speed swimming in Jurassic belemnitids. Biol. Lett. 12, 20150877 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
61.Seibel, B. A., Thuesen, E. V., Childress, J. J. & Gorodezky, L. A. Decline in pelagic cephalopod metabolism with habitat depth reflects differences in locomotory efficiency. Biol. Bull. 192, 262–278 (1997).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
62.Mattioli, E., Pittet, B., Petitpierre, L. & Mailliot, S. Dramatic decrease of pelagic carbonate production by nannoplankton across the Early Toarcian anoxic event (T-OAE). Glob. Planet. Change 65, 134–145 (2009).ADS 
Article 

Google Scholar 
63.Chamberlain, J. A. Locomotion in ancient seas: Constraint and opportunity in Cephalopod adaptive design. Geobios 15, 49–61 (1993).Article 

Google Scholar 
64.Rexfort, A. & Mutterlose, J. The role of biogeography and ecology on the isotope signature of cuttlefishes (Cephalopoda, Sepiidae) and the impact on belemnite studies. Palaeogeogr. Palaeoclimatol. Palaeoecol. 284, 153–163 (2009).Article 

Google Scholar 
65.Holland, S. M. The quality of the fossil record: a sequence stratigraphic perspective. Paleobiology 26, 148–168 (2000).Article 

Google Scholar 
66.Holland, S. M. The non-uniformity of fossil preservation. Philos. Trans. R. Soc. B Biol. Sci. 371, 2 (2016).
Google Scholar 
67.Korn, D. Impact of enviornmental perturbations o heterochronic develpments in Palaeozoic ammonoids. Evol. Chang. Heterochrony 245–260 (1995).68.Yacobucci, M. M. Plasticity of developmental timing as the underlying cause of high speciation rates in ammonoids. in Advancing research on living and fossil cephalopods 59–76 (Springer, Boston, MA, 1999).69.Landman, N. H. & Gyssant, J. R. Heterochrony and ecology in Jurassic and Cretaceous ammonites. Geobios 26, 247–255 (1993).Article 

Google Scholar 
70.McNamara, K. J. Heterochrony: the evolution of development. Evol. Educ. Outreach 5, 203–218 (2012).Article 

Google Scholar 
71.Dahlke, F. T., Wohlrab, S., Butzin, M. & Pörtner, H. O. Thermal bottlenecks in the life cycle define climate vulnerability of fish. Science 369, 65–70 (2020).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
72.Pörtner, H. O. & Farrell, A. P. Ecology: Physiology and climate change. Science https://doi.org/10.1126/science.1163156 (2008).Article 
PubMed 
PubMed Central 

Google Scholar 
73.Pimentel, M. S. et al. Impact of ocean warming on the early ontogeny of cephalopods: A metabolic approach. Mar. Biol. 159, 2051–2059 (2012).Article 

Google Scholar 
74.Komoroske, L. M. et al. Ontogeny influences sensitivity to climate change stressors in an endangered fish. Conserv. Physiol. 2, 1–13 (2014).Article 
CAS 

Google Scholar 
75.Pörtner, H. O., Bock, C. & Mark, F. C. Oxygen- & capacity-limited thermal tolerance: Bridging ecology & physiology. J. Exp. Biol. 220, 2685–2696 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
76.Harnik, P. G., Simpson, C. & Payne, J. L. Long-term differences in extinction risk among the seven forms of rarity. Proc. R. Soc. B Biol. Sci. 279, 4969–4976 (2012).Article 

Google Scholar 
77.Reddin, C. J., Kocsis, Á. T. & Kiessling, W. Climate change and the latitudinal selectivity of ancient marine extinctions. Paleobiology 45, 70–84 (2019).Article 

Google Scholar 
78.Dorey, N. et al. Ocean acidification and temperature rise: Effects on calcification during early development of the cuttlefish Sepia officinalis. Mar. Biol. 160, 2007–2022 (2013).CAS 
Article 

Google Scholar 
79.Sigwart, J. D. et al. Elevated pCO2 drives lower growth and yet increased calcification in the early life history of the cuttlefish Sepia officinalis (Mollusca: Cephalopoda) Julia. ICES J. Mar. Sci. 73, 970–980 (2016).Article 

Google Scholar 
80.Gutowska, M. A., Melzner, F., Pörtner, H. O. & Meier, S. Cuttlebone calcification increases during exposure to elevated seawater pCO2 in the cephalopod Sepia officinalis. Mar. Biol. 157, 1653–1663 (2010).CAS 
Article 

Google Scholar 
81.Kaplan, M. B., Mooney, T. A., McCorkle, D. C. & Cohen, A. L. Adverse effects of ocean acidification on early development of squid (Doryteuthis pealeii). PLoS ONE 8, e63714 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
82.Otjacques, E. et al. Cuttlefish buoyancy in response to food availability and ocean acidification. Biology (Basel). https://doi.org/10.3390/biology9070147 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
83.Neige, P. & Boletzky, S. Morphometrics of the shell of three Sepia species (Mollusca: Cephalopoda): Intra- and interspecific variation. Zool. Beitraege. 38, 137–156 (1997).
Google Scholar 
84.Rita, P., Weis, R., Duarte, L. V. & De Baets, K. Taxonomical diversity and palaeobiogeographical affinity of belemnites from the Pliensbachian-Toarcian GSSP (Lusitanian Basin, Portugal). Pap. Palaeontol. https://doi.org/10.1002/spp2.1343 (2020).Article 

Google Scholar 
85.MacArthur, R. H. Geographical ecology: patterns in the distribution of species. (Princeton University Press, 1972).86.Gaston, K. J. The structure and dynamics of geographic ranges. (Oxford University Press on Demand, 2003).87.Duarte, L. Sequence stratigraphy and depositional setting of the Pliensbachian and Toarcian marly limestones in the Lusitanian Basin Portugal. Ciências da Terra 16, 17–23 (2007).
Google Scholar 
88.Reddin, C. J., Nätscher, P. S., Kocsis, Á. T., Pörtner, H.-O. & Kiessling, W. Marine clade sensitivities to climate change conform across timescales. Nat. Clim. Chang. 10, 249–253 (2020).ADS 
Article 

Google Scholar 
89.Doyle, P. & Kelly, R. A. The Jurassic and Cretaceous belemnites of Kong Karls Land, Svalbard. (Norsk Polarinstitutt Oslo, 1988).90.Doyle, P. New records of dimitobelid belemnites from the cretaceous of james ross island Antarctica. Alcheringa 14, 159–175 (1990).Article 

Google Scholar 
91.Fedorov, A. et al. 3D Slicer as an image computing platform for the Quantitative Imaging Network. Magn. Reson. Imaging https://doi.org/10.1016/j.mri.2012.05.001 (2012).Article 
PubMed 
PubMed Central 

Google Scholar 
92.Rohatgi, A., Rehberg, S. & Stanojevic, Z. Webplotdigitizer: Version 4.1 of Webplotdigitizer. (2018).https://doi.org/10.5281/zenodo.1137880.93.Plate, T. & Heiberger, R. Package ‘ abind ’. (2016).94.Collyer, M. L. & Adams, D. C. RRPP: An r package for fitting linear models to high-dimensional data using residual randomization. Methods Ecol. Evol. 9, 1772–1779 (2018).Article 

Google Scholar 
95.Sherratt, E. Quick Guide to Geomorph v. 2.0. public.iastate.edu (2014).96.Torchiano, M. Package ‘ effsize ’. (2020).97.Gotelli, N. J., Dorazio, R. M., Ellison, A. M. & Grossman, G. D. Detecting temporal trends in species assemblages with bootstrapping procedures and hierarchical models. Philos. Trans. 365, 3621–3631 (2010).Article 

Google Scholar 
98.Hervé, M. Package ‘ RVAideMemoire ’. (2021).99.Mangiafico, S. Package ‘ rcompanion ’. (2021).100.Oksanen, J. et al. Package ‘vegan’ title community ecology package. Commun. Ecol. Packag. 2, 1–297 (2019).MathSciNet 

Google Scholar 
101.Pinheiro, J. et al. Package ‘nlme’. (2021).102.Dormann, C. F. et al. Collinearity: A review of methods to deal with it and a simulation study evaluating their performance. Ecography (Cop.) 36, 27–46 (2013).Article 

Google Scholar 
103.R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria (2019).104.Garnier, S., Ross, N., Rudis, B. & Sciaini, M. Package ‘viridis’. (2021). More

These authors contributed equally: Claire R. Peart, Sergio Tusso, Saurabh D. Pophaly.Science of Life Laboratories and Department of Evolutionary Biology, Uppsala University, Uppsala, SwedenClaire R. Peart, Sergio Tusso, Chi-Chih Wu, Aaron B. A. Shafer & Jochen B. W. WolfDivision of Evolutionary Biology, Faculty of Biology, LMU Munich, Planegg-Martinsried, GermanyClaire R. Peart, Sergio Tusso, Saurabh D. Pophaly, Fidel Botero-Castro & Jochen B. W. WolfMax Planck Institute for Plant Breeding Research, Cologne, GermanySaurabh D. PophalyLaboratorio de Ecología de Pinnípedos ‘Burney J. Le Boeuf’, Centro Interdisciplinario de Ciencias Marinas, Instituto Politécnico Nacional, Baja California Sur, MéxicoDavid Aurioles-GamboaDepartment of Natural Sciences, University of Houston-Downtown, Houston, TX, USAAmy B. BairdDepartment of Ecology and Conservation Biology, Texas A&M University, College Station, TX, USAJohn W. BickhamBritish Antarctic Survey, Natural Environment Research Council, Cambridge, UKJaume Forcada & Joseph I. HoffmanElephant Seal Research Group, Sea Lion Island, Falkland IslandsFilippo Galimberti & Simona SanvitoDepartment of Anatomy, University of Otago, Dunedin, New ZealandNeil J. GemmellDepartment of Animal Behaviour, Bielefeld University, Bielefeld, GermanyJoseph I. HoffmanNorwegian Polar Institute, Fram Centre, Tromsø, NorwayKit M. Kovacs & Christian LydersenDepartment of Environmental and Biological Sciences, University of Eastern Finland, Joensuu, FinlandMervi Kunnasranta & Tommi NymanNatural Resources Institute Finland (Luke), Joensuu, FinlandMervi KunnasrantaDepartment of Ecosystems in the Barents Region, Norwegian Institute of Bioeconomy Research, Svanhovd Research Station, Svanvik, NorwayTommi NymanLaboratory of Mammal Ecology, Universidade do Vale do Rio dos Sinos, São Leopoldo, BrazilLarissa Rosa de OliveiraNational Oceanic and Atmospheric Administration, National Marine Fisheries Service, Alaska Fisheries Science Center, Marine Mammal Laboratory, Seattle, WA, USAAnthony J. OrrThe Saimaa Ringed Seal Genome Project, Institute of Biotechnology, University of Helsinki, Helsinki, FinlandMia ValtonenForensic Science & Environmental Life Sciences, Trent University, Peterborough, Ontario, CanadaAaron B. A. Shafer More

1.De Grave, S. & Fransen, C. H. J. M. Carideorum Catalogus: The recent species of the dendrobranchiate, stenopodidean, procarididean and caridean shrimps (Crustacea: Decapoda). Zool. Med. Leiden 85(9), 195–588 (2011).
Google Scholar 
2.Shen, H., Braband, A. & Scholtz, G. Mitogenomic analysis of decapod crustacean phylogeny corroborates traditional views on their relationships. Mol. Phylogenet. Evol. 66(3), 776–789 (2013).PubMed 
Article 

Google Scholar 
3.Chace, F. A. Jr. & Kensley, B. The cardiac notch in decapods. J. Crustacean Biol. 12(3), 442–447 (1992).Article 

Google Scholar 
4.Holthuis, L.B., Fransen, C.H.J.M. & Van Achterberg, C. The recent genera of the Caridean and Stenopodidean shrimps (Crustacea, Decapoda): with an appendix on the order Amphionidacea. Nationaal Natuurhistorisch Museum, Leiden, pp. 6–328 (1993).5.Martin, J. W. & Davis, G. E. An updated classification of the recent Crustacea. Los Angeles: Natural History Museum of County. Contrib. Sci. 39, 1–124 (2001).
Google Scholar 
6.De Grave, S. et al. A classification of living and fossil genera of decapod crustaceans. Raffles Bull. Zool. Suppl. 21, 1–109 (2009).
Google Scholar 
7.Bracken, H. D., De Grave, S. & Felder, D. L. Phylogeny of the infraorder Caridea based on mitochondrial and nuclear genes (Crustacea: Decapoda). Decapod Crustacean Phylogenet. 18, 274–298 (2009).
Google Scholar 
8.Li, C. P., De Grave, S., Lei, H. C., Chan, T. Y. & Chu, K. H. Molecular systematics of caridean shrimps based on five nuclear genes: Implications for superfamily classification. Zool. Anz. 250, 270–279 (2011).Article 

Google Scholar 
9.De Grave, S., Li, C. P., Tsang, L. M., Chu, K. H. & Chan, T. Y. Unweaving hippolytoid systematics (Crustacea, Decapoda, Hippolytidae): Resurrection of several families. Zool. Scr. 43(5), 496–507 (2014).Article 

Google Scholar 
10.Baeza, J. A. Protandric simultaneous hermaphroditism in the shrimps Lysmata bahia and Lysmata intermedia. Invertebr. Biol. 127(2), 181–188 (2008).Article 

Google Scholar 
11.Baeza, J. A. & Bauer, R. T. Experimental test of social mediation of sex change in a protandric sequential hermaphrodite; the marine shrimp Lysmata wurdemanni (Crustacea: Caridea). Behav. Ecol. Sociobiol. 55, 544–550 (2004).Article 

Google Scholar 
12.Xu, Y., Song, L. S. & Li, X. Z. The molecular phylogeny of Caridea based on 16S rDNA sequences. Mar. Sci. 29(9), 36–41 (2005).CAS 

Google Scholar 
13.Baeza, J. A. Testing three models on the adaptive significance of protandric simultaneous hermaphroditism in a marine shrimp. Evolution 60, 1840–1850 (2006).PubMed 
Article 
PubMed Central 

Google Scholar 
14.Baeza, J. A. Molecular systematics of peppermint and cleaner shrimps: phylogeny and taxonomy of the genera Lysmata and Exhippolysmata (Crustacea: Caridea: Hippolytidae). Zool. J. Linn. Soc. Lond. 160, 254–265 (2010).Article 

Google Scholar 
15.Baeza, J. A. Molecular phylogeny of broken-back shrimps (genus Lysmata and allies): A test of the ‘Tomlinson-Ghiselin’ hypothesis explaining the evolution of hermaphroditism. Mol. Phylogenet. Evol. 69, 46–62 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
16.Baeza, J. A. & Prakash, S. An integrative taxonomic and phylogenetic approach reveals a complex of cryptic species in the ‘peppermint’ shrimp Lysmata wurdemanni sensu stricto. Zool. J. Linn. Soc. Lond. 185(4), 1018–1038 (2019).Article 

Google Scholar 
17.Boore, J. L. Animal mitochondrial genomes. Nucleic Acids Res. 27(8), 1767–1780 (1999).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
18.Liu, Q. N., Zhu, B. J., Dai, L. S., Wei, G. Q. & Liu, C. L. The complete mitochondrial genome of the wild silkworm moth, Actias selene. Gene 505(2), 291–299 (2012).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
19.Tan, M. H. et al. Comparative mitogenomics of the Decapoda reveals evolutionary heterogeneity in architecture and composition. Sci. Rep. 9, 10756 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
20.Gissi, C., Iannelli, F. & Pesole, G. Evolution of the mitochondrial genome of Metazoa as exemplified by comparison of congeneric species. Heredity 101(4), 301–320 (2008).CAS 
PubMed 
Article 

Google Scholar 
21.Rich, A. & Rajbhandary, U. L. Transfer RNA: Molecular structure, sequence, and properties. Annu. Rev. Biochem. 45(1), 805–860 (1976).CAS 
PubMed 
Article 

Google Scholar 
22.Staton, J. L., Daehler, L. L. & Brown, W. M. Mitochondrial gene arrangement of the horseshoe crab Limulus polyphemus L.: Conservation of major features among arthropod classes. Mol. Biol. Evol. 14(8), 867–874 (1997).CAS 
PubMed 
Article 

Google Scholar 
23.Miller, A. D., Murphy, N. P., Burridge, C. P. & Austin, C. M. Complete mitochondrial DNA sequences of the decapod crustaceans Pseudocarcinus gigas (Menippidae) and Macrobrachium rosenbergii (Palaemonidae). Mar. Biotechnol. 7(4), 339–349 (2005).CAS 
Article 

Google Scholar 
24.Ivey, J. L. & Santos, S. R. The complete mitochondrial genome of the Hawaiian anchialine shrimp Halocaridina rubra Holthuis, 1963 (Crustacea: Decapoda: Atyidae). Gene 394(1–2), 35–44 (2007).CAS 
PubMed 
Article 

Google Scholar 
25.Wang, Z. et al. Complete mitochondrial genome of Parasesarma affine (Brachyura: Sesarmidae): Gene rearrangements in Sesarmidae and phylogenetic analysis of the Brachyura. Int. J. Biol. Macromol. 118, 31–40 (2018).CAS 
PubMed 
Article 

Google Scholar 
26.Lü, Z. M. et al. Complete mitochondrial genome of Ophichthus brevicaudatus reveals novel gene order and phylogenetic relationships of Anguilliformes. Int. J. Biol. Macromol. 135, 609–618 (2019).PubMed 
Article 
CAS 

Google Scholar 
27.Shen, X., Li, X., Sha, Z., Yan, B. & Xu, Q. Complete mitochondrial genome of the Japanese snapping shrimp Alpheus japonicus (Crustacea: Decapoda: Caridea): Gene rearrangement and phylogeny within Caridea. Sci. China Life Sci. 55(7), 591–598 (2012).PubMed 
Article 

Google Scholar 
28.Wang, Q. et al. Characterization and comparison of the mitochondrial genomes from two Alpheidae species and insights into the phylogeny of Caridea. Genomics 112(1), 65–70 (2020).CAS 
PubMed 
Article 

Google Scholar 
29.Shen, X. et al. The complete mitochondrial genome of the ridgetail white prawn Exopalaemon carinicauda Holthuis, 1950 (Crustacean: Decapoda: Palaemonidae) revealed a novel rearrangement of tRNA genes. Gene 437(1–2), 1–8 (2009).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
30.Moritz, C. & Brown, W. M. Tandem duplications in animal mitochondrial DNAs: Variation in incidence and gene content among lizards. Proc. Natl. Acad. Sci. 84(20), 7183–7187 (1987).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
31.Lavrov, D. V., Boore, J. L. & Brown, W. M. Complete mtDNA sequences of two millipedes suggest a new model for mitochondrial gene rearrangements: Duplication and nonrandom loss. Mol. Biol. Evol. 19(2), 163–169 (2002).CAS 
PubMed 
Article 

Google Scholar 
32.Lunt, D. H. & Hyman, B. C. Animal mitochondrial DNA recombination. Nature 387(6630), 247 (1997).ADS 
CAS 
PubMed 
Article 

Google Scholar 
33.Dowton, M. & Campbell, N. J. H. Intramitochondrial recombination-is it why some mitochondrial genes sleep around?. Trends Ecol. Evol. 16(6), 269–271 (2001).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
34.Sun, S., Cheng, J., Sun, S. & Sha, Z. Complete mitochondrial genomes of two deep-sea pandalid shrimps, Heterocarpus ensifer and Bitias brevis: Insights into the phylogenetic position of Pandalidae (Decapoda:Caridea). J. Oceanol. Limnol. 38(3), 816–825 (2020).CAS 
Article 

Google Scholar 
35.Tan, M. H., Gan, H. M., Lee, Y. P., Poore, G. C. & Austin, C. M. Digging deeper: new gene order rearrangements and distinct patterns of codons usage in mitochondrial genomes among shrimps from the Axiidea, Gebiidea and Caridea (Crustacea: Decapoda). Peer J. 5, e2982 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
36.Aljanabi, S. M. & Martinez, I. Universal and rapid salt-extraction of high quality genomic DNA for PCR-based techniques. Nucleic Acids Res. 25, 4692–4693 (1997).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
37.Dierckxsens, N., Mardulyn, P. & Smits, G. NOVOPlasty: De novo assembly of organelle genomes from whole genome data. Nucleic Acids Res. 45(4), e18 (2017).PubMed 

Google Scholar 
38.Bernt, M. et al. MITOS: Improved de novo metazoan mitochondrial genome annotation. Mol. Phylogenet. Evol. 69, 313–319 (2013).Article 

Google Scholar 
39.Altschul, S. F. et al. Gapped BLAST and PSI-BLAST: A new generation of protein database search programs. Nucleic Acids Res. 25(17), 3389–3402 (1997).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
40.Kumar, S., Stecher, G. & Tamura, K. MEGA7: Molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol. Biol. Evol. 33, 1870–1874 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
41.Perna, N. T. & Kocher, T. D. Patterns of nucleotide composition at fourfold degenerate sites of animal mitochondrial genomes. J. Mol. Evol. 41, 353–358 (1995).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
42.Grant, J. R. & Stothard, P. The CGView Server: A comparative genomics tool for circular genomes. Nucleic Acids Res. 36, 181–184 (2008).Article 
CAS 

Google Scholar 
43.Librado, P. & Rozas, J. DnaSP v5: A software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452 (2009).CAS 
Article 

Google Scholar 
44.Nielsen, R. Statistical tests of selective neutrality in the age of genomics. Heredity 86, 641–647 (2001).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
45.Katoh, K., Misawa, K., Kuma, K. & Miyata, T. MAFFT: A novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 30(14), 3059–3066 (2002).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
46.Castresana, J. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol. Biol. Evol. 17(4), 540–552 (2000).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
47.Xia, X. & Xie, Z. DAMBE: Software package for data analysis in molecular biology and evolution. J. Hered. 92(4), 371–373 (2001).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
48.Nguyen, L. T., Schmidt, H. A., von Haeseler, A. & Minh, B. Q. IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol. Biol. Evol. 32(1), 268–274 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
49.Kalyaanamoorthy, S., Minh, B. Q., Wong, T. K. F., von Haeseler, A. & Jermiin, L. S. ModelFinder: Fast model selection for accurate phylogenetic estimates. Nat. Methods 14, 587–589 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
50.Ronquist, F. & Huelsenbeck, J. P. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574 (2003).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
51.Nylander, J. A., Ronquist, F., Huelsenbeck, J. P. & Nieves-Aldrey, J. Bayesian phylogenetic analysis of combined data. Syst. Biol. 53(1), 47–67 (2004).PubMed 
Article 
PubMed Central 

Google Scholar  More

Case dataMedically attended influenza B cases in New Zealand were identified from samples taken from patients with influenza-like illness (ILI) attended by a network of general practitioners recruited for surveillance (2430 cases with an identified influenza B lineage) and from non-surveillance hospital samples (1606 cases with an identified lineage) analyzed by regional diagnostic laboratories and by the World Health Organization (WHO) National Influenza Centre at the Institute for Environmental Science and Research (ESR). Briefly, general practice surveillance operates from May to September, with participating practices collecting nasopharyngeal or throat swabs from the first ILI patient examined on each Monday, Tuesday, and Wednesday. ILI is defined as an “acute respiratory tract infection characterized by an abrupt onset of at least two of the following: fever, chills, headache, and myalgia”38. A subset of the New Zealand data (cases from 2002 to 2013) was previously compiled by Vijaykrishna et al.28 along with cases from Australia reported to the WHO Collaborating Centre for Reference and Research on Influenza in Melbourne.Statistical model of influenza B susceptibility based on infection historyFor lineage V (B/Victoria), we modeled the number of cases in people born in birth year b observed in season y as a multinomial draw with probabilities given by:$${theta }_{V}(b,y)=D(b,y)beta (b,y){Z}_{V}(b,y)rho (b,y)$$
(1)
with an analogous equation defining the multinomial distribution θY(b, y) for lineage Y (B/Yamagata). D(b, y) is the fraction of the population that was born in year b as of observation season y. Z(b, y) is the susceptibility to lineage V during season y of a person born in year b relative to that of an unexposed person. β(b, y) is a baseline probability of infection with influenza B that captures differences in transmission associated with age (thus depending on b and y) and is equal to β1 if people born in year b are in preschool during season y (0–5 years old), β2 if they are school-age children or teenagers (6–17 years old), or β3 if they are 18 or older. ρ(b,y) is an age-specific factor equal to a parameter ρ if people are  0. Letting α1 and α2 be the instantaneous attack rates for preschoolers and school-age children:$${P}_{mathrm{N}}(A)=left{begin{array}{ll}{e}^{-{alpha }_{1}(A-m)}frac{(1-{e}^{-{alpha }_{1}})}{{alpha }_{1}},hfill&,{text{if}},A; le; {A}_{mathrm{s}}\ {e}^{-{alpha }_{1}({A}_{mathrm{s}}-m)}{e}^{-{alpha }_{2}(A-{A}_{mathrm{s}})}frac{(1-{e}^{-{alpha }_{2}})}{{alpha }_{2}},&,{text{if}},A; > ; {A}_{mathrm{s}}end{array}right.$$
(26)
where As is the age at which children start going to school (4 years old in the Netherlands69). It is noteworthy that for school-age children (the equation for A  > As on the bottom), the correction term for uncertainty in sampling is not necessary for the time spent in preschool (assumed to be exactly As years), only for the time after preschool (A − As).Handling cases with missing lineage informationWe assumed cases with missing lineage information in 2002 (n = 61), 2011 (n = 312), and 2019 (n = 206) belonged to B/Victoria, as 99% or more of identified cases in those seasons were B/Victoria (86/87 cases in 2002, 276/280 cases in 2011, and 552/552 cases in 2019) as were 94%, 92%, and 92%, respectively, of isolates from sequence databases (for Australia and New Zealand combined). We assumed cases with missing lineage information belonged to B/Yamagata in 2013 (n = 37), 2014 (n = 77), and 2017 (n = 87), when the majority of identified cases were B/Yamagata (268/272, 131/138, and 473/489, respectively), as were 99%, 94%, and 84%, respectively, of isolates in sequence databases. Unidentified cases in other seasons were disregarded, because both lineages were present at higher frequencies among identified cases. Removing unidentified cases altogether in all seasons led to similar parameter estimates.Sequence divergence analysisTo estimate the amount of evolution within and between lineages, we analyzed all complete HA and NA sequences from human influenza B isolates available on GISAID in July 2019. The set of isolates used in this analysis differs from the set used to estimate lineage frequencies, because we required isolates to have complete sequences (although not all sequences listed as complete on GISAID were in fact complete). Two isolates collected in 2000 (B/Hong Kong/548/2000 and B/Victoria/504/2000) were deposited as B/Victoria but our BLAST assignment indicated they were in fact B/Yamagata (their low divergence from B/Yamagata strains was a clear outlier). NA sequences from isolates B/Kanagawa/73 and B/Ann Arbor/1994 were only small fragments (99 and 100 amino acids long) poorly aligned with other sequences and were thus excluded. We also excluded NA sequences from B/Yamagata isolates B/Catalonia/NSVH100773835/2018 and B/Catalonia/NSVH100750997/2018, because they were extremely diverged (60% and 38%) from the reference strain B/Yamagata/16/88 and aligned poorly with other sequences.To compare sequence diversity within and between lineages over time, we aligned sequences using MAFFT v. 7.31070 and calculated percent amino acid differences in pairs of sequences from the same lineage and in pairs with one sequence from each. For each year, we sampled 100 sequences from each lineage (or used all sequences if 100 or fewer were available) to limit the number of pairwise calculations. To estimate how much B/Yamagata and B/Victoria evolved since the late 1980s, we calculated percent amino acid differences between each B/Yamagata and B/Victoria sequence, and the corresponding HA and NA sequences of reference strains B/Yamagata/16/88 and B/Victoria/2/87. Unlike in the analysis of pairwise divergence within each time point, we used all sequences from each lineage in each year. We excluded sites in which one or both sequences had gaps or ambiguous amino acids.To compare HA and NA divergence between influenza B lineages with divergence between influenza A subtypes, we downloaded complete HA and NA sequences from H3N2 and H1N1 isolated since 1977 and available on GISAID in August 2019. Homologous sites in the HA of H3N2 and H1N1 are difficult to identify by conventional sequence alignment, and instead we used the algorithm by Burke and Smith71 implemented on the Influenza Research Database website72. Both H3N2 and H1N1 sequences were aligned with the reference H3N2 sequence A/Aichi/2/68. We verified that this method matched sites on the stalk and head of the H1N1 HA with sites on the stalk and head of H3N2 HA by comparing the resulting alignment with the alignment in Supplementary Fig. S2 of Kirkpatrick et al.73. To limit the total number of influenza A sequences analyzed, we randomly selected 100 H3N2 and 100 H1N1 sequences for years in which more than 100 sequences were available, and used all available sequences for the remaining years. Isolates A/Canterbury/58/2000, A/Canterbury/87/2000, and A/Canterbury/55/2000 were excluded, because both H1N1-like and H3N2-like sequences were available under the same isolate name on GISAID.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

1.IPBES. Global Assessment Report on Biodiversity and Ecosystem Services of the Intergovernmental Science—Policy Platform on Biodiversity and Ecosystem Services (eds Brondizio, E. S. et al.) (IPBES Secretariat, 2019).2.Millennium Ecosystem Assessment. Ecosystems and Human Well-Being: Biodiversity Synthesis ed Ma (World Resources Institute, 2005). http://www.loc.gov/catdir/toc/ecip0512/2005013229.html. Accessed June 2019.3.Hooper, D. U. et al. A global synthesis reveals biodiversity loss as a major driver of ecosystem change. Nature 486(7401), 105 (2012).ADS 
CAS 
Article 

Google Scholar 
4.Vellend, M. et al. Homogenization of forest plant communities and weakening of species–environment relationships via agricultural land use. J. Ecol. 95(3), 565–573. https://doi.org/10.1111/j.1365-2745.2007.01233.x (2007).Article 

Google Scholar 
5.Karp, D. S. et al. Intensive agriculture erodes β-diversity at large scales. Ecol. Lett. 15(9), 963–970. https://doi.org/10.1111/j.1461-0248.2012.01815.x (2012).Article 
PubMed 

Google Scholar 
6.Anderson, M. J. et al. Navigating the multiple meanings of β diversity: A roadmap for the practicing ecologist. Ecol. Lett. 14(1), 19–28. https://doi.org/10.1111/j.1461-0248.2010.01552.x (2011).ADS 
Article 
PubMed 

Google Scholar 
7.Socolar, J. B., Gilroy, J. J., Kunin, W. E. & Edwards, D. P. How should beta-diversity inform biodiversity conservation?. Trends Ecol. Evol. 31(1), 67–80. https://doi.org/10.1016/j.tree.2015.11.005 (2016).Article 
PubMed 

Google Scholar 
8.Mori, A. S., Isbell, F. & Seidl, R. β-Diversity, community assembly, and ecosystem functioning. Trends Ecol. Evol. 33(7), 549–564 (2018).Article 

Google Scholar 
9.Vellend, M. Conceptual synthesis in community ecology. Q. Rev. Biol. 85(2), 183–206 (2010).Article 

Google Scholar 
10.Wang, S., Lamy, T., Hallett, L. M. & Loreau, M. Stability and synchrony across ecological hierarchies in heterogeneous metacommunities: Linking theory to data. Ecography (Cop) 42(6), 1200–1211. https://doi.org/10.1111/ecog.04290 (2019).Article 

Google Scholar 
11.Olden, J. D. Biotic homogenization: A new research agenda for conservation biogeography. J. Biogeogr. 33(12), 2027–2039. https://doi.org/10.1111/j.1365-2699.2006.01572.x (2006).Article 

Google Scholar 
12.Loreau, M., Mouquet, N. & Gonzalez, A. Biodiversity as spatial insurance in heterogeneous landscapes. Proc. Natl. Acad. Sci. 100(22), 12765–12770 (2003).ADS 
CAS 
Article 

Google Scholar 
13.Harrison, S. Species Diversity, Spatial Scale, and Global Change (Sinauer Sunderland, 1993).
Google Scholar 
14.Sax, D. F. & Gaines, S. D. Species diversity: From global decreases to local increases. Trends Ecol. Evol. 18(11), 561–566 (2003).Article 

Google Scholar 
15.Hillebrand, H. & Matthiessen, B. Biodiversity in a complex world: Consolidation and progress in functional biodiversity research. Ecol. Lett. 12(12), 1405–1419 (2009).Article 

Google Scholar 
16.Magurran, A. E. & McGill, B. J. Biological Diversity: Frontiers in Measurement and Assessment (Oxford University Press, 2010).
Google Scholar 
17.Usseglio, P. Quantifying reef fishes: Bias in observational approaches. In Ecology of Fishes on Coral Reefs (ed Mora, C.) 270–273 (Cambridge University Press, 2015). https://www.cambridge.org/core/books/ecology-of-fishes-on-coral-reefs/quantifying-reef-fishes-bias-in-observational-approaches/660760F9E62CC61DEB48C8124AD44CDC. Accessed June 2019.18.Caldwell, Z. R., Zgliczynski, B. J., Williams, G. J. & Sandin, S. A. Reef Fish survey techniques: Assessing the potential for standardizing methodologies. PLoS One 11(4), e0153066. https://doi.org/10.1371/journal.pone.0153066 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
19.Worm, B. et al. Impacts of biodiversity loss on ocean ecosystem services. Science 314(5800), 787–790 (2006).ADS 
CAS 
Article 

Google Scholar 
20.Barbier, E. B. Marine ecosystem services. Curr. Biol. 27(11), R507–R510 (2017).CAS 
Article 

Google Scholar 
21.Goodwin, K. D. et al. DNA sequencing as a tool to monitor marine ecological status. Front. Mar. Sci. 4, 107. https://doi.org/10.3389/fmars.2017.00107 (2017).Article 

Google Scholar 
22.Deiner, K. et al. Environmental DNA metabarcoding: Transforming how we survey animal and plant communities. Mol. Ecol. 26(21), 5872–5895. https://doi.org/10.1111/mec.14350 (2017).Article 
PubMed 

Google Scholar 
23.Taberlet, P., Coissac, E., Pompanon, F., Brochmann, C. & Willerslev, E. Towards next-generation biodiversity assessment using DNA metabarcoding. Mol. Ecol. 21(8), 2045–2050. https://doi.org/10.1111/j.1365-294X.2012.05470.x (2012).CAS 
Article 
PubMed 

Google Scholar 
24.Creer, S. et al. The ecologist’s field guide to sequence-based identification of biodiversity. Methods Ecol. Evol. 7(9), 1008–1018. https://doi.org/10.1111/2041-210X.12574 (2016).Article 

Google Scholar 
25.Stat, M. et al. Ecosystem biomonitoring with eDNA: Metabarcoding across the tree of life in a tropical marine environment. Sci. Rep. 7(1), 12240. https://doi.org/10.1038/s41598-017-12501-5 (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
26.Bakker, J. et al. Environmental DNA reveals tropical shark diversity in contrasting levels of anthropogenic impact. Sci. Rep. 7(1), 16886. https://doi.org/10.1038/s41598-017-17150-2 (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
27.Port, J. A. et al. Assessing vertebrate biodiversity in a kelp forest ecosystem using environmental DNA. Mol. Ecol. 25(2), 527–541. https://doi.org/10.1111/mec.13481 (2016).CAS 
Article 
PubMed 

Google Scholar 
28.Andruszkiewicz, E. A. et al. Biomonitoring of marine vertebrates in Monterey Bay using eDNA metabarcoding. PLoS One 12(4), e0176343. https://doi.org/10.1371/journal.pone.0176343 (2017).MathSciNet 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
29.Yamamoto, S. et al. Environmental DNA metabarcoding reveals local fish communities in a species-rich coastal sea. Sci. Rep. 7, 40368. https://doi.org/10.1038/srep40368 (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
30.O’Donnell, J. L. et al. Spatial distribution of environmental DNA in a nearshore marine habitat. PeerJ 5, e3044. https://doi.org/10.7717/peerj.3044 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
31.Jeunen, G.-J. et al. Environmental DNA (eDNA) metabarcoding reveals strong discrimination among diverse marine habitats connected by water movement. Mol. Ecol. Resour. 19(2), 426–438. https://doi.org/10.1111/1755-0998.12982 (2019).CAS 
Article 
PubMed 

Google Scholar 
32.Stat, M. et al. Combined use of eDNA metabarcoding and video surveillance for the assessment of fish biodiversity. Conserv. Biol. 33(1), 196–205 (2019).Article 

Google Scholar 
33.West, K. M. et al. eDNA metabarcoding survey reveals fine-scale coral reef community variation across a remote, tropical island ecosystem. Mol. Ecol. 29(6), 1069–1086. https://doi.org/10.1111/mec.15382 (2020).CAS 
Article 
PubMed 

Google Scholar 
34.Graham, H. M. Effects of local deforestation on the diversity and structure of Southern California giant kelp forest food webs. Ecosystems 7(4), 341–357. https://doi.org/10.1007/s10021-003-0245-6 (2004).Article 

Google Scholar 
35.Miller, R. J. et al. Giant kelp, Macrocystis pyrifera, increases faunal diversity through physical engineering. Proc R Soc B Biol Sci 285(1874), 20172571 (2018).Article 

Google Scholar 
36.Lamy, T. et al. Scale-specific drivers of kelp forest communities. Oecologia 186(1), 217–233 (2018).ADS 
Article 

Google Scholar 
37.Vergés, A. et al. Long-term empirical evidence of ocean warming leading to tropicalization of fish communities, increased herbivory, and loss of kelp. Proc. Natl. Acad. Sci. 113(48), 13791–13796 (2016).Article 

Google Scholar 
38.Steneck, R. S. et al. Kelp forest ecosystems: Biodiversity, stability, resilience and future. Environ. Conserv. 29(04), 436–459 (2003).Article 

Google Scholar 
39.Nekola, J. C. & White, P. S. The distance decay of similarity in biogeography and ecology. J. Biogeogr. 26(4), 867–878. https://doi.org/10.1046/j.1365-2699.1999.00305.x (1999).Article 

Google Scholar 
40.Claisse, J. T. et al. Biogeographic patterns of communities across diverse marine ecosystems in southern California. Mar. Ecol. 39(S1), e12453. https://doi.org/10.1111/maec.12453 (2018).MathSciNet 
Article 

Google Scholar 
41.Jerde, C. L., Wilson, E. A. & Dressler, T. L. Measuring global fish species richness with eDNA metabarcoding. Mol. Ecol. Resour. 19(1), 19–22. https://doi.org/10.1111/1755-0998.12929 (2019).Article 
PubMed 

Google Scholar 
42.Sigsgaard, E. E. et al. Seawater environmental DNA reflects seasonality of a coastal fish community. Mar. Biol. 164(6), 128. https://doi.org/10.1007/s00227-017-3147-4 (2017).CAS 
Article 

Google Scholar 
43.Nickols, K. J., Wilson White, J., Largier, J. L. & Gaylord, B. Marine population connectivity: Reconciling large-scale dispersal and high self-retention. Am. Nat. 185(2), 196–211. https://doi.org/10.1086/679503 (2015).Article 
PubMed 

Google Scholar 
44.Nickols, K. J., Gaylord, B. & Largier, J. L. The coastal boundary layer: Predictable current structure decreases alongshore transport and alters scales of dispersal. Mar. Ecol. Prog. Ser. 464, 17–35 (2012).ADS 
Article 

Google Scholar 
45.Sassoubre, L. M., Yamahara, K. M., Gardner, L. D., Block, B. A. & Boehm, A. B. Quantification of environmental DNA (eDNA) shedding and decay rates for three marine fish. Environ. Sci. Technol. 50(19), 10456–10464. https://doi.org/10.1021/acs.est.6b03114 (2016).ADS 
CAS 
Article 
PubMed 

Google Scholar 
46.Collins, R. A. et al. Persistence of environmental DNA in marine systems. Commun. Biol. 1(1), 185. https://doi.org/10.1038/s42003-018-0192-6 (2018).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
47.Andruszkiewicz Allan, E., Zhang, W. G., Lavery, C. A. & Govindarajan, F. A. Environmental DNA shedding and decay rates from diverse animal forms and thermal regimes. Environ. DNA 3(2), 492–514. https://doi.org/10.1002/edn3.141 (2021).Article 

Google Scholar 
48.Hansen, B. K., Bekkevold, D., Clausen, L. W. & Nielsen, E. E. The sceptical optimist: Challenges and perspectives for the application of environmental DNA in marine fisheries. Fish Fish. 19(5), 751–768. https://doi.org/10.1111/faf.12286 (2018).Article 

Google Scholar 
49.Weltz, K. et al. Application of environmental DNA to detect an endangered marine skate species in the wild. PLoS One 12(6), e0178124. https://doi.org/10.1371/journal.pone.0178124 (2017).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
50.Fram, J. P. et al. Physical pathways and utilization of nitrate supply to the giant kelp, Macrocystis pyrifera. Limnol. Oceanogr. 53(4), 1589–1603. https://doi.org/10.4319/lo.2008.53.4.1589 (2008).ADS 
CAS 
Article 

Google Scholar 
51.Jackson, G. A. & Winant, C. D. Effect of a kelp forest on coastal currents. Cont. Shelf. Res. 2(1), 75–80 (1983).ADS 
Article 

Google Scholar 
52.Grant, W. D. & Madsen, O. S. The continental-shelf bottom boundary layer. Annu. Rev. Fluid Mech. 18(1), 265–305. https://doi.org/10.1146/annurev.fl.18.010186.001405 (1986).ADS 
MathSciNet 
Article 
MATH 

Google Scholar 
53.Leary, P. R. et al. “Internal tide pools” prolong kelp forest hypoxic events. Limnol. Oceanogr. 62(6), 2864–2878. https://doi.org/10.1002/lno.10716 (2017).ADS 
CAS 
Article 

Google Scholar 
54.Gaylord, B. et al. Spatial patterns of flow and their modification within and around a giant kelp forest. Limnol. Oceanogr. 52(5), 1838–1852 (2007).ADS 
Article 

Google Scholar 
55.Lafferty, K. D., Benesh, K. C., Mahon, A. R., Jerde, C. L. & Lowe, C. G. Detecting Southern California’s white sharks with environmental DNA. Front. Mar. Sci. 5, 355. https://doi.org/10.3389/fmars.2018.00355 (2018).Article 

Google Scholar 
56.Miya, M. et al. MiFish, a set of universal PCR primers for metabarcoding environmental DNA from fishes: Detection of more than 230 subtropical marine species. R. Soc. Open Sci. 2(7), 150088 (2015).ADS 
CAS 
Article 

Google Scholar 
57.Hyde, J. R. & Vetter, R. D. The origin, evolution, and diversification of rockfishes of the genus Sebastes (Cuvier). Mol. Phylogenet. Evol. 44(2), 790–811 (2007).CAS 
Article 

Google Scholar 
58.Min, M. A., Barber, P. H. & Gold, Z. MiSebastes: An eDNA metabarcoding primer set for rockfishes (genus Sebastes). bioRxiv. (2020). http://biorxiv.org/content/early/2020/10/30/2020.10.29.360859.abstract. Accessed January 2021.59.Gold, Z., Sprague, J., Kushner, D. J., Zerecero Marin, E. & Barber, P. H. eDNA metabarcoding as a biomonitoring tool for marine protected areas. PLoS One 16(2), e0238557. https://doi.org/10.1371/journal.pone.0238557 (2021).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
60.Civade, R. et al. Spatial representativeness of environmental DNA metabarcoding signal for fish biodiversity assessment in a natural freshwater system. PLoS One 11(6), e0157366 (2016).Article 

Google Scholar 
61.Berry, T. E. et al. Marine environmental DNA biomonitoring reveals seasonal patterns in biodiversity and identifies ecosystem responses to anomalous climatic events. PLoS Genet. 15(2), e1007943. https://doi.org/10.1371/journal.pgen.1007943 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
62.Ausubel, J. H., Stoeckle, M. Y. & Gaffney, P. Final Report of the 1st US National Conference on Marine Environmental DNA (eDNA). (2019).
63.Reed, D. C. SBC LTER: Reef: Annual time series of biomass for kelp forest species, ongoing since 2000. Environ. Data Initiat. https://doi.org/10.6073/pasta/23965abf42954f345cfd6642fe3c4810 (2018).64.O’Donnell, J. L., Kelly, R. P., Lowell, N. C. & Port, J. A. Indexed PCR primers induce template-specific bias in large-scale DNA sequencing studies. PLoS One 11(3), e0148698 (2016).Article 

Google Scholar 
65.Zhang, J., Kobert, K., Flouri, T. & Stamatakis, A. PEAR: A fast and accurate Illumina Paired-End reAd mergeR. Bioinformatics 30(5), 614–620 (2014).CAS 
Article 

Google Scholar 
66.Rognes, T., Flouri, T., Nichols, B., Quince, C. & Mahé, F. VSEARCH: A versatile open source tool for metagenomics. PeerJ 4, e2584 (2016).Article 

Google Scholar 
67.Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17(1), 10–12 (2011).Article 

Google Scholar 
68.Mahé, F., Rognes, T., Quince, C., de Vargas, C. & Dunthorn, M. Swarm: Robust and fast clustering method for amplicon-based studies. PeerJ 2, e593 (2014).Article 

Google Scholar 
69.Huson, D. H. et al. MEGAN community edition—interactive exploration and analysis of large-scale microbiome sequencing data. PLoS Comput. Biol. 12(6), 1–12 (2016).Article 

Google Scholar 
70.McMurdie, P. J. & Holmes, S. Phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One 8(4), e61217. https://doi.org/10.1371/journal.pone.0061217 (2013).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
71.Crist, T. O., Veech, J. A., Gering, J. C. & Summerville, K. S. Partitioning species diversity across landscapes and regions: A hierarchical analysis of alpha, beta, and gamma diversity. Am. Nat. 162(6), 734–743 (2003).Article 

Google Scholar 
72.Muggeo, V. M. R. Estimating regression models with unknown break-points. Stat. Med. 22(19), 3055–3071 (2003).Article 

Google Scholar 
73.Legendre, P., Borcard, D. & Roberts, D. W. Variation partitioning involving orthogonal spatial eigenfunction submodels. Ecology 93(5), 1234–1240. https://doi.org/10.1890/11-2028.1 (2012).Article 
PubMed 

Google Scholar 
74.Silva, A. R., Dias, C. T. S., Cecon, P. R. & Rêgo, E. R. An alternative procedure for performing a power analysis of Mantel’s test. J. Appl. Stat. 42(9), 1984–1992. https://doi.org/10.1080/02664763.2015.1014894 (2015).MathSciNet 
Article 
MATH 

Google Scholar 
75.Dufrêne, M. & Legendre, P. Species assemblages and indicator species: The need for a flexible asymetrical approach. Ecol. Monogr. 67(3), 345–366 (1997).
Google Scholar 
76.Team, R. C. R: A language and environment for statistical computing. (2018). https://www.r-project.org/. Accessed June 2018.77.Oksanen, J. et al. Package ‘vegan.’ Community Ecol Packag version:2. (2015). More

Study systemWe studied a D. magna population (OHZ) from a small, shallow man-made pond in Oud-Heverlee, Belgium (50°50´ N– 4°39′ E). This pond was constructed for pisciculture in 1970 and has a detailed record of fish-stocking densities for 16 years (Fig. 1a). Dormant stages of D. magna were sampled from three depths of a sediment core, corresponding to three time periods that varied in the level of fish-predation pressure: (1) the pre-fish period (1970–1972), during which no fish were stocked in the pond; (2) the high-fish period (1976–1979), a period with high fish-predation pressure due to intensive fish stocking; (3) the reduced-fish period (1988–1990), with relaxed fish predation pressure due to a reduction in fish stocking (Fig. 1a)9,10,37. This archive was previously sampled using a standard Plexiglas corer with inner diameter of 5.2 cm10. Dating of the sedimentary archive could not be completed with traditional radioisotope analysis, but was based on dry weight and organic matter content under the assumption of a constant sedimentation rate since the establishment of the pond10. The cores contained the full sediment archive, including the transition to the mineral sediment. Sediment cores were aligned using the patterns of Daphnia dormant egg abundance and changes in size of the dormant egg cases as described in Cousyn et al.10. The dormant stages were hatched in the laboratory and taking advantage of the parthenogenetic reproduction mode of D. magna as long as conditions are favorable, we started up clonal lines. The resulting clonal lines are each genetically unique, as dormant stages in D. magna are the result of sexual reproduction. Our approach was thus to sequence the full genome of a random sample of 12 individuals from each of three depth layers of a sediment core representing populations that occurred in three periods with distinct fish-predation pressure.In addition to the reconstruction of temporal genome dynamics, we used twelve regional populations of D. magna distributed along strong environmental gradients of fish-predation pressure in the region. Six populations (DANA, U2, TER1, MO, KNO15, and TER2) were sampled from fishless ponds, while six populations (ZW4, LRV, ZW3, OHN, OM2, and OM3) were sampled from ponds that harbored fish (Supplementary Table 1). These genotypes were hatched from dormant eggs isolated from the upper 2–3 cm of sediment of the study ponds.Whole-genome sequencingTo reconstruct the genomic history, we resequenced the 36 D. magna lines resurrected from the OHZ pond and validated it with additional whole-genome resequencing of 144 D. magna genotypes spread across twelve spatial populations along a fish gradient in the region (Supplementary Table 1). Twelve individuals from each temporal subpopulation of the sediment core and 8–17 individuals per population in the spatial survey were used for genomic DNA extraction using the Nucleo Spin Tissue extraction kit (Macherey-Nagel, Germany), with overnight incubation at 56 °C and following the manufacturer’s instructions. We quantified DNA using PicoGreen reagent (Life Technologies) on a DTX 880 spectrofluorometer (Beckman Coulter). For each sample, 1 µg of gDNA was normalized in a final volume of 50 µl of Tris Buffer, pH 8.5, and sheared using an E220 Focused Ultrasonicator in conjunction with a microTube plate (Covaris) in accordance with the manufacturer’s recommendations. Sheared genomic DNA was assayed on a 2200 TapeStation (Agilent) with High Sensitivity DNA Screentapes to determine the distribution of sheared fragments. The sheared genomic DNA was then prepared into Illumina compatible DNA Sequencing (DNASeq) 100-bp paired-end libraries utilizing NextFlex chemistry (Bio Scientific). Following library construction, libraries were assayed on a 2200 TapeStation (Agilent) with High Sensitivity DNA Screentapes to determine the final library size. Libraries were quantified using the Illumina Library Quantification Kit (Kapa Biosystems) and normalized to an average concentration of 2 nM prior to pooling. Genomic DNA quantification and normalization, shearing setup, library construction, library quantification, library normalization, and library pooling were performed utilizing a Biomek FXP dual-hybrid automated liquid handler (Beckman Coulter). C-Bot (TruSeq PE Cluster Kit v3, Illumina) was used for cluster generation and the Illumina HiSeq2500 platform (TruSeq SBS Kit v3 reagent kit) for paired-end sequencing with 100-bp read length following the manufacturer’s instructions.Short-read mapping and variant callingThe paired end reads (100 × 2) of each individual were first analyzed using FASTQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) for quality checks. Subsequently, low-quality base trimming and adapter cleaning was performed using the Trimmomatic software38. Here, parameter values to remove adapter sequences were chosen for seedMismatches (2), palindromeClipThreshold (30), and simpleClipThreshold (10). The minimum phred quality required to keep a base was set to 28, and the minimum read length to 50 bp. Furthermore, the cleaned reads were mapped to the D. magna genome version 2.439 using Bowtie240 software with very-sensitive parameter settings (-D 20 -R 3 -N 0 -L 20 -i S,1,0.50) and insert size between 200 and 700 bp. The mapped reads were then marked for duplicates using the MarkDuplicates feature of Picard tools (http://broadinstitute.github.io/picard/) to avoid PCR duplicates. The resulting sorted BAM files were then used for variant calling using FreeBayes41. FreeBayes41 is a haplotype-based variant caller that calls SNPs, indel, and complex variants. Minimum base quality was set to 30 with minimum coverage of four reads. We obtained more than 3 × 106 raw variants (3 441 615) for the OHZ temporal subpopulations and 6 × 106 raw variants for the spatial populations.Only biallelic SNPs supported by at least four reads and sequenced in at least 90% of individuals were retained after filtering. The draft genome of Daphnia magna consists of thousands of scaffolds and contigs. To remove repetitive and paralogous regions in the genome, we used the 293 scaffolds greater than 5 kb that altogether represent 84% of the sequenced genome. Further SNP filtering was performed based on D. magna gene models, such that each polymorphic SNP contained within genic regions could be unambiguously assigned to only one gene locus, thereby removing uncertainties attributed to sequence reads mapping to paralogs and to overlapping genes coding on alternative strands of DNA. Finally, SNPs at frequencies below 5% (pooled subpopulations) were removed retaining a total of 724,321 SNPs (mean coverage 20 reads per SNP/individual; 99.6% SNPs with missing values less than 5%) for the temporal analysis of the OHZ population and 748,511 SNPs for the spatial populations. These SNPs were used for downstream analyses.Population differentiationWe calculated both genome-wide and locus-specific levels of genetic differentiation (FST; Weir & Cockerham 198442) using the diffCalc function of the diveRsity43 package in R44. These calculations were performed for each pair of temporal subpopulations (i.e., pairwise FST) in the temporal setting (OHZ) and for six random pairs of nonfish and fish populations in the spatial survey.To calculate allele frequencies in the temporal analysis, we used vcfglxgt function of vcflib (https://github.com/vcflib/vcflib) to set genotypes that are most likely to be true based on maximum genotype likelihood. We then identified the significant differences in allele frequencies between temporal subpopulations over time (P value < 0.01) using a modified chi-square test developed by RS Waples (Waples 1989)11 and implemented in the TAFT software45 (hereafter referred to as Waples test) that accounts for effective population size (Ne), yielding 30,669 significant SNPs (4.23% of total OHZ SNPs) by comparing the prefish and high-fish temporal subpopulation and yielding 11,257 SNPs (1.55% of total OHZ SNPs) for the high-fish and reduced-fish temporal subpopulation comparison; 1771 SNPs showed significant allele-frequency changes in both transitions, most of them also showing a significant reversal in the second transition. To determine whether the observed number of SNPs that showed a significant reversal in allele frequencies is higher than expected by chance, we estimated the null distribution by randomly permuting the temporal subpopulation labels (i.e., prefish, high fish, and reduced fish) and alleles per locus (724,321 SNPs), and recalculating the number of reversals based on Buffalo and Coop 202012.Estimation of effective population size (N e)Effective population sizes (Ne) were calculated from θ = 4Neµ, across the whole genome and with a mutation rate per generation of 4 × 10−946 and a generation time of one year (Daphnia undergoes 10–15 asexual generations and one sexual generation per year), where Ɵ is Watterson’s diversity index and µ is mutation rate. Watterson Ɵ was calculated using the folded SFS option in ANGSD software47 and found to be stable, i.e., near 0.03 across the three subpopulations (prefish, high fish, and reduced fish). The calculated Ne was found to be ~1.66 million in the prefish temporal subpopulation and ~1.72 million for the high-fish and reduced-fish temporal subpopulations. Similarly, for spatial populations, the value ranges from ~1.06 to 1.45 million (Supplementary Table 2).Detecting genomic islands of differentiationFor each scaffold, and for each pairwise comparison among temporal subpopulations and six independent fish and no-fish replicate pairs in the spatial survey, a hidden Markov model (HMM) was used to distinguish genomic regions of high, moderate, or low differentiation among (sub)populations. We used a similar approach as used earlier by Soria-Carrasco et al. (2014)48. In brief, for each of these three levels of genetic differentiation (i.e., the hidden states), a Gaussian distribution of log10(FST + 1) was assumed with the mean and variance initialized as those of the log10(FST + 1) values within each respective level. We then used the Baum–Welch algorithm49 to refine the Gaussian model for each state and the transition matrix among the states. Direct transition from the low to the high state was not allowed. Hidden states were then estimated from the data and we estimated the parameters by the Viterbi algorithm using the R package HiddenMarkov50. A high differentiating island between genomes is defined to contain at least three consecutive SNPs categorized as high-state SNPs by HMM, yielding 6111 and 2879 islands of genomic differentiation between the prefish vs high-fish (mean length: 2428 bp) and the high-fish vs reduced-fish comparison (mean length: 1713 bp), respectively. Similarly, for six independent spatial fish vs no-fish comparisons, the number of islands of differentiation ranged between 4136 and 7493 (with range of mean length 1879–3290 bp), depending upon the comparison (Supplementary Table 3).Functional annotation and enrichment analysisWe investigated the function of the outlier SNPs (P values smaller than 0.01) in the comparison among temporal subpopulations (prefish vs high fish and high fish vs reduced fish) and in HMM-based high-differentiation islands in spatial comparisons. Transcriptome-based functional annotation was performed using the Daphnia magna genome version 2.439. The pathway enrichment analysis was performed using the orthologous genes of D. magna in the D. pulex genome51 based on OrthoDB gene families52 and the KEGG pathway database53. Out of ~29,000 annotated genes of D. magna, 17,400 genes have 17,832 orthologs in D. pulex. However, due to the fragmented status of the D. magna genome assembly, manual curation for high-quality gene models resulted in a total of 12,264 D. magna genes used in our study, of which 2402 genes are annotated to KEGG pathways. Ortholog mapping is not unique. A given gene from the source species, here D. magna, can map to a single, multiple, or no ortholog in the target species, here D. pulex. This can bias statistical tests when referencing to D. pulex genomics resources. We used the number of nonunique mappings for each D. magna gene on the KEGG pathways of D. pulex to weight-adjust the confusion matrix for Fisher’s exact test to obtain the correct P values. Significant pathways are defined as those with FDR corrected (Benjamini–Hochberg method) P values smaller than 0.05. The data analysis was performed using Python packages (NumPy v1.17.454, SciPy v1.4.155, statsmodels v0.11.1, and plotly v4.8.1).Rarefaction analysisRarefaction analyses were used to determine the rate at which outlier SNPs accumulate in the temporal subpopulation or in the set of regional populations as a function of sample size, i.e., the number of individuals sampled from a given population or group of populations. These analyses were performed separately for the prefish as well as for the reduced-fish temporal subpopulation, in both cases to assess the number of individuals needed to accumulate a given percentage of the SNPs that were suggested to be important for the evolution in response to fish through an outlier analysis of the prefish to high-fish transition. With the rarefaction analysis on the prefish population, we estimate the minimum number of individuals that are needed to reach sufficient genetic variation to enable the observed level of adaptation to fish in this population. The rarefaction analysis on the reduced-fish temporal subpopulation was to assess whether the level of genetic polymorphism declined or increased following the period of strong selection by fish. We thus aimed to evaluate how much evolutionary potential a certain number of individuals from the oldest (i.e., before the introduction of fish) as well as the youngest temporal subpopulation (i.e., after a wave of selection) represent. In the first set of analyses, we used the 1109 SNPs belonging to the divergent SNPs that showed significant changes in allele frequencies in the prefish to high-fish transition and also a significant reversal during the high-fish to reduced-fish and were potentially under positive selection (i.e. excluding hitchhiked SNPs) in the OHZ population. This group of SNPs represent polymorphisms that are presumably adaptive or at least contribute to adaptive allelic variants and hence contribute to the adaptive potential of the temporal subpopulations. The analyses were performed by rarefying the genotype matrix of all 12 individuals from either the prefish or the reduced-fish population to all possible (i.e., 4095) subsets of samples of 1–12 genotypes. For each of these subsets we then calculated the average proportion of polymorphic SNPs. These values were plotted against sample size to generate rarefaction curves. Similarly, rarefaction analyses were performed for the 1003 SNPs belonging to the outlier SNPs that were potentially under positive selection (i.e., excluding hitchhiked SNPs) in the OHZ population and that were also present as SNPs in the full spatial dataset (90.4% of the total number of 1109 SNPs). In this case, rarefaction curves were plotted by randomly resampling (1000 times) 1–30 individuals from the total of the 111 individuals of the Leuven regional populations in the spatial data set (i.e., the cluster of populations that represents a sample of nearby populations (within a radius of 10 km) to the focal OHZ population).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More