Ecology

1.Bayne, B. et al. The proposed dropping of the genus Crassostrea for all Pacific cupped oysters and its replacement by a new genus Magallana: a dissenting view. J. Shellfish Res. 36, 545–547 (2017).Article 

Google Scholar 
2.Mann, R. Some biochemical and physiological aspects of growth and gametogenesis in Crassostrea gigas and Ostrea edulis grown at sustained elevated temperatures. J. Mar. Biol. Assoc. UK 59, 95–110 (1979).CAS 
Article 

Google Scholar 
3.Humphreys, J., Herbert, R. J., Roberts, C. & Fletcher, S. A reappraisal of the history and economics of the Pacific oyster in Britain. Aquaculture 428, 117–124 (2014).Article 

Google Scholar 
4.Ellis, T., Gardiner, R., Gubbins, M., Reese, A. & Smith, D. Aquaculture statistics for the UK, with a focus on England and Wales 2012. Centre for Environment Fisheries & Aquaculture Science (Cefas) Weymouth (2015).5.Herbert, R. J. et al. Ecological impacts of non-native Pacific oysters (Crassostrea gigas) and management measures for protected areas in Europe. Biodivers. Conserv. 25, 2835–2865 (2016).Article 

Google Scholar 
6.Reise, K., Buschbaum, C., Büttger, H., Rick, J. & Wegner, K. M. Invasion trajectory of Pacific oysters in the northern Wadden Sea. Mar. Biol. 164, 68 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
7.Geburzi, J. C. & McCarthy, M. L. How do they do it? Understanding the success of marine invasive species. In YOUMARES 8—Oceans Across Boundaries: Learning from each other, 109–124 (Springer, 2018).8.Herbert, R., Roberts, C., Humphreys, J. & Fletcher, S. The Pacific oyster (Crassostrea gigas) in the UK: Economic, legal and environmental issues associated with its cultivation, wild establishment and exploitation. Report for the Shellfish Association of Great Britain (2012).9.Fabioux, C., Huvet, A., Le Souchu, P., Le Pennec, M. & Pouvreau, S. Temperature and photoperiod drive Crassostrea gigas reproductive internal clock. Aquaculture 250, 458–470 (2005).Article 

Google Scholar 
10.Diederich, S., Nehls, G., Van Beusekom, J. E. & Reise, K. Introduced Pacific oysters (Crassostrea gigas) in the northern Wadden Sea: Invasion accelerated by warm summers?. Helgol. Mar. Res. 59, 97 (2005).ADS 
Article 

Google Scholar 
11.Mills, S.R.A. Population structure and ecology of wild Crassostrea gigas (Thunberg, 1793) on the south coast of England. Ph.D. thesis, University of Southampton (2016).12.Dutertre, M., Beninger, P. G., Barillé, L., Papin, M. & Haure, J. Rising water temperatures, reproduction and recruitment of an invasive oyster, Crassostrea gigas, on the French Atlantic coast. Mar. Environ. Res. 69, 1–9 (2010).CAS 
PubMed 
Article 

Google Scholar 
13.Chávez-Villalba, J. et al. Broodstock conditioning of the oyster Crassostrea gigas: Origin and temperature effect. Aquaculture 214, 115–130 (2002).Article 

Google Scholar 
14.Rico-Villa, B., Pouvreau, S. & Robert, R. Influence of food density and temperature on ingestion, growth and settlement of Pacific oyster larvae, Crassostrea gigas. Aquaculture 287, 395–401 (2009).Article 

Google Scholar 
15.Li, G. & Hedgecock, D. Genetic heterogeneity, detected by PCR-SSCP, among samples of larval Pacific oysters (Crassostrea gigas) supports the hypothesis of large variance in reproductive success. Can. J. Fish. Aquat. Sci. 55, 1025–1033 (1998).CAS 
Article 

Google Scholar 
16.Hedge, L. H. & Johnston, E. L. Colonisation of the non-indigenous Pacific oyster Crassostrea gigas determined by predation, size and initial settlement densities. PLoS ONE9 (2014).17.Maurer, D. et al. Reproduction de l’huître creuse dans le Bassin d’Arcachon. Année 2015. Ifremer Report (2016).18.Quayle, D.B. Pacific oyster culture in British Columbia (Department of Fisheries and Oceans, 1988).19.Rico-Villa, B. et al. A flow-through rearing system for ecophysiological studies of Pacific oyster Crassostrea gigas larvae. Aquaculture 282, 54–60 (2008).Article 

Google Scholar 
20.Kheder, R. B., Moal, J. & Robert, R. Impact of temperature on larval development and evolution of physiological indices in Crassostrea gigas. Aquaculture 309, 286–289 (2010).Article 

Google Scholar 
21.Kennedy, V. S. & Breisch, L. L. Maryland’s Oysters: Research and Management Vol. 81 (University of Maryland College Park, Maryland, 1981).
Google Scholar 
22.Helm, M. Cultured aquatic species information programme—Crassostrea gigas. Cultured aquatic species fact sheets. FAO Inland Water Resources and Aquaculture Service (2007).23.Child, A. & Laing, I. Comparative low temperature tolerance of small juvenile European, Ostrea edulis L., and Pacific oysters, Crassostrea gigas Thunberg. Aquacul. Res. 29, 103–113 (1998).Article 

Google Scholar 
24.Strand, A., Waenerlund, A. & Lindegarth, S. High tolerance of the Pacific oyster (Crassostrea gigas, Thunberg) to low temperatures. J. Shellfish Res. 30, 733–735 (2011).Article 

Google Scholar 
25.Rinde, E. et al. Increased spreading potential of the invasive Pacific oyster (Crassostrea gigas) at its northern distribution limit in Europe due to warmer climate. Mar. Freshw. Res. 68, 252–262 (2017).ADS 
Article 

Google Scholar 
26.Wrange, A.-L. et al. Massive settlements of the Pacific oyster, Crassostrea giga, in Scandinavia. Biol. Invasions 12, 1145–1152 (2010).Article 

Google Scholar 
27.Spencer, B., Edwards, D., Kaiser, M. & Richardson, C. Spatfalls of the non-native Pacific oyster, Crassostrea gigas, in British waters. Aquat. Conserv. Mar. Freshw. Ecosyst. 4, 203–217 (1994).Article 

Google Scholar 
28.England, N. Pacific oyster survey of the North East Kent European marine sites. Natural England Commissioned Report NECR016 (2009).29.Smith, I. P., Guy, C. & Donnan, D. Pacific oysters, Crassostrea gigas, established in Scotland. Aquat. Conserv. Mar. Freshw. Ecosyst. 25, 733–742 (2015).Article 

Google Scholar 
30.Cook, E. J. et al. Impacts of climate change on non-native species. Mar. Clim. Change Impact Partnersh. Sci. Rev. 155–166 (2013).31.Cook, E., Beveridge, C., Lamont, P., O’Higgins, T. & Wilding, T. Survey of wild Pacific oyster Crassostrea gigas in Scotland. In Scottish Aquaculture Research Forum Report SARF099 (2014).32.Kochmann, J. Into the wild: documenting and predicting the spread of Pacific oysters (Crassostrea gigas) in Ireland. Ph.D. thesis, University College Dublin (2012).33.Syvret, M., Fitzgerald, A. & Hoare, P. Development of a Pacific oyster aquaculture protocol for the UK: Technical report. Sea Fish Industry Authority, FIFG Project No. 7 (2008).34.d’Auriac, M. B. A. et al. Rapid expansion of the invasive oyster Crassostrea gigas at its northern distribution limit in Europe: Naturally dispersed or introduced? PLoS ONE, 12 (2017).35.Dame, R. F. & Prins, T. C. Bivalve carrying capacity in coastal ecosystems. Aquat. Ecol. 31, 409–421 (1997).Article 

Google Scholar 
36.Leguerrier, D., Niquil, N., Petiau, A. & Bodoy, A. Modeling the impact of oyster culture on a mudflat food web in Marennes-Oléron Bay (France). Mar. Ecol. Prog. Ser. 273, 147–162 (2004).ADS 
Article 

Google Scholar 
37.Forrest, B. M., Keeley, N. B., Hopkins, G. A., Webb, S. C. & Clement, D. M. Bivalve aquaculture in estuaries: Review and synthesis of oyster cultivation effects. Aquaculture 298, 1–15 (2009).Article 

Google Scholar 
38.Ferreira, J. G. et al. Ecological carrying capacity for shellfish aquaculture: Sustainability of naturally occurring filter-feeders and cultivated bivalves. J. Shellfish Res. 37, 709–726 (2018).Article 

Google Scholar 
39.Jordan-Cooley, W. C., Lipcius, R. N., Shaw, L. B., Shen, J. & Shi, J. Bistability in a differential equation model of oyster reef height and sediment accumulation. J. Theor. Biol. 289, 1–11 (2011).MathSciNet 
PubMed 
MATH 
Article 

Google Scholar 
40.Lipcius, R. N. et al. Modeling quantitative value of habitats for marine and estuarine populations. Front. Mar. Sci. 6, 280 (2019).Article 

Google Scholar 
41.Enríquez-Díaz, M., Pouvreau, S., Chávez-Villalba, J. & Le Pennec, M. Gametogenesis, reproductive investment, and spawning behavior of the Pacific giant oyster Crassostrea gigas: Evidence of an environment-dependent strategy. Aquacult. Int. 17, 491–506 (2009).Article 

Google Scholar 
42.Wood, L. E. et al. Unaided dispersal risk of Magallana gigas into and around the UK: Combining particle tracking modelling and environmental suitability scoring. Biological Invasions, 1–20 (2021).43.Hily, C. Prolifération de l’huître creuse du Pacifique Crassotrea gigas sur les côtes manche-atlantique françaises: bilan, dynamique, conséquences écologiques, économiques et ethnologiques, expériences et scénarios de gestion. Rapport LITEAU, 20 (2009).44.McKnight, W. & Chudleigh, I. J. Pacific oyster Crassostrea gigas control within the inter-tidal zone of the North East Kent Marine Protected Areas, UK. Conserv. Evid. 12, 28–32 (2015).
Google Scholar 
45.Brown, J. & Hartwick, E. A habitat suitability index model for suspended tray culture of the Pacific oyster, Crassostrea gigas Thunberg.. Aquacult. Res. 19, 109–126 (1988).Article 

Google Scholar 
46.Diederich, S. High survival and growth rates of introduced Pacific oysters may cause restrictions on habitat use by native mussels in the Wadden Sea. J. Exp. Mar. Biol. Ecol. 328, 211–227 (2006).Article 

Google Scholar 
47.Moran, A. & Manahan, D. Physiological recovery from prolonged ‘starvation’ in larvae of the Pacific oyster Crassostrea gigas. J. Exp. Mar. Biol. Ecol. 306, 17–36 (2004).CAS 
Article 

Google Scholar 
48.Calvo, G. W., Luckenbach, M. W. & Burreson, E. M. A comparative field study of Crassostrea gigas and Crassostrea virginica in relation to salinity in Virginia. Special Report in Applied Marine Science and Ocean Engineering, 349 (1999).49.Petton, B., Boudry, P., Alunno-Bruscia, M. & Pernet, F. Factors influencing disease-induced mortality of Pacific oysters, Crassostrea gigas. Aquacul. Environ. Interact. 6, 205–222 (2015).Article 

Google Scholar 
50.Li, L. et al. Divergence and plasticity shape adaptive potential of the Pacific oyster. Nat. Ecol. Evol. 2, 1751–1760 (2018).PubMed 
Article 

Google Scholar 
51.Ferreira, J., Duarte, P. & Ball, B. Trophic capacity of Carlingford Lough for oyster culture-analysis by ecological modelling. Aquat. Ecol. 31, 361–378 (1997).Article 

Google Scholar 
52.Cognie, B., Haure, J. & Barillé, L. Spatial distribution in a temperate coastal ecosystem of the wild stock of the farmed oyster Crassostrea gigas (Thunberg). Aquaculture 259, 249–259 (2006).Article 

Google Scholar 
53.Enríquez-Díaz, M., Pouvreau, S., Chávez-Villalba, J. & Le Pennec, M. Gametogenesis, reproductive investment, and spawning behavior of the Pacific giant oyster Crassostrea gigas: evidence of an environment-dependent strategy. Aquacult. Int. 17, 491 (2009).Article 

Google Scholar 
54.Ben-Horin, T. et al. Intensive oyster aquaculture can reduce disease impacts on sympatric wild oysters. Aquacul. Environ. Interact. 10, 557–567 (2018).Article 

Google Scholar 
55.Mailleret, L. & Lemesle, V. A note on semi-discrete modelling in the life sciences. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 367, 4779–4799 (2009).ADS 
MathSciNet 
MATH 
Article 

Google Scholar 
56.Powell, E., Klinck, J., Hofmann, E. & Ray, S. Modeling oyster populations. IV: Rates of mortality, population crashes and management. Fish. Bull. 92, 347–373 (1994).
Google Scholar 
57.Wilson, R. A stage-structured oyster population model for reef restoration. Undergraduate Honors Theses Paper, 1403 (2019).58.Guo, X., Hedgecock, D., Hershberger, W. K., Cooper, K. & Jr, S. K. A. Genetic determinants of protandric sex in the Pacific oyster, Crassostrea gigas Thunberg. Evolution 52, 394–402 (1998).59.Morris, D. et al. Cefas coastal temperature network (2016).60.Pouvreau, S. et al. Velyger database: The oyster larvae monitoring French project. SEANOE 10, 41888 (2016).
Google Scholar 
61.Dhoop, T. & Thompson, C. Directional waverider metadata, supplement for QC data download from Realtime Data page. Channel Coastal Observatory (2019).62.Collins, M. et al. Long-term climate change: projections, commitments and irreversibility. In Climate Change 2013-The Physical Science Basis: Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, 1029–1136 (Cambridge University Press, 2013).63.Pastor, D. Reproductive biology of Crassostrea gigas. Ph.D. thesis, University of Southampton (2010).64.Benton, T. G. & Grant, A. Elasticity analysis as an important tool in evolutionary and population ecology. Trends Ecol. Evol. 14, 467–471 (1999).CAS 
PubMed 
Article 

Google Scholar 
65.Grant, A. & Benton, T. G. Elasticity analysis for density-dependent populations in stochastic environments. Ecology 81, 680–693 (2000).Article 

Google Scholar 
66.Caswell, H. & Gassen, N. S. The sensitivity analysis of population projections. Demogr. Res. 33, 801–840 (2015).Article 

Google Scholar 
67.R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria (2019).68.Soetaert, K., Petzoldt, T. & Setzer, R. W. Solving differential equations in R: Package deSolve. J. Stat. Softw. 33, 1–25 (2010).
Google Scholar 
69.Inkscape Project. Inkscape. More

The experiment underlying the study provides a diversity gradient of 1–60 plant species, established in assemblages randomly chosen from a pool of species typical of Arrhenatheretum grasslands. Recently sown on fertile arable soil and maintained by weeding, this experiment is a highly artificial system that fails to meet the definition of semi-natural grasslands7. Four years after establishment, a management intensity gradient of one to four annual cuts and three fertilization levels was established in subplots randomly assigned to the 1–60-species plots. Data presented in this study were collected in the following year.Intensive management was thus imposed on plant species typical of Arrhenaterethum meadows, a plant community characterized by two annual cuts8. The potential effect size of increased management intensity is thus underestimated by applying the management to a plant community not adapted to it. More importantly, it is unlikely that the species-richness of high-diversity plots could be maintained under increased management intensity over longer periods. In fact, 22% of these subplots managed at very high intensity had to be excluded for missing or insufficient yield after only two years, indicating that their species did not persist under high defoliation frequency and fertilizer levels, even when competitors were excluded by weeding.While the discussion hardly addresses this crucial trade-off between management intensity and plant diversity, Schaub et al.6 do indicate that repeated resowing is likely to be necessary to maintain high diversity under increased management intensities. In contrast to permanent grasslands, whose species composition is shaped by site conditions and management, species selection in (re-)sown grasslands is a conscious choice. To be advantageous, mixtures have to show larger yields than the most productive monoculture, so-called transgressive overyielding. Transgressive overyielding is one of the reasons why mixtures, especially grass-clover mixtures, are frequently used in sown grasslands. A European-scale experiment demonstrated that four-species mixtures showed transgressive overyielding at a wide range of sites under intensive agricultural management9,10. Although Schaub et al.6 generally quantify the diversity effects in comparison to monocultures, they argue that grasslands with the high-diversity characteristic of semi-natural grasslands have benefits not only over monocultures but over low-diversity grasslands, such as the 1–8 species standard mixtures shown in Fig. 6 of their paper. However, their results fail to demonstrate that their high-diversity plots show any transgressive overyielding even over monocultures, not to speak of low-diversity mixtures. As species assemblages of the experiment are randomly drawn from the species pool, monocultures and low-diversity mixtures cannot be expected to include the most productive species or species combinations and thus cannot be used to assess transgressive overyielding. When transgressive overyielding was quantified for one- to eight-species plots of the same experiment under extensive management in 2003, it decreased with species number. While two-species mixtures showed a mean transgressive overyielding of 5%, eight-species mixtures were only 70% as productive as the corresponding best monoculture, on average11.Accordingly, the experimental design fails to capture the real trade-offs faced by grassland managers, either in permanent or in sown grassland. It cannot answer if high levels of diversity and the associated biodiversity benefits can be maintained under intensive management for a longer period than just a few years. Neither can it show a productivity benefit of high-diversity grassland assemblages compared to species-poor mixtures, or even monocultures, when in practice the sown species are deliberately chosen rather than randomly drawn from a species pool. While the underlying biodiversity experiment has made valuable contributions to our fundamental understanding of plant diversity effects on ecosystem functioning, it thus cannot be used to derive direct management recommendations for managed grassland. More

Study areaThis study was conducted using tree census data collected from 119 forest inventory plots (73 tropical, 46 subtropical) situated across a latitudinal range of 7.1°N (Colombia) to 27.8°S (Argentina), a longitudinal range of 79.5° to −63.8° W, and an elevation range of 500–3511 m asl (Fig. 1). The mean annual temperature (MAT) of plots ranged from 7.3 to 23.8 °C (mean = 16.7 ± 4.1 °C; mean ± SD) and mean annual precipitation (MAP) of the plots ranged from 608 to 4313 mm y−1 (mean = 1405.0 ± 623.9 mm y−1) (External Databases 1). The number of plots sampled in each country was: Argentina = 46, Bolivia = 26, Peru = 16, Ecuador = 21, and Colombia = 10 (Fig. 1). The 119 forest plots ranged in size from 0.32 to 1.28 ha and represent a cumulative sample area of 104.4 ha (horizontal areas corrected for slope) that containe more than 63,000 trees with a diameter at breast height (DBH, 1.3 m) ≥10 cm (External Database 1). Ninety-four of the plots (79.0%) were ≥1 ha in size. Neither secondary forests nor plantations were included. However, only seven of the plots (five in Argentina and two in Bolivia) were located in forests >100 km2 in extent41, which suggests that at least the edges and borders of some plots could have experienced some degree of disturbance or degradation. All plots were censused at least twice between 1991 and 2017 (census intervals ranged between 2 and 9 years).In each plot, we tagged, mapped, measured, and collected vouchers of all trees and palms (DBH ≥ 10 cm). DBH was measured 50 cm above buttresses or aerial roots when present (where the stem was cylindrical). During the second or subsequent set of censuses, DBH growth, recruitment, and mortality were recorded. In cases where the recorded DBH growth of the second census was less than −0.1 cm y−1 or greater than 7.5 cm y−1, the DBH of the second census was augmented/reduced in order to match these minimum/maximum values42. To homogenize and validate species names of palms and trees recorded in each country and plot, we submitted the combined list from all plots to the Taxonomic Name Resolution Service (TNRS; http://tnrs.iplantcollaborative.org/) version 3.0. Any species with an unassigned TNRS accepted name or with a taxonomic status of ‘no opinion’, ‘illegitimate’, or ‘invalid’ was manually reviewed. Families and genera were changed in accordance with the new species names. If a full species name was not provided or could not be found, the genus and/or family name from the original file was retained.Aboveground carbon stocksThe aboveground biomass (AGB) of each tree was estimated using the allometric equation proposed by Chave et al43., defined as: AGB = 0.0673 × (WD × DBH2 × H)0.976 where AGB (kg) is the estimated aboveground biomass, DBH (cm) is the diameter of the tree at breast height, H (m) is the estimated total height, and WD (g cm−3) is the stem wood density. To estimate WD, we assigned the WD values available in the literature44 to each species found in each plot. In cases where we could not assign a WD value at the species level, we used the average value at the genus- or family level. For unidentified individuals, we used the average WD value of all other species in the plot. Tree height (H) was estimated (see below) based on the heights measured on a subset of the individual stems in each plot using digital hypsometers or clinometers. The estimated AGB of each tree was then converted to units of aboveground carbon (AGC) by applying a conversion factor of 1 kg AGB = 0.456 kg C45. The AGC per ha was then determined by converting kg to Mg, summing the values for all trees in a plot, and extrapolating or interpolating to a sample area of 1 ha.Estimates of AGB and AGC are highly dependent on tree height. Unfortunately, tree height was difficult or impossible to measure on all stems due to physical and logistical constraints. Therefore, we estimated the height of each stem based on allometric relationships between DBH and tree height that we developed for each plot based on height and DBH measurements taken on a subset of individuals. Although the AGB/AGC estimates are only for trees with DBH ≥ 10, we used trees with DBH ≥ 5 cm to construct the H:DBH models when possible in order to be as comparable as possible with the existing pantropical H:DBH models46. In total, 44,442 trees had their heights measured in the field and were employed to construct the H:DBH models. The percentage of trees with direct field measurements of H (DBH ≥ 5 cm) in each country was: Argentina = 19%, Bolivia = 98%, Peru = 96%, Ecuador = 97%, and Colombia = 46%. In Argentina, 32 of 46 plots did not have any field measurements of H, while all plots in all other countries had field measurements of H for at least a subset of trees.We tested and compared the expected effects of using H:DBH models constructed using the local (plot), country, or pantropical (regional) level data. To select the best model to estimate H from DBH at the plot and country level, we used the function modelHD available in the BIOMASS package for R47. We chose the best allometric model from four candidate models (two log-log polynomial models, the three-parameter Weibull model, and a two-parameter Michaelis-Menten model (Supplementary Table 7)) by selecting the model with the lowest RSE and bias (Supplementary Table 8). At the regional level, we used a pantropical model46. The use of country and pantropical H:DBH allometries underestimates tree heights in the lowlands and overestimates tree heights in highlands, thereby homogenizing AGB estimates along elevational gradients10,48 (Supplementary Figs. 11, 12, 13). Using plot level allometries eliminates this problem. However, in the 32 plots in Argentina where we had no information about tree height, we used the country-level H:DBH model developed with the data available in the remaining 14 plots to estimate the height of each tree, which could have homogenized the AGC estimates along the Argentinian elevational gradient (Supplementary Figs. 11, 12, 13).Aboveground carbon dynamicsThe AGC dynamics of each plot was estimated from the annualized values of AGC mortality, AGC productivity (AGC change due to recruitment + growth), and AGC net change3. The calculations of the separate AGC dynamic components was performed as follows: (i) AGC mortality (Mg ha−1 y−1) = the sum of the AGC of all individuals that died between censuses divided by the time between measurements. (ii) AGC recruitment (Mg C ha−1 y−1) = the sum of the AGC of individuals that recruited into DBH ≥ 10 cm between censuses divided by the time between measurements. However, for each tree recruited (DBH ≥ 10 cm), we subtracted the corresponding AGC associated with a tree of 9.99 cm (i.e. just below the detection limit) in order to avoid overestimations of the overall increase in AGC due to recruitment49. (iii) AGC growth (Mg ha−1 y−1) = the sum of the increase in AGC of all individuals with DBH ≥ 10 cm that survived between censuses divided by the time between censuses. (iv) AGC net change (Mg ha−1 y−1) = the difference between AGC stock in the last census (AGCfinal) and AGC stock in the first census (AGC1) divided by the elapsed time (t; in years) between measurements [(AGC net change = AGCfinal − AGC1)/t]. We recognize that these methods exclude C stored in soils or in belowground tissues9,48; however, quantifying just aboveground C stocks and fluxes provides valuable information about the overall status of these forests as net C sinks or sources.ClimateClimate variables at each plot location were extracted from the CHELSA28 bioclimatic rasters at a resolution of 30-arcsec (~1 km2 at the equator). The climate variables extracted were: Mean Annual Temperature (MAT), Mean Diurnal Range (MDR), Isothermality (Isoth), Temperature Seasonality (TS), Maximum Temperature of Warmest Month (MaxTWarmM), Minimum Temperature of Coldest Month (MinTCM), Temperature Annual Range (TAR), Mean Temperature of Wettest Quarter (MeanTWarmQ), Mean Temperature of Driest Quarter (MeanTDQ), Mean Temperature of Warmest Quarter (MeanTWetQ), Mean Temperature of Coldest Quarter (MeanTCQ), Mean Annual Precipitation (MAP), Precipitation of Wettest Month (PWetM), Precipitation of Driest Month (PDM), Precipitation Seasonality (PS), Precipitation of Wettest Quarter (PWetQ), Precipitation of Driest Quarter (PDQ), Precipitation of Warmest Quarter (PWarmQ), Precipitation of Coldest Quarter (PCQ). We separated all variables associated with temperature (°C) from those associated with precipitation (mm y−1) and applied a Principal Component Analysis (PCA) to the 11 variables associated with temperature (PCAtemp) and a separate PCA to the eight variables associated with precipitation (PCAprec). The first two principal components of both PCAtemp and PCAprec (four PCA axes in total) were selected for use in subsequent analyses. Plot elevations were estimated based on their coordinates and the SRTM 1 ArcSec Global V3 (https://lta.cr.usgs.gov) 30 m resolution digital elevation model (DEM).PCAtemp1 (Supplementary Fig. 1a) explained 53.0% of the total variance of the temperature variables and had high loading from Isothermality and Maximum Temperature of Warmest Month, which was primarily associated with changes in elevation (r = −0.97, p  More

Study areaThe study was conducted along the Agi-gawa River, a tributary of the Kiso-gawa River system in central Japan (35°23 42″–35°26 49″N, 137°25 12″–137°28 01″E; Fig. 1), with the Agi-gawa Dam (110 km from the river mouth, 418 m a.s.l.). The Agi-gawa River is a 3rd to 4th-order river with a naturally sand-rich bed derived from weathered granite that characterizes the local geology36. The Agi-gawa Dam (35°25 32″N, 137°25 55″E) had begun operations in 1990; it is a 102 m high rockfill dam with a catchment area of 82 km2, a storage capacity of 4.8 × 107 m3, a mean depth of ~ 45 to 50 m at the dam site, and a hydraulic residence time of 71 days. Although three small sub-dams at the upstream end of the impoundment trap particulates, the sediment speed in the reservoir has been 1,000,000 m3 for 24 years. The dam serves multiple purposes, including flood control, industrial and urban water supply, and the maintenance of baseflow. Further information on the Agi-gawa Dam is available in Katano et al.37.Figure 1The study area shows six study reaches in three stream segments along the Agi-gawa River and Iinuma-gawa Stream, Gifu Prefecture, Japan. Gray circles denote reaches, which are numbered from upstream to downstream within each segment: UD1 and UD2 are upstream of the dam, DD1 and DD2 are downstream of the dam, whilst TR1 and TR2 are in the tributary. The two black circles denote the sediment replenished reaches (S1 and S2). The three small rectangles at the upstream ends of the impoundment are sub-dams, constructed to reduce the inputs of particulates to the impoundment. This map is based on the Digital Topographic Map 25,000 published by Geospatial Information Authority of Japan.Full size imageSediment replenishment and sampling sitesSediment replenishment was undertaken 0.8 and 1.8 km downstream of the Agi-gawa Dam (S1 and S2, Fig. 1) on February 16 and 27, 2005. A total of 1,200 m3 of sediment (D50 ≈ 0.6 mm; mainly sand) was mined from the upstream sub-dams and transported to S1 and S2. We estimate that this constituted 0.086% of the annual sedimentation in the Agi-gawa Dam (e.g., in 2007, replenished sediment per year × 100/sedimentation in the reservoir). The sediment (800 and 400 m3) was replenished at high-flow banks in both sites. The replenished sediment was gradually washed during the high flows at the end of June (visual observation by dam administrators) (Fig. 2). We confirmed that this replenished sediment remained on both banks in March, and no sediments remained on both banks in early July.Figure 2(a) Precipitation (mm·d) (b) mean inflow to the impoundment per day (m3·s-1); and (c) mean outflow from the Agi-gawa Dam per day (m3·s-1). The vertical broken line indicates the study period. Note that the y-axes for (b) and (c) have a logarithmic scale.Full size imageField sampling was conducted twice between March 15 and 18, 2005, prior to sediment flushing and between August 22 and 24, 2005, following sediment flushing [7 weeks after the end of the sediment drift out (Fig. 2)]. The later sampling date was scheduled to investigate the continuous effects (i.e., not immediate effects) of replenished sediment on the riverbed environment and macroinvertebrate assemblages before the replenished sediment had completely been transported further downstream from S1 and S2.Three study segments (length: 1–2 km each) were selected: (1) upstream of the dam and impounded area (UD); (2) downstream of the dam (DD); and (3) in the tributary (TR). These sites were along a 6.0 km stretch of the Agi-gawa River and a 1.0 km stretch of the Iinuma-gawa Stream (catchment area = 24 km2); the latter is a tributary that flows into the Agi-gawa River 2.7 km, downstream of the dam (Fig. 1, Table 1). Each segment contained two study reaches (six reaches in total), and each study reach was 160 m long with several pool–riffle sequences; all reaches were > 300 m apart. DD1 and DD2 were located immediately downstream of the sediment-displaced banks (S1 and S2; 100 m upstream of DD1 and DD2, respectively). Measurements at the two reaches within the same segment were completed on the same day, and the reaches were surveyed in an upstream direction. The dominant land use along the study area was paddy fields, with sparse riparian forest.Table 1 General characteristics of the three study segments and two seasons.Full size tableAlthough the most suitable reference site for DD is the DD prior to dam construction, we were unable to investigate the site prior to dam construction. Therefore, we treated the reference sites as sites that were less affected by the dam than DD on the present day. Katano et al.37 indicated that the difference between the TR and UD sites was smaller than that between DD and UD/TR sites in terms of biota and geology. However, UD was characterized by a wider channel and higher discharge than TR, due to differences in their catchment areas (Table 1). As we did not have a definitive reference, we treated both UD and TR as reference sites (see “Statistical analysis” section). Therefore, how DD in March and DD in August is different from UD and TR can be interpreted as the effect of sediment reduction.Physical environment and water qualitySix riffles were selected at each study reach, and a sampling location (50 × 50 cm quadrat) was established in the mid-channel area of each riffle. Prior to invertebrate sampling, physical environmental variables were measured.Substrate coarseness was measured by gently floating a Plexiglas observation box (50 × 50 × 10 cm deep) divided into four grid squares (25 × 25 cm) on the surface water such that the grid had projected onto the streambed. The size of the substrate material was coded based on the intermediate-axis length: 1 = sand (particles  16 mm) and sieved through a 0.25 mm mesh sieve. Sieved samples and substrate material smaller than pebbles were mixed in a container and preserved in 5% formalin in the field.The material in each container was later divided into two size fractions using 1-and 0.25 mm mesh sieves. To simplify the sorting process, all material retained in the 0.25 mm sieve was mixed and divided into 2n subsamples (maximum n = 32) using a splitter (Idea Co., Tokyo, Japan), following the method described by Vinson and Hawkins43. All macroinvertebrates in subsamples in the 1 mm sieve were counted and identified to the lowest taxonomic level possible, usually to genus or species level using the taxonomic keys of Kawamura and Ueno44, Merritt and Cummins34, Kathman and Brinkhurst45, Kawai and Tanida35, and Torii46.Macroinvertebrate taxa were also classified into five functional feeding groups (FFGs) according to Kawamura and Ueno44, Merritt and Cummins34, Kathman and Brinkhurst45, Kawai and Tanida (2005)35, and Torii46. FFGs were defined as collector-filterers, collector-gatherers, predators, scrapers, and shredders. If a species belonged to ≥ 2 FFGs, the number of individuals was apportioned across the FFGs. We also counted the number of burrowers (#burrowers), inorganic case-bearing caddisflies (#ICB), and net-spinners (#net spinners) of macroinvertebrate assemblages according Kawamura and Ueno44, Merritt and Cummins34, Kathman and Brinkhurst45, Kawai and Tanida35, and Torii46 (see Supplementary Table S1). This classification was carried out as such life-habit traits are important for surviving in a regulated river containing reduced quantities of sand and gravel on the riverbed37. The Chironomidae family was excluded in the life-habit analysis as they consist of various life forms. Once all invertebrates were removed, dry mass (mg m−2) and ash-free dry mass (AFDM, mg m−2) of benthic coarse particulate organic matter (BCPOM,  > 1 mm), and benthic fine particulate organic matter (BFPOM,  0.25 mm) were obtained by drying in an oven at 60 °C for 1 day and combusting in a muffle furnace at 550 °C for 4 h. BCPOM and BFPOM were calculated based on the difference between the dry mass and the AFDM.The total number of invertebrate individuals and the AFDM of BFPOM in each sample were estimated by multiplying by the corresponding 2n value. The number of taxa and density of invertebrates in each sample were calculated as the sum of the values in both size fractions. Additionally, we determined Shannon’s diversity index (H), Simpson’s evenness index, and the percentage of Ephemeroptera, Plecoptera, and Trichoptera (%EPT)47. A sample from UD2 in March had been lost and therefore could not be included in the analyses.Periphyton was sampled from cobbles adjacent to each sampling location. Periphyton was removed from a 5 × 5 cm area on the upper surface of each cobble with a toothbrush. Each sample was placed in a separate container with 250 mL of water. Within 24 h of sample collection, a subsample of the well-mixed content in each container was filtered using a glass-fiber filter (GF/C; Whatman Co., Maidstone, UK). Each filter was placed in a separate vial with 20 mL of 99.5% ethanol and stored in a dark refrigerator at 4 °C for 24 h. The extracted pigments were measured using a spectrophotometer (U-1800; Shimadzu Co., Kyoto, Japan), following the method of Lorenzen48.Analysis of case materials of an inorganic case-bearing caddisflyWe compared the particle size structure of replenished sediment, riverbed sediment, and case materials for case-bearing caddisfly. The replenished sediment was directly sampled in a 1 L polyethylene jar at the upstream replenished bank (S1) on March 16, 2005 (Fig. 1). Riverbed sediment was sampled at two stations; 100 m upstream of S1, and 100 m upstream of DD1 between August 22 and 24, 2005. At each station of the river, a metallic narrow cup (200 mL) with a lid was pushed into a vacancy between the cobbles, which had been randomly selected, and fine sediments (up to small gravel) in the vacancies were sampled by closing the lid underwater. Sampling was carried out three times (i.e., three different vacancies in the cobbles), and subsamples were pooled for measurement. The replenished and riverbed sediment was combusted at 550 °C for 2 h in a muffle furnace to remove organic contamination. Combusted samples were separated with eight sieves with a mesh size range of 0.075–9.5 mm (JIS A 1204). Each fraction was weighed, and the grain size accumulation curve of each type of sediment and its D50 were obtained.In a macroinvertebrate sample at DD1 between August 22 and 24, 2005, ten individuals from two case-bearing caddisfly larvae, Glossosoma sp. and Gumaga orientalis, which were prevalent at DD1 during this period (see Results), were randomly selected from the formalin-fixed sample. The case was carefully removed from the larvae and combusted as described above for the replenished and riverbed sediment. The number of case material grains was measured using a dissection microscope.Statistical analysesWe described results based on two main assumptions: (1) the DD in March is the dam-affected reach (cf. unregulated reaches UD and TR), and (2) the changes in DD from March to August were mainly a result of sediment replenishment. In the statistical analyses, the p criterion (⍺) was set at 0.05.To consider the effects of the segment, replicate reach, and season on variables, nested multivariate analysis of variance (MANOVA) was used to test whether any measured variables at the riffle scale differed between segments (UD, DD, and TR). Three segments and two replicate reaches were nested within each season (March and August) and segment (i.e., Season/Segment/Reach), whereby measurements within each reach were treated as subsamples. In the MANOVA, we also consider the interactions of the variances to interpret the interactions among the sampling segments and seasons to consider the independent effects on the factors.To perform MANOVA, we assumed that temporal variability was greater than spatial variability within each reach for variables measured over 24 h (e.g., water quality), and the opposite would hold true for variables measured only once (e.g., macroinvertebrates). Therefore, subsamples within each reach were either spatially or temporally replicated, depending on the variable type. Temporal replicates (four samples collected every 6 h) were treated as a repeated factor (time factor). A nested MANOVA was used for variables quantified once at each location (e.g., macroinvertebrates), and nested repeated-measures MANOVA (rm-MANOVA) were used for variables quantified over a 24 h period at each reach (e.g., water quality). When a significant difference was detected by MANOVA with non-significant interactions, each variable was tested separately with a nested ANOVA for variable groups once at each location or the nested rm-ANOVA for repeated-measured variables, as appropriate for the particular variable. The risk of inflating Type 1 errors for the ANOVA was reduced using Bonferroni adjustments.These MANOVA and ANOVA tests were conducted with R version 3.6.049. The residuals of each variable in each MANOVA and ANOVA model were verified using the Shapiro–Wilk normality test prior to analyses, and normality was improved using arcsine(x) or log (x + 1) transformation when appropriate.Tukey’s multiple comparison test in a one-way ANOVA model (Season/Segment/Reach) was used for comparisons between segments. Any significant changes in values for variables from UD to DD were interpreted as the effects of the dam based on the assumption that conditions in UD and DD were similar prior to dam construction; this was because replenished sediment had not been supplied in March (see before). However, UD may be unsuitable as a reference site compared with TR as the former may be at least partly affected by the dam. This may particularly be the case for benthic invertebrates, such as the interruption of the upstream flight of adult females50. Therefore, UD and TR were treated as reference sites for reservoir and tributary effects, respectively. This was because both were unaffected by the dam, and sediment replenishment as tributaries may function as sites for resource recovery for the dam-affected mainstem of the river37,51,52, despite differing watershed areas. Therefore, the similarity of variables between the TR and UD sites was statistically confirmed such that they could be treated as reference sites. As such, the recovery from March to August could reliably demonstrate the effect of sediment replenishment. For example, although the value at DD differed from that at TR and/or UD in March, it was similar to that at UD and/or TR in August.Multivariate analyses were conducted using the R “vegan” package version 2.5.6 to compare invertebrate assemblage structures between segments. Bray–Curtis coefficients based on species abundance were used to calculate a dissimilarity matrix, and dissimilarities between UD and DD, and between TR and DD in each season were tested using two-way ANOVA and Tukey post-hoc tests.Macroinvertebrate assemblage organization in relation to environmental gradients was analyzed using redundancy analysis (RDA) with the “rda” function of “vegan” package. This was because the preliminary analysis using detrended correspondence analysis (DCA) showed that the gradient lengths of DCA were More

Experimental study and protein domain analysisInsectsHosts and parasitoids were maintained as previously described25. Five Leptopilina heterotoma (Hymenoptera: Figitidae) populations were used for experiments: a population from Japan (Sapporo), two populations from the United Kingdom (1: Whittlesford; 2: Great Shelford) and two populations from Belgium (1: Wilsele; 2: Eupen). Information on collection sites, including GPS coordinates, can be found in25.Determination of host fat contentD. simulans and D. melanogaster hosts were allowed to lay eggs during 24 h in glass flasks containing ~ 50 mL standard medium25. After two days, developing larvae were sieved and ~ 200 were larvae placed in a Drosophila tube containing ~ 10 mL medium. Seven days after egg laying, newly formed pupae were frozen at – 18 °C, after which fat content was determined as described in25, where dry weight before and after neutral fat extraction was used to calculate absolute fat amount (in μg) for each host. The host pupal stage was chosen for estimating fat content, because at this point the host ceases to feed, while the parasitoid starts consuming the entire host36. All data were analysed using R Project version 3.4.360. Fat content of hosts was compared using a one-way ANOVA with host species as fixed factor.Manipulation of host fat contentTo generate leaner D. melanogaster hosts, we adapted our standard food medium25 to contain 100 times less (0.5 g) sugar per litre water. Manipulating sugar content did not alter the structure of the food medium, thus maintaining similar rearing conditions, with the exception of sugar content. Fat content of leaner and fatter D. melanogaster hosts was determined and analysed as described above.Fat synthesis quantification of wasp populationsMated female L. heterotoma were allowed to lay eggs on host fly larvae collected as described above with ad libitum access to honey as a food source until death. Honey consists of sugars and other carbohydrates that readily induce fat synthesis. After three weeks, adult offspring emergence was monitored daily and females were haphazardly placed in experimental treatments. Females were either killed at emergence (to measure teneral lipid reserves) or after feeding for 7 days on honey. Wasps were frozen at − 18 °C after completion of experiments. Fat content was determined as described above for hosts. The ability for fat synthesis was then determined by comparing mean fat levels of recently emerged compared to fed individuals, similar to procedures described in15,25,28. An increase in fat levels after feeding is indicative of active fat synthesis; equal or lower fat levels suggest fat synthesis did not take place. Each population tested on D. melanogaster or D. simulans represented an independent dataset that was analysed separately, as in Visser et al. 201825, because we are interested in the response of each population on each host species. We used T-tests when data was normally distributed and variances equal, log-transformed data for non-normal data, and a Welch’s t-test when variances were unequal. We corrected for multiple testing using Benjamini and Hochberg’s False Discovery Rate61.Fat synthesis quantification using a familial design and GC–MS analysesTo tease apart the effect of wasp genotype and host environment, we used a split-brood design where the offspring of each mother developed on lean D. simulans or fat D. melanogaster hosts in two replicated experiments (experiment 1 and 2). In both experiments, mothers were allowed to lay eggs in ~ 200 2nd to 3rd instar host larvae of one species for four days, after which ~ 200 host larvae of the other species were offered during four days. The order in which host larvae were presented was randomized across families. Following offspring emergence, daughters were allocated into two treatment groups: a control where females were fed a mixture of honey and water (1:2 w/w) or a treatment group fed a mixture of honey and deuterated water (Sigma Aldrich) (1:2 w/w; stable isotope treatment) for 7 days. Samples were prepared for GC–MS as described in 28. Incorporation of up to three deuterium atoms can be detected, but percent incorporation is highest when only 1 deuterium atom is incorporated. As incorporation of a single atom unequivocally demonstrates active fat synthesis, we only analysed percent incorporation (in relation to the parent ion) for the abundance of the m + 1 ion. Percent incorporation was determined for five fatty acids, C16:1 (palmitoleic acid), C16:0 (palmitate), C18:2 (linoleic acid), C18:1 (oleic acid), and C18:0 (stearic acid), and the internal standard C17:0 (margaric acid). Average percent incorporation for C17:0 was 19.4 (i.e. baseline incorporation of naturally occurring deuterium) and all values of the internal standard remained within 3 standard deviations of the mean (i.e. 1.6). Percent incorporation of control samples was subtracted from treatment sample values to correct for background levels of deuterium (i.e. only when more deuterium is incorporated in treatment compared to controls fatty acids are actively being synthesized). For statistical analyses, percent incorporation was first summed for C16:1, C16:0, C18:2, C18:1 and C18:0 to obtain overall incorporation levels, as saturated C16 and C18 fatty acids are direct products of the fatty acid synthesis pathway (that can subsequently be desaturated).Data (presented in Fig. 1) was analysed by means of a linear mixed effects model (GLMM, lme4 package) with host (lean D. simulans and fat D. melanogaster) and experiment (conducted twice) as fixed effect, family nested within population (Japan, United Kingdom 1 and 2, Belgium 1 and 2) as random factor, and percentage of incorporation of stable isotopes as dependent variable (log transformed; n = 138). Non-significant terms (i.e., experiment) were sequentially removed from the model to obtain the minimal adequate model as reported in Table 2. When referring to “families,” we are referring to the comparison of daughters of singly inseminated females, which (in these haplodiploid insects) share 75% of their genome.Identification of functional acc and fas genes in distinct parasitoid speciesTo obtain acc and fas nucleotide sequences for L. clavipes, G. legneri, P. maculata and A. bilineata, we used D. melanogaster mRNA ACC transcript variant A (NM_136498.3 in Genbank) and FASN1-RA (FBtr0077659 in FlyBase) and blasted both sequences against transcripts of each parasitoid (using the blast function available at http://www.parasitoids.labs.vu.nl62,63). Each nucleotide sequence was then entered in the NCBI Conserved Domain database64 to determine the presence of all functional protein domains. All sequences were then translated using the Expasy translate tool (https://web.expasy.org/translate/), where the largest open reading frame was selected for further use and confirming no stop codons were present. Protein sequences were then aligned using MAFFT v. 7 to compare functional amino acid sequences between all species (Supplementary files 1 and 2)65.Simulation studyWe consider the general situation where phenotypic plasticity is only sporadically adaptive and ask the question whether and under what circumstances plasticity can remain functional over long evolutionary time periods when the regulatory processes underlying plasticity are gradually broken down by mutations. We consider a regulatory mechanism that switches on or off a pathway (like fat synthesis) in response to environmental conditions (e.g., host fat content).Fitness considerationsWe assume that the local environment of an individual is characterized by two factors: fat content F and nutrient content N, where nutrients represent sugars and other carbohydrates that can be used to synthesize fat. Nutrients are measured in units corresponding to the amount of fat that can be synthesized from them. We assume that fitness (viability and/or fecundity) is directly proportional to the amount of fat stored by the individual. When fat synthesis is switched off, this amount is equal to F, the amount of fat in the environment. When fat synthesis is switched on, the amount of fat stored is assumed to be (N – c + (1 – k)F). This expression reflects the following assumptions: (i) fat is synthesized from the available nutrients, but this comes at a fitness cost c; (ii) fat can still be absorbed from the environment, but at a reduced rate ((1 – k)). It is adaptive to switch on fat synthesis if (N – c + (1 – k)F) is larger than F, or equivalently if (F < tfrac{1}{k}(N - c)).The right-hand side of this inequality is a straight line, which is illustrated by the blue line in Fig. 4. The three boxes in Fig. 4 illustrate three types of environmental conditions. Red box low-fat environments. Here, (F < tfrac{1}{k}(N - c)) is always satisfied, implying that fat synthesis should be switched on constitutively. Yellow box high-fat environments. Here, (F > tfrac{1}{k}(N – c)), implying that fat synthesis should be switched off constitutively.

Orange box intermediate-fat environments. Here, fat synthesis should be plastic and switched on if for the given environment (N, F) the fat content is below the blue line and switched off otherwise.

Figure 4Environmental conditions encountered by the model organisms. For a given combination of environmental nutrient content N and environmental fat content F, it is adaptive to switch on fat synthesis if (N, F) is below the blue line (corresponding to (F < tfrac{1}{k}(N - c))) and to switch it off otherwise. The three boxes illustrate three types of environment: a low-fat environment (red) where fat synthesis should be switched on constitutively; a high-fat environment (yellow) where fat synthesis should be switched off constitutively; and an intermediate-fat environment (orange) where a plastic switch is selectively favoured.Full size imageThe simulations reported here were all run for the parameters (k = tfrac{1}{2}{text{ and }}c = tfrac{1}{4}). We also investigated many other combinations of these parameters; in all cases, the results were very similar to those reported in Fig. 3.Gene regulatory networks (GRN)In our model, the switching device was implemented by an evolving gene regulatory network (as in van Gestel and Weissing66). The simulations shown in Fig. 3 of the main text are based on the simplest possible network that consists of two receptor nodes (sensing the fat and the nutrient content in the local environment, respectively) and an effector node that switches on fat synthesis if the combined weighted input of the two receptor nodes exceeds a threshold value T and switches it off otherwise. Hence, fat synthesis is switched on if (w_{F} F + w_{N} N > T) (and off otherwise). The GRN is characterized by the weighing factors (w_{F} {text{ and }}w_{N}) and the threshold T. These parameters are transmitted from parents to offspring, and they evolve subject to mutation and selection. We also considered alternative network structures (all with two receptor nodes and one effector node, but with a larger number of evolvable weighing factors67, and obtained very similar results, see below).For the simple GRN described above, the switching device is 100% adaptive when the switch is on (i.e., (w_{F} F + w_{N} N > T)) if (F < tfrac{1}{k}(N - c)) and off otherwise. A simple calculation yields that this is the case if: (w_{N} > 0{, }w_{F} = – k{kern 1pt} w_{N} {text{ and }}T = c{kern 1pt} w_{N}).Evolution of the GRNFor simplicity, we consider an asexual haploid population with discrete, non-overlapping generations and fixed population size (N = 10,000). Each individual has several gene loci, each locus encoding one parameter of the GRN. In case of the simple network described above, there are three gene loci, each with infinitely many alleles. Each individual harbours three alleles, which correspond to the GRN parameters (w_{F} {, }w_{N} {text{ and }}T), and hence determine the functioning of the genetic switch. In the simulations, each individual encounters a randomly chosen environment ((N{, }F)). Based on its (genetically encoded) GRN, the individual decides on whether to switch on or off fat synthesis. If synthesis is switched on, the individual’s fitness is given by (N – c + (1 – k)F); otherwise its fitness is given by F. Subsequently, the individuals produce offspring, where the number of offspring produced is proportional to the amount of fat stored by an individual. Each offspring inherits the genetic parameters of its parent, subject to mutation. With probability μ (per locus) a mutation occurs. In such a case the parental value (in case of a simple network: the parent’s allelic value (w_{F} {, }w_{N} {text{ or }}T)) is changed to a mutated value ((w_{F} { + }delta {, }w_{N} { + }delta {text{ or }}T + delta)), where the mutational step size δ is drawn from a normal distribution with mean zero and standard deviation σ. In the reported simulations, we chose (mu = 0.001) and (sigma = 0.1). The speed of evolution is proportional to (mu cdot sigma^{2}), implying that the rate of change in Fig. 3 (both the decay of plasticity and the rate of regaining adaptive plasticity) are positively related to μ and σ.Preadaptation of the GRNsStarting with a population with randomly initialized alleles for the GRN parameters, we first let the population evolve for 10,000 generations in the intermediate-fat environment (the orange box in Fig. 4). In all replicate simulations, a “perfectly adapted switch” (corresponding to (w_{N} > 0{, }w_{F} = – k{kern 1pt} w_{N} {text{ and }}T = c{kern 1pt} w_{N})) evolved, typically within 1,000 generations. Still, the evolved GRNs differed across replicates, as they evolved different values of (w_{N} > 0). These evolved networks were used to seed the populations in the subsequent “decay” simulations.Evolutionary decay of the GRNsFor the decay experiments reported in Fig. 3 of the main text, we initiated a large number of monomorphic replicate populations with one of the perfectly adapted GRNs from the preadaptation phase. These populations were exposed for an extended period of time (1,000,000 generations) to a high-fat environment (the yellow box in Fig. 4), where all preadapted GRNs switched off fat synthesis. However, in some scenarios, the environmental conditions changed back sporadically (with probability q) to the intermediate-fat environment (the orange box in Fig. 4), where it is adaptive to switch on fat metabolism in 50% of the environmental conditions (when (N, F) is below the blue line in Fig. 4). In Fig. 3, we report on the changing rates (q = 0.0) (no changing back; red), (q = 0.001) (changing back once every 1,000 generations; purple), and (q = 0.01) (changing back once every 100 generations; pink). When such a change occurred, the population was exposed to the intermediate-fat environment for t generations (Fig. 3 is based on t = 3).Throughout the simulation, the performance of the network was monitored every 100 generations as follows: 100 GRNs were chosen at random from the population, and each of these GRNs was exposed to 100 randomly chosen environmental conditions from the intermediate-fat environment (orange box in Fig. 4). From this, we could determine the average percentage of “correct” decisions (where the network should be switched on if and only if (F < tfrac{1}{k}(N - c)). 1.0 means that the GRN is still making 100% adaptive decisions; 0.5 means that the GRN only makes 50% adaptive decision, as would be expected by a random GRN or a GRN that switches the pathway constitutively on or off. This measure for performance in the “old” intermediate-fat environment was determined for 100 replicate simulations per scenario and plotted in Fig. 3 (mean ± standard deviation).Evolving robustness of the GRNsThe simulations in Fig. 3 are representative for all networks and parameters considered. Whenever (q = 0.0), the performance of the regulatory switch eroded in evolutionary time, but typically at a much lower rate in case of the more complex GRNs. Whenever (q = 0.01), the performance of the switch went back to levels above 90% and even above 95% for the more complex GRNs. Even for (q = 0.001), a sustained performance level above 75% was obtained in all cases.Intriguingly, in the last two scenarios the performance level first drops rapidly (from 1.0 to a much lower level, although this drop is less pronounced in the more complex GRNs) and subsequently recovers to reach high levels again. Apparently, the GRNs have evolved a higher level of robustness, a property that seems to be typical for evolving networks8. For the simple GRN studied in Fig. 3, this outcome can be explained as follows. The initial network was characterized by the genetic parameters (w_{N} > 0{, }w_{F} = – k{kern 1pt} w_{N} {text{ and }}T = c{kern 1pt} w_{N}) (see above), where (w_{N}) was typically a small positive number. In the course of evolutionary time, the relation between the three evolving parameters remained approximately the same, but (w_{N}) (and with it the other parameters) evolved to much larger values. This automatically resulted in an increasingly robust network, since mutations with a given step size distribution affect the performance of a network much less when the corresponding parameter is large in absolute value.Costs of plasticityPhenotypically plastic organisms can incur different types of costs68. In our simple model, we only consider the cost of phenotype-environment mismatching, that is, the costs of expressing the ‘wrong’ phenotype in a given environment. When placed in a high-fat environment, the preadapted GRNs in our simulations take the ‘right’ decision to switch off fat metabolism. Accordingly, they do not face any costs of mismatching. Yet, the genetic switch rapidly decays (as indicated in Fig. 3 by the rapid drop in performance when tested in an intermediate-fat environment), due to the accumulation of mutations.It is not unlikely that there are additional fitness costs of plasticity, such as the costs for the production and maintenance of the machinery underlying plasticity68. In the presence of such constitutive costs, plasticity will be selected against when organisms are living in an environment where only one phenotype is optimal (as in the high- and low-fat environments in Fig. 4). This would obviously affect the evolutionary dynamics in Fig. 3, but the size of the effect is difficult to judge, as the constitutive costs of plasticity are notoriously difficult to quantify. In case of the simple switching device considered in our model, we consider the constitutive costs of plasticity as marginal, but these costs might be substantial in other scenarios. More

A reef manta ray visits a cleaning station at Lady Elliot Island, Australia. Credit: A. O. Armstrong et al./Ecol. Evol. (CC BY 4.0)

Ecology
08 April 2021
What manta rays remember: the best spots to get spruced up

Giant fish preserve a mental map of where cleaning fish provide the highest-quality pest removal.

Share on Twitter
Share on Twitter

Share on Facebook
Share on Facebook

Share via E-Mail
Share via E-Mail

Even sea creatures need pampering. Manta rays make regular visits to ‘cleaning stations’, where small fish rid the rays of skin parasites at the coral-reef equivalent of a day spa. Now it seems that rays can identify and remember spots where they have received quality cleaning.Cleaning stations are often centred on corals inhabited by cleaner shrimp or fish. To understand how these stations influence rays’ movements, Asia Armstrong at the University of Queensland in St Lucia, Australia, and her colleagues tracked 34 reef manta rays (Mobula alfredi) off the coast of eastern Australia for about 1.5 years.The highest density of rays was found at places where cleaning fish called blue-streak cleaner wrasses (Labroides dimidiatus) were most abundant. Rays typically visited cleaning stations during the day, when cleaner wrasses are most active, and favoured stations close to foraging regions.Rays are thought to prefer stations that provide superior cleaning — where the cleaners don’t bite them, for example. The rays’ behaviour suggests that they have a mental map of spots that offer both high-quality cleaning and proximity to foraging grounds.

Ecol. Evol. (2021)

Ecology More

1.Fushimi, H. Production of juvenile marine finfish for stock enhancement in Japan. Aquaculture 200, 33–53. https://doi.org/10.1016/S0044-8486(01)00693-7 (2001).Article 

Google Scholar 
2.Masuda, R. & Tsukamoto, K. Stock enhancement in Japan: review and perspective. Bull. Mar. Sci. 62, 337–358. https://www.ingentaconnect.com/content/umrsmas/bullmar/1998/00000062/00000002/art00005 (1998).3.Sekino, M., Saitoh, K., Yamada, T., Hara, M. & Yamashita, Y. Genetic tagging of released Japanese flounder (Paralichthys olivaceus) based on polymorphic DNA markers. Aquaculture 244, 49–61. https://doi.org/10.1016/j.aquaculture.2004.11.006 (2005).CAS 
Article 

Google Scholar 
4.Arnold, W. S. Bivalve enhancement and restoration strategies in Florida, USA. Hydrobiologia 465, 7–19. https://doi.org/10.1023/A:1014596909319 (2001).Article 

Google Scholar 
5.Castell, L. L., Naviti, W. & Nguyen, F. Detectability of cryptic juvenile Trochus niloticus Linnaeus in stock enhancement experiments. Aquaculture 144, 91–101. https://doi.org/10.1016/S0044-8486(96)01320-8 (1996).Article 

Google Scholar 
6.McCormick, T. B., Herbinson, K., Mill, T. S. & Altick, J. A review of abalone seeding, possible significance and a new seeding device. Bull. Mar. Sci. 55, 680–693. https://www.ingentaconnect.com/content/umrsmas/bullmar/1994/00000055/f0020002/art00035 (1994).7.Zohar, Y. et al. The Chesapeake Bay blue crab (Callinectes sapidus): A multidisciplinary approach to responsible stock replenishment. Rev. Fish. Sci. 16, 24–34. https://doi.org/10.1080/10641260701681623 (2008).Article 

Google Scholar 
8.Funge-Smith, S., Briggs, M. & Miao, W. Regional Overview of Fisheries and Aquaculture in Asia and the Pacific 2012 (RAP Publication (FAO), 2012).
Google Scholar 
9.Mao, Y. et al. Chapter 4. Bivalve production in China. In Goods and Services of Marine Bivalves (eds Smaal, A. C. et al.) 51–72 (Springer, 2019).
Google Scholar 
10.Ran, Z. et al. Fatty acid and sterol changes in razor clam Sinonovacula constricta (Lamarck 1818) reared at different salinities. Aquaculture 473, 493–500. https://doi.org/10.1016/j.aquaculture.2017.03.017 (2017).CAS 
Article 

Google Scholar 
11.Suzuki, T., Inoue, K. & Ozawa, T. Environmental degradation in Ise and Mikawa Bays after 1960’s as viewed from intertidal molluscan community. Boll. Nagoya Univ. Museum 22, 31–64. https://doi.org/10.18999/bulnum.022.04 (2006).Article 

Google Scholar 
12.Nakamura, T. et al. Marine reservoir effect deduced from 14C dates on marine shells and terrestrial remains at archeological sites in Japan. Nucl. Instrum. Methods Phys. Res. B 259, 453–459. https://doi.org/10.1016/j.nimb.2007.01.186 (2007).ADS 
CAS 
Article 

Google Scholar 
13.Ministry of the Environment, Japan. Ariake Sea and Yatsushiro Sea Comprehensive Survey Evaluation Committee Report. https://www.env.go.jp/council/20ari-yatsu/rep061221/all.pdf (2006).14.Ito, S., Eguchi, T. & Kawahara, I. Rearing experiment on planktonic larvae of the JackKnife clam, Sinonovacula constricta. Boll. Saga Prefect. Ariake Fish. Res. Dev. Cent. 20, 49–53. https://agriknowledge.affrc.go.jp/RN/2010813185.pdf (2001).15.Ohkuma, H., Ymagachi, T., Kawahara, I., Eguchi, T. & Ito, S. A study on the development of techniques for mass production of seeds of jackknife clam, Sinonovacula constricta. Boll. Saga Prefect. Ariake Fish. Res. Dev. Cent. 22, 47–54. https://agriknowledge.affrc.go.jp/RN/2030813210.pdf (2004).16.Tsukuda, M. et al. Variation in the Distribution of the JackKnife Clam, Sinonovacula constricta, on the Muddy Tidal flat of Ariake Sound off Saga Prefecture. Boll. Saga Prefect. Ariake Fish. Res. Dev. Cent. 28, 47–49. https://agriknowledge.affrc.go.jp/RN/2010925685.pdf (2017).17.Holman, L. E., de la Garcia, S. D., Onoufriou, A., Hillestad, B. & Johnston, I. A. A workflow used to design low density SNP panels for parentage assignment and traceability in aquaculture species and its validation in Atlantic salmon. Aquaculture 476, 59–64. https://doi.org/10.1016/j.aquaculture.2017.04.001 (2017).CAS 
Article 

Google Scholar 
18.Li, Y. H. & Wang, H. P. Advances of genotyping-by-sequencing in fisheries and aquaculture. Rev. Fish Biol. Fish. 27, 535–559. https://doi.org/10.1007/s11160-017-9473-2 (2017).Article 

Google Scholar 
19.Robledo, D., Palaiokostas, C., Bargelloni, L., Martínez, P. & Houston, R. Applications of genotyping by sequencing in aquaculture breeding and genetics. Rev. Aquac 10, 670–682. https://doi.org/10.1111/raq.12193 (2018).Article 
PubMed 

Google Scholar 
20.You, X., Shan, X. & Shi, Q. Research advances in the genomics and applications for molecular breeding of aquaculture animals. Aquaculture 526, 735357. https://doi.org/10.1016/j.aquaculture.2020.735357 (2020).CAS 
Article 

Google Scholar 
21.Arthington, A. H. Ecological and genetic impacts of introduced and translocated freshwater fishes in Australia. Can. J. Fish. Aquat. Sci. 48, 33–43. https://doi.org/10.1139/f91-302 (1991).Article 

Google Scholar 
22.Habtemariam, B. T., Arias, A., García-Vázquez, E. & Borrell, Y. J. Impacts of supplementation aquaculture on the genetic diversity of wild Ruditapes decussatus from northern Spain. Aquacult. Environ. Interact. 6, 241–254. https://doi.org/10.3354/aei00128 (2015).Article 

Google Scholar 
23.Williams, S. L. & Orth, R. J. Genetic diversity and structure of natural and transplanted eelgrass populations in the Chesapeake and Chincoteague Bays. Estuaries 21, 118–128. https://doi.org/10.2307/1352551 (1998).Article 

Google Scholar 
24.Yamakawa, A. Y. & Imai, H. PCR-RFLP typing reveals a new invasion of Taiwanese Meretrix (Bivalvia: Veneridae) to Japan. Aquat. Invasions 8, 407–415. https://doi.org/10.3391/ai.2013.8.4.04 (2013).Article 

Google Scholar 
25.Niu, D. H., Feng, B. B., Liu, D. B., Zhong, Y. M., Shen, H. D. & Li, J. L. Significant genetic differentiation among ten populations of the razor clam Sinonovacula constricta along the coast of china revealed by a microsatellite analysis. Zool. Stud. 51, 406–414. http://zoolstud.sinica.edu.tw/Journals/51.3/406.pdf. (2012).26.DeFaveri, J., Shikano, T., Ghani, N. I. A. & Merilä, J. Contrasting population structures in two sympatric fishes in the Baltic Sea basin. Mar. Biol. 159, 1659–1672. https://doi.org/10.1007/s00227-012-1951-4 (2012).Article 

Google Scholar 
27.Jeffery, N. W. et al. RAD sequencing reveals genomewide divergence between independent invasions of the European green crab (Carcinus maenas) in the Northwest Atlantic. Ecol. Evol. 7, 2513–2524. https://doi.org/10.1002/ece3.2872 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
28.Kato, D. et al. Evaluation of the population structure and phylogeography of the Japanese Genji firefly, Luciola cruciata, at the nuclear DNA level using RAD-Seq analysis. Sci. Rep. 10, 1533. https://doi.org/10.1038/s41598-020-58324-9 (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
29.Yokota, M., Harada, Y. & Iizuka, M. Genetic drift in a hatchery and the maintenance of genetic diversity in hatchery-wild systems. Fish Sci. 69, 101–109. https://doi.org/10.1046/j.1444-2906.2003.00593.x (2003).CAS 
Article 

Google Scholar 
30.Harada, Y., Yokota, M. & Iizuka, M. Genetic risk of domestication in artificial fish stocking and its possible reduction. Popul. Ecol. 40, 311–324. https://doi.org/10.1007/BF02763463 (1998).Article 

Google Scholar 
31.Etter, P. D., Bassham, S., Hohenlohe, P. A., Johnson, E. A. & Cresko, W. A. SNP discovery and genotyping for evolutionary genetics using RAD sequencing, in: Molecular methods for evolutionary genetics. Methods Mol. Biol. Humana Press 772, 157–178. https://doi.org/10.1007/978-1-61779-228-1_9 (2011).CAS 
Article 

Google Scholar 
32.Rochette, N. C. & Catchen, J. M. Deriving genotypes from RAD-seq short-read data using stacks. Nat. Protoc. 12, 2640–2659. https://doi.org/10.1038/nprot.2017.123 (2017).CAS 
Article 
PubMed 

Google Scholar 
33.Ran, Z. et al. Chromosome-level genome assembly of the razor clam Sinonovacula constricta (Lamarck, 1818). Mol. Ecol. Resour. 19, 1647–1658. https://doi.org/10.1111/1755-0998.13086 (2019).CAS 
Article 
PubMed 

Google Scholar 
34.Li, H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv: 1303.3997. https://arxiv.org/abs/1303.3997 (2013).35.Zheng, X. et al. A high-performance computing toolset for relatedness and principal component analysis of SNP data. Bioinformatics 28, 3326–3328. https://doi.org/10.1093/bioinformatics/bts606 (2012).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
36.Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27, 2156–2158. https://doi.org/10.1093/bioinformatics/btr330 (2011).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
37.Alexander, D. H., Novembre, J. & Lange, K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655–1664. https://pubmed.ncbi.nlm.nih.gov/19648217/ (2009).38.Purcell, S. et al. PLINK: A tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575. https://doi.org/10.1086/519795 (2007).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
39.Eklund, A. R package beeswarm: the bee swarm plot, an alternative to stripchart. https://cran.r-project.org/web/packages/beeswarm/beeswarm.pdf (2016).40.Stamatakis, A. RAxML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313. https://doi.org/10.1093/bioinformatics/btu033 (2014).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
41.Chifman, J. & Kubatko, L. Quartet inference from SNP data under the coalescent model. Bioinformatics 30, 3317–3324. https://doi.org/10.1093/bioinformatics/btu530 (2014).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
42.Swofford, D. L. PAUP*. Phylogenetic analysis using parsimony (*and other methods). Version 4 (Sinauer Associates, 2002).
Google Scholar  More

Species threat abatement and restoration (STAR) metricWe developed and analysed a STAR metric that evaluates the potential benefit for threatened species of actions to reduce threats and restore habitat. Like the Red List Index7,8, STAR is derived from existing data in the IUCN Red List and is intended to help address an urgent need. STAR is spatially explicit, enabling identification of specific threat abatement and habitat restoration opportunities in particular places, which, if implemented, could reduce species extinction risk to levels that would exist without ongoing human impact. Abatement of threats to species encompasses reduction in threat intensity and/or action to mitigate the impacts of threats. Positive population and/or distribution changes, along with the resulting reduction of species extinction risk, have been documented in response to threat abatement13. STAR assumes that, for the great majority of species (see Supplementary Discussion), complete alleviation of threats would reduce extinction risk through halting the decline and/or permitting sufficient recovery in population and distribution, such that the species could be downlisted to the IUCN Red List category of Least Concern. We recognize that complete threat reduction is difficult, incremental conservation gains will need to be tracked at the species level14 and species recovery will vary across a species’ range14.For each species, a global STAR threat abatement (START) score is defined. This varies from zero for species of Least Concern to 100 for Near Threatened, 200 for Vulnerable, 300 for Endangered and 400 for Critically Endangered species (using established weighting ratios7,8). The sum of START values across all species represents the global threat abatement effort needed for all species to become Least Concern. START scores can be disaggregated spatially, based on the area of habitat (AOH) currently available for each species in a particular location (as a proxy for population proportion). This shows the potential contribution of conservation actions in that location to reducing the extinction risk for all species globally. The local START score can be further disaggregated by threat, based on the known contribution of each threat to the species’ risk of extinction (see Methods). This quantifies how actions that abate a specific threat at a particular location contribute to the global abatement of extinction risk for all species.The STAR metric also includes a complementary habitat restoration component to reflect the potential benefits to species of restoring lost habitat. During the United Nations Decade on Ecosystem Restoration (2021–2030), restoration efforts are likely to expand. The STAR restoration component applies a similar logic to the STAR threat abatement component, but for habitat that has been lost and is potentially restorable (that is, restorable AOH). The STAR restoration component does not make assumptions about the extent of habitat restoration required for individual species, but instead quantifies the potential contribution that habitat restoration activities could make to reducing species’ extinction risk. For a particular species at a particular location, the STAR restoration (STARR) score reflects the proportion that restorable habitat at the location represents of the global area of remaining habitat for that species. Importantly, a multiplier is applied to STARR scores to reflect the slower and lower success rate in delivering benefits to species from restored habitat compared with conserved existing habitat15. Again, STARR scores can be disaggregated by threat and summed across species within the location.STAR is intended to provide a metric to underpin the establishment of science-based targets as explicit contributions from individual actors towards the post-2020 biodiversity framework, by allowing assessment of actions and locations according to their potential ability to deliver towards international conservation targets. Individual spatially based START and STARR scores, for all species present in a particular location or country, represent a proportion of the global opportunity to reduce species’ extinction risk through threat abatement and restoration, respectively. For each species, the total START score could be achieved by the complete abatement of all threats in remaining habitat, or an equivalent value of the STAR metric can be achieved by a combination of threat abatement in the remaining habitat and restoration of lost habitat (with concomitant threat abatement therein). The metric can support establishment of science-based targets by a range of actors across spatial scales. By enabling governments and non-state actors to quantify their potential contributions, STAR, along with other tools, could facilitate achievement of global policy goals, notably the species component of the Sustainable Development Goals and the expected post-2020 Global Biodiversity Framework.STAR uses existing publicly available datasets: species’ extinction risk categories and threats available from the IUCN Red List6 (or, for country endemics not yet assessed globally, from national red lists); and species’ AOH estimated using species’ ranges, habitat associations, and elevation limits, along with digital elevation models and current and historical land cover maps (here, we used backcast land cover maps of the distribution of habitat pre-human impact, as in ref. 16). To demonstrate the utility of STAR, we calculated global STAR scores for the groups of terrestrial vertebrate species that are comprehensively assessed on the IUCN Red List (that is, threatened and Near Threatened species of amphibians, birds and mammals globally; n = 5,359).Potential to reduce species extinction riskGlobally, the greatest contribution that could be made to reduce the extinction risk of these groups is tackling threats from annual and perennial non-timber crop production, which account for 24.5% of the global START score (Fig. 1). A further 16.4% is contributed by logging and wood harvesting. There are likely to be specific targets for reducing agriculture and forestry threats in the post-2020 framework3, and applying STAR quantifies the large potential contribution that mitigating these threats could make to the goal for species conservation. Appropriate activities to deliver on such targets range along a continuum from land sharing through to land sparing17.Fig. 1: Contribution to the global START score of different threats and the potential contribution of habitat restoration.The total global START score represents the global threat abatement effort needed for all Near Threatened and threatened (Vulnerable, Endangered and Critically Endangered, according to the IUCN Red List) amphibian, bird and mammal species to be reclassified as Least Concern. This score can be disaggregated by threat type, based on the known contribution of each threat to species’ risk of extinction. The STARR score quantifies the potential contribution that habitat restoration activities could make to reducing overall species’ extinction risk. The total START score could thus be achieved by the complete abatement of all threats in existing natural habitat, or through a combination of threat abatement in existing habitat and restoration of lost habitat (with concomitant threat abatement therein).Full size imageSTAR can be used in combination with existing policy and planning tools to quantify the potential contribution of action targets towards species conservation outcomes. The proposed post-2020 framework includes an action target for the protection of sites of particular importance to biodiversity3. Key Biodiversity Areas11, which include Important Bird and Biodiversity Areas18 and Alliance for Zero Extinction sites19, correspond to such sites. Key Biodiversity Areas so far cover 8.8% of the terrestrial surface (www.keybiodiversityareas.org; identification is ongoing), but already capture 47% of the global START score for the vertebrate groups analysed. They represent large proportions of some national START scores: >70% in Mexico and Venezuela and >50% in Madagascar, Ecuador, the Philippines and Tanzania.START scores can also support target setting at national and sub-national scales to help meet international policy goals. The control and eradication of invasive species forms one of the CBD’s proposed post-2020 action targets3. New Zealand has already set a Predator Free 2050 goal that aims to eradicate three invasive mammal species by 2050. New Zealand contributes 0.8% to the global START score for the three vertebrate groups included in this study. Achieving the Predator Free 2050 goal would contribute 30% of the total START score for New Zealand, amounting to 0.2% of the global START score.All countries contribute towards the global START score, but scores are highly skewed, with a few countries having high START scores and most having low scores for the vertebrate groups analysed (Fig. 2a and Extended Data Fig. 1). The highest-scoring countries are located in biodiverse regions with many threatened endemic species20: Indonesia contributes 7.1% of the global START score, Colombia 7.0%, Mexico 6.1%, Madagascar 6.0% and Brazil 5.2%. These top five countries contribute 31.3% of the global START score. In contrast, the lowest-scoring 88 countries together contribute only 1% of the global START score. This does not imply that these low-scoring countries have negligible species conservation responsibilities; the global decline in even common species indicates that all countries must act to reverse the degradation of nature and restore the diversity and abundance of species and integrity of ecosystems21, as well as preventing extinctions at a national scale. Moreover, most countries have a Red List Index22, or an equivalent, quantifying their progress or failure in reducing the global extinction risk of assessed species relative to their national responsibility for global species conservation. STAR provides a means to guide the reduction of extinction risk and so assist all countries in meeting national species conservation targets.Fig. 2: Global distribution of START and STARR scores.a,b, Global STAR scores for amphibians, birds and mammals at a 50-km grid cell resolution for START scores (a) and STARR scores (b). Each species has a global START score, weighted relative to their extinction risk. This global START score can be disaggregated spatially, based on the AOH currently available for each species in a particular location. The total START score per grid cell (a) is thus the sum of the individual species’ START scores per grid cell across all Near Threatened and threatened species of amphibians, birds and mammals included in this study. The global STARR score per species reflects the potential contribution that habitat restoration activities could make to reducing species’ extinction risk, and is spatially disaggregated based on the availability of restorable habitat. Thus, the total STARR score per grid cell (b) is the sum of the individual species’ STARR scores per grid cell across all species included in this study. For the legends in a and b, each range excludes the lower bound and includes the upper bound.Full size imageAt the global level, we estimated that an equivalent to 55.9% of the global START score for vertebrates could, theoretically, be achieved by restoring lost habitat within the current range (Fig. 1). Ecosystem restoration objectives have been identified in many national biodiversity strategies for the CBD, as well as in many countries’ commitments under the Bonn Challenge, and as part of Nationally Determined Contributions under the United Nations Framework Convention on Climate Change. The STAR metric has the potential to support restoration initiatives alongside species conservation targets by quantifying the potential benefit to particular species of restoring habitat in specific places23 (Fig. 2b). Restoration may be particularly important for some species, including those assessed under Red List sub-criteria D/D1 (with a very small population) or Bac (with a small range with severe fragmentation, plus extreme fluctuations). For species uniquely assessed under these criteria (2.8% of those included in this study), threat abatement alone is unlikely to eliminate extinction risk, so this might need to be complemented by restoration in order to achieve Least Concern status (see Supplementary Discussion). Moreover, depending on habitat loss and threat type, restoration of habitat may be beneficial for a larger proportion of threatened species.Application of STAR at the landscape scaleWe tested the landscape-scale application of the STAR metric in the southern part of Bukit Tigapuluh landscape, in central Sumatra, Indonesia (Fig. 3a). The Bukit Tigapuluh Sustainable Landscape and Livelihoods Project is a sustainable commercial rubber initiative. The study area (approximately 88,000 ha) includes a 5-km buffer (which is set aside to support local livelihoods, wildlife conservation areas and forest protection and restoration) and two ecosystem restoration areas (which form a conservation management zone that protects the Bukit Tigapuluh National Park from encroachment).Fig. 3: STAR results for the Bukit Tigapuluh Sustainable Landscape and Livelihoods Project.The Bukit Tigapuluh Sustainable Landscape and Livelihoods Project is a sustainable commercial rubber initiative. The study area (approximately 88,000 ha) includes a 5-km buffer, which is set aside to support local livelihoods, wildlife conservation areas and forest protection and restoration, and two ecosystem restoration areas, which form a conservation management zone that protects the Bukit Tigapuluh National Park from encroachment. a, Mapped START scores in areas with remaining forest (green) and STARR scores in areas where forest has been lost (purple) at the 30-m grid cell resolution. b, START scores per threat for the top five highest-scoring threats across the study area (the concession, 5-km buffer and ecosystem restoration areas combined).Full size imageThe total START score for the study area represents 0.2% of the START score for Sumatra, 0.04% of the START score for Indonesia and 0.003% of the global START. The major threats are from annual and perennial non-timber crops, logging and wood harvesting, and the collection of terrestrial animals (Fig. 3b). The proximate causes of these pressures in the project area are rubber cultivation, oil palm cultivation, industrial logging, subsistence wood cutting and hunting. STAR analysis shows that areas with the greatest potential to contribute to species conservation through threat mitigation are in remaining natural habitat close to the national park, with a small area of high potential also to the west, where the relatively small distribution of the orbiculus leaf-nosed bat (Hipposideros orbiculus) overlaps the site (Fig. 3a). Additionally, due to recent forest loss, 47% of the START score for the study area could be achieved through habitat restoration (that is, STARR). Investment in these management actions has the potential to deliver these quantified contributions to national and global biodiversity targets.Operationalization and future developmentThe STAR metric makes use of the best available data, producing results that are relevant to policy and practice. However, there is scope for future refinement as the underlying data improve. Here, the STAR metric covers amphibians, birds and mammals globally, constituting a well-studied but small proportion of taxonomic diversity (see Extended Data Figs. 2 and 3 for variation among taxa). STAR can be expanded to other taxonomic groups, including freshwater and marine species, as data become available (reptiles, cacti, cycads, conifers, freshwater fish and reef-building corals are among the groups imminently available for incorporation). Global application of STAR will require comprehensive assessment of taxonomic groups, testing of the transferability of the STAR metric assumptions among taxa as Red List coverage expands, and further development of methods to calculate AOH. AOH calculation does not currently capture spatial variation in species’ population density, which will be important for many species14; such data have not been gathered on a global scale yet and could be incorporated as available.The completeness of threat data in the IUCN Red List is uneven but is continually improving. The STAR metric does not currently reflect spatial variation in threat magnitude within species’ ranges; more broadly, there is a lack of information on the spatial distribution of threats24. Most species included in this study have relatively small ranges; the total current AOH is More