Philippa Kaur

1.
Giljum, S., Dittrich, M., Lieber, M. & Lutter, S. Global patterns of material flows and their socio-economic and environmental implications: A MFA study on all countries world-wide from 1980 to 2009. Resources 3, 319–339 (2014).
Article  Google Scholar 
2.
IRP, U. Global Resources Outlook 2019: Natural Resources for the Future we Want. A Report of the International Resource Panel. Report No. DTI/2226/NA (United Nations Environment Programme, 2019).

3.
Krausmann, F., Schandl, H., Eisenmenger, N., Giljum, S. & Jackson, T. Material flow accounting: Measuring global material use for sustainable development. Ann. Rev. Env. Resour. 42, 647–675 (2017).
Article  Google Scholar 

4.
Calvo, G., Mudd, G., Valero, A. & Valero, A. Decreasing ore grades in global metallic mining: A theoretical issue or a global reality? Resources 5 (2016).

5.
Prior, T., Giurco, D., Mudd, G., Mason, L. & Behrisch, J. Resource depletion, peak minerals and the implications for sustainable resource management. Glob. Environ. Change 22, 577–587 (2012).
Article  Google Scholar 

6.
West, J. Decreasing metal ore grades. J. Ind. Ecol. 15, 165–168 (2011).
Article  Google Scholar 

7.
Mudd, G. M. Global trends in gold mining: Towards quantifying environmental and resource sustainability. Resour. Policy 32, 42–56 (2007).
Article  Google Scholar 

8.
Sonter, L. J., Moran, C. J., Barrett, D. J. & Soares-Filho, B. S. Processes of land use change in mining regions. J. Clean. Prod. 84, 494–501 (2014).
Article  Google Scholar 

9.
Werner, T., Bebbington, A. & Gregory, G. Assessing impacts of mining: Recent contributions from GIS and remote sensing. Extract. Indus. Soc. 6, 993–1012 (2019).
Article  Google Scholar 

10.
Kobayashi, H., Watando, H. & Kakimoto, M. A global extent site-level analysis of land cover and protected area overlap with mining activities as an indicator of biodiversity pressure. J. Clean. Prod. 84, 459–468 (2014).
Article  Google Scholar 

11.
Sonter, L. J., Ali, S. H. & Watson, J. E. M. Mining and biodiversity: key issues and research needs in conservation science. Proc. Biol. Sci. 285 (2018).

12.
Islam, K., Vilaysouk, X. & Murakami, S. Integrating remote sensing and life cycle assessment to quantify the environmental impacts of copper-silver-gold mining: A case study from laos. Resour. Conserv. Recy. 154, 104630 (2020).
Article  Google Scholar 

13.
Butt, N. et al. Biodiversity risks from fossil fuel extraction. Science 342, 425–426 (2013).
ADS  CAS  Article  Google Scholar 

14.
Murguía, D. I., Bringezu, S. & Schaldach, R. Global direct pressures on biodiversity by large-scale metal mining: Spatial distribution and implications for conservation. J. Eenviron. Manage. 180, 409–420 (2016).
Google Scholar 

15.
Endl, A., Tost, M., Hitch, M., Moser, P. & Feiel, S. Europe’s mining innovation trends and their contribution to the sustainable development goals: Blind spots and strong points. Resour. Policy 101440 (2019).

16.
Bruckner, M., Fischer, G., Tramberend, S. & Giljum, S. Measuring telecouplings in the global land system: A review and comparative evaluation of land footprint accounting methods. Ecol. Econ. 114, 11–21 (2015).
Article  Google Scholar 

17.
Schaffartzik, A. et al. Trading land: A review of approaches to accounting for upstream land requirements of traded products. J. Ind. Ecol. 19, 703–714 (2015).
Article  Google Scholar 

18.
USGS – United States Geological Survey. Mineral resources online spatial data, https://mrdata.usgs.gov/ (2018).

19.
S&P Global Market Intelligence. SNL metals and mining database, https://www.spglobal.com/marketintelligence/en/campaigns/metals-mining (2018).

20.
Murguía, D. I. & Bringezu, S. Measuring the specific land requirements of large-scale metal mines for iron, bauxite, copper, gold and silver. Prog. Ind. Ecol. 10, 264–285 (2016).
Article  Google Scholar 

21.
Werner, T. T. et al. Global-scale remote sensing of mine areas and analysis of factors explaining their extent. Glob. Environ. Change 60 (2020).

22.
Mountrakis, G., Im, J. & Ogole, C. Support vector machines in remote sensing: A review. ISPRS J. Photogramm. 66, 247–259 (2011).
Article  Google Scholar 

23.
Belgiu, M. & Dragu, L. Random forest in remote sensing: A review of applications and future directions. ISPRS J. Photogramm. 114, 24–31 (2016).
Article  Google Scholar 

24.
Zhu, X. X. et al. Deep learning in remote sensing: A comprehensive review and list of resources. IEEE Geosc. Rem. Sen. M. 5, 8–36 (2017).
Article  Google Scholar 

25.
Wulder, M. A., Coops, N. C., Roy, D. P., White, J. C. & Hermosilla, T. Land cover 2.0. Int. J. Remote Sens. 39, 4254–4284 (2018).
ADS  Article  Google Scholar 

26.
Zhu, Z. et al. Benefits of the free and open Landsat data policy. Remote Sens. Environ. 224, 382–385 (2019).
ADS  Article  Google Scholar 

27.
Petropoulos, G. P., Partsinevelos, P. & Mitraka, Z. Change detection of surface mining activity and reclamation based on a machine learning approach of multi-temporal Landsat TM imagery. Geocarto Int. 28, 323–342 (2013).
Article  Google Scholar 

28.
LaJeunesse Connette, K. J. et al. Assessment of mining extent and expansion in Myanmar based on freely-available satellite imagery. Remote Sens. 8 (2016).

29.
Yu, L. et al. Monitoring surface mining belts using multiple remote sensing datasets: A global perspective. Ore Geol. Rev. 101, 675–687 (2018).
Article  Google Scholar 

30.
Vasuki, Y. et al. The spatial-temporal patterns of land cover changes due to mining activities in the darling range, western australia: A visual analytics approach. Ore Geol. Rev. 108, 23–32 (2019).
Article  Google Scholar 

31.
Mukherjee, J., Mukherjee, J., Chakravarty, D. & Aikat, S. A novel index to detect opencast coal mine areas from Landsat 8 OLI/TIRS. IEEE J-STARS 12, 891–897 (2019).
Google Scholar 

32.
Waldrop, M. M. News Feature: What are the limits of deep learning? PNAS 116, 1074–1077 (2019).
CAS  Article  Google Scholar 

33.
EOX IT Services GmbH. Sentinel-2 cloudless (contains modified Copernicus sentinel data 2017 and 2018), https://s2maps.eu (2018).

34.
Pebesma, E. Simple Features for R: Standardized Support for Spatial Vector Data. R J. 10, 439–446 (2018).
Article  Google Scholar 

35.
Gutschlhofer, J. & Maus, V. Web application for mining area polygonization version 1.2. Zenodo https://doi.org/10.5281/zenodo.3691743 (2020).

36.
Lesiv, M. et al. Characterizing the spatial and temporal availability of very high resolution satellite imagery in Google Earth and Microsoft Bing maps as a source of reference data. Land 7 (2018).

37.
Bradshaw, A. Restoration of mined lands—using natural processes. Ecol. Eng. 8, 255–269 (1997).
Article  Google Scholar 

38.
EUROSTAT. Countries, 2016 – administrative units – dataset (generalised dataset derived from eurogeographics and UN-FAO GI data), https://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/countries/ (2018).

39.
Amatulli, G. et al. A suite of global, cross-scale topographic variables for environmental and biodiversity modeling. Sci. Data 5, 180040 (2018).
Article  Google Scholar 

40.
Maus, V. et al. Global-scale mining polygons (version 1). Pangaea https://doi.org/10.1594/PANGAEA.910894 (2020).

41.
Marazuela, M., Vázquez-Suñé, E., Ayora, C., García-Gil, A. & Palma, T. The effect of brine pumping on the natural hydrodynamics of the Salar de Atacama: The damping capacity of salt flats. Sci. Total Environ. 654, 1118–1131 (2019).
ADS  CAS  Article  Google Scholar 

42.
Liu, W., Agusdinata, D. B. & Myint, S. W. Spatiotemporal patterns of lithium mining and environmental degradation in the Atacama Salt Flat, Chile. Int. J. Appl. Earth Obs. 80, 145–156 (2019).
Article  Google Scholar 

43.
Hansen, K. Brazil’s Carajás mines, NASA Earth Observatory, https://earthobservatory.nasa.gov/images/144457/brazils-carajas-mines (2018).

44.
Mining Technology. Batu Hijau copper-gold mine, Indonesia, https://www.mining-technology.com/projects/batu/ (2020).

45.
Shen, L. & Gunson, A. J. The role of artisanal and small-scale mining in China’s economy. J. Clean. Prod. 14, 427–435 (2006).
Article  Google Scholar 

46.
Shen, L., Dai, T. & Gunson, A. J. Small-scale mining in China: Assessing recent advances in the policy and regulatory framework. Resour. Policy 34, 150–157 (2009).
Article  Google Scholar 

47.
Potere, D. Horizontal positional accuracy of Google Earth’s high-resolution imagery archive. Sensors 8, 7973–7981 (2008).
Article  Google Scholar 

48.
Vajsová B & Åstrand, P. J. New sensors benchmark report on Sentinel-2A sensor over Maussane test site for CAP purposes. Report No. EUR 27674EN (Publications Office of the European Union, 2015).

49.
Cochran, W. G. Sampling Techniques. Series in Probability and Statistics (Wiley, 1977), 3 edn.

50.
Olofsson, P. et al. Good practices for estimating area and assessing accuracy of land change. Remote Sens. Environ. 148, 42–57 (2014).
ADS  Article  Google Scholar 

51.
OGC – Open Geospatial Consortium. GeoPackage Encoding Standard, https://www.geopackage.org/ (2005).

52.
OGC – Open Geospatial Consortium. Geographic tagged image file format (GeoTIFF), https://www.ogc.org/standards/geotiff (2019).

53.
The Internet Society. RFC 4180: Common format and MIME type for comma-separated values (CSV). https://tools.ietf.org/html/rfc4180 (2005).

54.
Wang, J.-F., Zhang, T.-L. & Fu, B.-J. A measure of spatial stratified heterogeneity. Ecol. Indic. 67, 250–256 (2016).
Article  Google Scholar 

55.
Brunsdon, C., Fotheringham, A. S. & Charlton, M. E. Geographically weighted regression: A method for exploring spatial nonstationarity. Geogr. Anal. 28, 281–298 (1996).
Article  Google Scholar 

56.
Brunsdon, C., Fotheringham, S. & Charlton, M. Geographically weighted regression. J. R. Stat. Soc., Ser. D Stat. 47, 431–443 (1998).
Article  Google Scholar 

57.
QGIS Development Team. QGIS geographic information system, version 3.12.0. Open Source Geospatial Foundation, https://www.qgis.org (2020).

58.
R Core Team. R: A language and environment for statistical computing, version 3.6.1. Foundation for Statistical Computing, Vienna, Austria, https://www.R-project.org (2019).

59.
Python Core Team. Python: A dynamic, open source programming language, version 2.7.17. Python Software Foundation, https://www.python.org (2019).

60.
OGC – Open Geospatial Consortium. Web map service interface standard (WMS), https://www.ogc.org/standards/wms (2020).

61.
GNU general public license, version 3. Free Software Foundation, https://www.gnu.org/licenses/gpl-3.0.en.html (2019).

62.
GDAL/OGR contributors. GDAL/OGR geospatial data abstraction software library, version 2.4.2. Open Source Geospatial Foundation, https://gdal.org (2019).

63.
Chang, W., Cheng, J., Allaire, J., Xie, Y. & McPherson, J. Shiny: Web Application Framework for R, version 1.3.2, https://CRAN.R-project.org/package=shiny (2019)

64.
The PostgreSQL Global Development Group. PostgreSQl: an open source object-relational database system, version 11.6, https://www.postgresql.org/ (2019).

65.
PostGIS Team. PostGIS: a spatial database extender for PostgreSQL object relational database, version 2.5.4. Open Source Geospatial Foundation, https://postgis.net (2019). More

Adams HD, Collins AD, Briggs SP, Vennetier M, Dickman LT, Sevanto SA et al. (2015) Experimental drought and heat can delay phenological development and reduce foliar and shoot growth in semiarid trees. Glob Change Biol 21:4210–4220
Google Scholar 

Basler D, Körner C (2012) Photoperiod sensitivity of bud burst in 14 temperate forest tree species. Agric For Meteorol 165:73–81
Google Scholar 

Bhalerao R, Keskitalo J, Sterky F, Erlandsson R, Björkbacka H, Birve SJ et al. (2003) Gene expression in autumn leaves. Plant Physiol 131:430–442
PubMed  PubMed Central  Google Scholar 

Bradshaw HD, Stettler RF (1995) Molecular genetics of growth and development in populus. IV. Mapping QTLs with large effects on growth, form, and phenology traits in a forest tree. Genetics 139:963–973
CAS  PubMed  Google Scholar 

Brelsford CC, Nybakken L, Kotilainen TK, Robson TM (2019) The influence of spectral composition on spring and autumn phenology in trees. Tree Physiol 39:925–950
CAS  PubMed  Google Scholar 

Christensen K, Frederiksen H, Vaupel JW, McGue M (2003) Age trajectories of genetic variance in physical functioning: a longitudinal study of Danish twins aged 70 years and older. Behav Genet 33:125–136
PubMed  Google Scholar 

Chuine I (2010) Why does phenology drive species distribution? Philos Trans R Soc B 365:3149–3160
Google Scholar 

Class B, Brommer JE (2020) Can dominance genetic variance be ignored in evolutionary quantitative genetic analyses of wild populations? Evolution 74:1540–1550
PubMed  Google Scholar 

Cong N, Shen M, Piao S (2016) Spatial variations in responses of vegetation autumn phenology to climate change on the Tibetan Plateau. J Plant Ecol 10:744–752.
Google Scholar 

Costa e Silva J, Borralho NMG, Potts BM (2004) Additive and non-additive genetic parameters from clonally replicated and seedling progenies of Eucalyptus globulus. Theor Appl Genet 108:1113–1119
PubMed  Google Scholar 

Cronk QCB (2005) Plant eco‐devo: the potential of poplar as a model organism. N Phytologist 166:39–48
CAS  Google Scholar 

Evans LM, Slavov GT, Rodgers-Melnick E, Martin J, Ranjan P, Muchero W et al. (2014) Population genomics of Populus trichocarpa identifies signatures of selection and adaptive trait associations. Nat Genet 46:1089
CAS  PubMed  Google Scholar 

Fabbrini F, Gaudet M, Bastien C, Zaina G, Harfouche A, Beritognolo I et al. (2012) Phenotypic plasticity, QTL mapping and genomic characterization of bud set in black poplar. BMC Plant Biol 12:47
CAS  PubMed  PubMed Central  Google Scholar 

Falconer DS, Mackay TFC (1996) Introduction to quantitative genetics. 4th edn. Longman

Fracheboud Y, Luquez V, Björkén L, Sjödin A, Tuominen H, Jansson S (2009) The control of autumn senescence in European Aspen. Plant Physiol 149:1982–1991
CAS  PubMed  PubMed Central  Google Scholar 

Hadfield JD (2010) MCMC methods for multi-response generalized linear mixed models: the MCMCglmm R package. J Stat Softw 33:1–22.
Google Scholar 

Horvath DP, Anderson JV, Chao WS, Foley ME (2003) Knowing when to grow: signals regulating bud dormancy. Trends Plant Sci 8:534–540
CAS  PubMed  Google Scholar 

Howe GT, Aitken SN, Neale DB, Jermstad KD, Wheeler NC, Chen THH (2003) From genotype to phenotype: unraveling the complexities of cold adaptation in forest trees. Can J Bot 81:1247–1266
CAS  Google Scholar 

Howe GT, Saruul P, Davis J, Chen THH (2000) Quantitative genetics of bud phenology, frost damage, and winter survival in an F2 family of hybrid poplars. Theor Appl Genet 101:632–642
Google Scholar 

Ingvarsson PK, Garcia MV, Luquez V, Hall D, Jansson S (2008) Nucleotide polymorphism and phenotypic associations Within and Around the phytochrome B2 Locus in European Aspen (Populus tremula, Salicaceae). Genetics 178:2217
CAS  PubMed  PubMed Central  Google Scholar 

Jansson S, Douglas CJ (2007) Populus: a model system for plant biology. Annu Rev Plant Biol 58:435–458
CAS  PubMed  Google Scholar 

Karacic A, Verwijst T, Weih M (2003) Above-ground woody biomass production of short-rotation populus plantations on agricultural land in Sweden. Scand J For Res 18:427–437
Google Scholar 

Kennedy BW, Schaeffer LR (1989) Genetic evaluation under an animal model when identical genotypes are represented in the population. J Anim Sci 67:1946–1955
Google Scholar 

Keskitalo J, Bergquist G, Gardeström P, Jansson S (2005) A cellular timetable of autumn senescence. Plant Physiol 139:1635–1648
CAS  PubMed  PubMed Central  Google Scholar 

Lagercrantz ULF (2009) At the end of the day: a common molecular mechanism for photoperiod responses in plants? J Exp Bot 60:2501–2515
CAS  PubMed  Google Scholar 

Leuchner M, Menzel A, Werner H (2007) Quantifying the relationship between light quality and light availability at different phenological stages within a mature mixed forest. Agric For Meteorol 142:35–44
Google Scholar 

Luquez V, Hall D, Albrectsen BR, Karlsson J, Ingvarsson P, Jansson S (2008) Natural phenological variation in aspen (Populus tremula): the SwAsp collection. Tree Genet Genomes 4:279–292
Google Scholar 

Lynch M, Walsh B (1998) Genetics and analysis of quantitative traits. Sinauer Sunderland

MacKenzie CM, Primack RB, Miller-Rushing AJ (2018) Local environment, not local adaptation, drives leaf-out phenology in common gardens along an elevational gradient in Acadia National Park, Maine. Am J Bot 105:986–995
Google Scholar 

Marron N, Storme V, Dillen SY, Bastien C, Ricciotti L, Salani F et al. (2010) Genomic regions involved in productivity of two interspecific poplar families in Europe. 2. Biomass production and its relationships with tree architecture and phenology. Tree Genet Genomes 6:533–554
Google Scholar 

McKown AD, Guy RD, Klápště J, Geraldes A, Friedmann M, Cronk QC et al. (2014a) Geographical and environmental gradients shape phenotypic trait variation and genetic structure in Populus trichocarpa. N. Phytologist 201:1263–1276
CAS  Google Scholar 

McKown AD, Klápště J, Guy RD, Geraldes A, Porth I, Hannemann J et al. (2014b) Genome‐wide association implicates numerous genes underlying ecological trait variation in natural populations of Populus trichocarpa. N. Phytologist 203:535–553
CAS  Google Scholar 

McKown AD, Klápště J, Guy RD, El-Kassaby YA, Mansfield SD (2018) Ecological genomics of variation in bud-break phenology and mechanisms of response to climate warming in Populus trichocarpa. N. Phytologist 220:300–316
CAS  Google Scholar 

Michelson I, Ingvarsson Pär K, Robinson Kathryn M, Edlund E, Eriksson Maria E, Nilsson O et al. (2017) Autumn senescence in aspen is not triggered by day length. Physiol Plant 162:123–134
PubMed  Google Scholar 

Olson M, Levsen N, Soolanayakanahally Raju Y, Guy Robert D, Schroeder William R, Keller Stephen R et al. (2013) The adaptive potential of Populus balsamifera L. to phenology requirements in a warmer global climate. Mol Ecol 22:1214–1230
CAS  PubMed  Google Scholar 

Pliura A, Suchockas V, Sarsekova D, Gudynaitė V (2014) Genotypic variation and heritability of growth and adaptive traits, and adaptation of young poplar hybrids at northern margins of natural distribution of Populus nigra in Europe. Biomass Bioenergy 70:513–529
Google Scholar 

Polgar CA, Primack RB (2011) Leaf‐out phenology of temperate woody plants: from trees to ecosystems. N Phytologist 191:926–941
Google Scholar 

Porth I, El‐Kassaby YA (2015) Using Populus as a lignocellulosic feedstock for bioethanol. Biotechnol J 10:510–524
CAS  PubMed  Google Scholar 

Porth I, Klápště J, McKown AD, La Mantia J, Guy RD, Ingvarsson PK et al. (2015) Evolutionary quantitative genomics of Populus trichocarpa. PLoS ONE 10:e0142864
PubMed  PubMed Central  Google Scholar 

Porth I, Klápště J, McKown AD, La Mantia J, Hamelin RC, Skyba O et al. (2014) Extensive functional pleiotropy of REVOLUTA substantiated through forward genetics. Plant Physiol 164:548–554
CAS  PubMed  Google Scholar 

Rohde A, Bastien C, Boerjan W (2011b) Temperature signals contribute to the timing of photoperiodic growth cessation and bud set in poplar. Tree Physiol 31:472–482
PubMed  Google Scholar 

Rohde A, Storme V, Jorge V, Gaudet M, Vitacolonna N, Fabbrini F et al. (2011a) Bud set in poplar – genetic dissection of a complex trait in natural and hybrid populations. N Phytologist 189:106–121
CAS  Google Scholar 

Sannigrahi P, Ragauskas AJ, Tuskan GA (2010) Poplar as a feedstock for biofuels: a review of compositional characteristics. Biofuels Bioprod Bioref 4:209–226
CAS  Google Scholar 

Singh RK, Svystun T, AlDahmash B, Jönsson AM, Bhalerao RP (2017) Photoperiod- and temperature-mediated control of phenology in trees—a molecular perspective. N Phytologist 213:511–524
CAS  Google Scholar 

Slavov GT, Leonardi S, Burczyk J, Adams WT, Strauss Sh, DiFazio SP (2009) Extensive pollen flow in two ecologically contrasting populations of Populus trichocarpa. Mol Ecol 18:357–373

Triozzi PM, Ramos-Sánchez JM, Hernández-Verdeja T, Moreno-Cortés A, Allona I, Perales M (2018) Photoperiodic regulation of shoot apical growth in poplar. Front Plant Sci 9:1030–1030
PubMed  PubMed Central  Google Scholar 

Vitasse Y, Delzon S, Dufrêne E, Pontailler J-Y, Louvet J-M, Kremer A et al. (2009) Leaf phenology sensitivity to temperature in European trees: Do within-species populations exhibit similar responses? Agric For Meteorol 149:735–744
Google Scholar 

Walsh B, Lynch M (2018) Evolution and selection of quantitative traits. Oxford University Press

Weih M (2004) Intensive short rotation forestry in boreal climates: present and future perspectives. Can J For Res 34:1369–1378
Google Scholar 

Wilson AJ, Kruuk LEB, Coltman DW (2005) Ontogenetic patterns in heritable variation for body size: using random regression models in a wild ungulate population. Am Naturalist 166:E177–E192
Google Scholar 

Yu Q, Tigerstedt PMA, Haapanen M (2001) Growth and phenology of hybrid aspen clones (Populus tremula L. x Populus tremuloides Michx.)

Zanewich KP, Pearce DW, Rood SB (2018) Heterosis in poplar involves phenotypic stability: cottonwood hybrids outperform their parental species at suboptimal temperatures. Tree Physiol 38:789–800
CAS  PubMed  Google Scholar 

Zhou L, Bawa R, Holliday JA (2014) Exome resequencing reveals signatures of demographic and adaptive processes across the genome and range of black cottonwood (Populus trichocarpa). Mol Ecol 23:2486–2499
CAS  PubMed  Google Scholar  More

1.
Schminke, H. K. Entomology for the copepodologist. J. Plankton Res. 29, i149–i162, https://doi.org/10.1093/plankt/fbl073 (2007).
Article  Google Scholar 
2.
Le Borgne, R. Equivalences between the measures of biovolumes, dry weight, ashfree dry weight, carbon, nitrogen and phosphorus of the mesozooplankton of the tropical Atlantic. Cah. O.R.S.T.O.M. Ser. Oceanogr. 13, 179–196 (1975).
Google Scholar 

3.
Corral Estrada, J. Contribucion al conocimiento del plancton de Canarias. Estudio cuantitativo, sistematico y observaciones ecologicas de nos Copépodos epipelagicos en la zona de Santa Cruz de Tenerife en el curso de un ciclo anual. Publnes Fac. Cienc. Madrid, Seccion de Biol., (A) 129 Doct. thesis, Univ. Madrid (1970).

4.
Lovegrove, T. In Some contemporary studies in marine science (ed Barnes, H.) 429-467 (Allen and Unwin, 1966).

5.
Harris, R., Wiebe, P., Lenz, J., Skjoldal, H. R. & Huntley, M. ICES zooplankton methodology manual. (Elsevier, 2000).

6.
McEnnulty, F. R. et al. The Australian Zooplankton Biomass Database (1932–2019). Australian Ocean Data Network https://doi.org/10.26198/5c4170d42ab24 (2019).

7.
Davies, C. H. et al. Over 75 years of zooplankton data from Australia. Ecology 95, 3229–3229, https://doi.org/10.1890/14-0697.1 (2014).
Article  Google Scholar 

8.
Davies, C. H. et al. A database of marine phytoplankton abundance, biomass and species composition in Australian waters (vol 3, 160043, 2016). Sci. Data 3, https://doi.org/10.1038/sdata.2016.111.

9.
Davies, C. H. et al. A database of chlorophyll a in Australian waters. Sci. Data 5, https://doi.org/10.1038/sdata.2018.18

10.
Skerratt, J. H. et al. Simulated nutrient and plankton dynamics in the Great Barrier Reef (2011–2016). J. Mar. Syst. 192, 51–74, https://doi.org/10.1016/j.jmarsys.2018.12.006 (2019).
Article  Google Scholar 

11.
Bernardi, D. In A manual on methods for the assessment of secondary productivity in freshwaters (eds Downing, J.A & Rigler, F.H.) 59–86 (Blackwell Scientific, Oxford, 1984).

12.
Omori, M. & Ikeda, T. Methods in marine zooplankton ecology. (Wiley, 1984).

13.
Richardson, A. J. et al. Using continuous plankton recorder data. Prog. Oceanogr. 68, 27–74, https://doi.org/10.1016/j.pocean.2005.09.011 (2006).
ADS  Article  Google Scholar 

14.
Everett, J. D. et al. Modeling What We Sample and Sampling What We Model: Challenges for Zooplankton Model Assessment. Front. Mar. Sci. 4, https://doi.org/10.3389/fmars.2017.00077 (2017).

15.
Herman, A. W., Beanlands, B. & Phillips, E. F. The next generation of Optical Plankton Counter: The Laser-OPC. J. Plankton Res. 26, 1135–1145, https://doi.org/10.1093/plankt/fbh095 (2004).
Article  Google Scholar 

16.
Bivand, R. & Lewin-Koh, N. maptools: Tools for Handling Spatial Objects. R package version 0.9–5, https://CRAN.R-project.org/package=maptools. (2019).

17.
Amante, C. & Eakins, B. W. ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis. NOAA Technical Memorandum NESDIS. NGDC-24, https://doi.org/10.7289/V5C8276M.19 (NOAA., National Geophysical Data Center, 2009).

18.
Pante, E. & Simon-Bouhet, B. marmap: A Package for Importing, Plotting and Analyzing Bathymetric and Topographic Data in R. Plos One 8, https://doi.org/10.1371/journal.pone.0073051 (2013).

19.
Wiebe, P. H. & Benfield, M. C. From the Hensen net toward four-dimensional biological oceanography. Prog. Oceanogr. 56, 7–136, https://doi.org/10.1016/s0079-6611(02)00140-4 (2003).
ADS  Article  Google Scholar 

20.
Smith, P. E., Counts, R. C. & Clutter, R. I. Changes in the filtering efficiency of plankton nets due to clogging under tow. Cons. Perm. Int. p. l’Expl. d. l. Mer 32, 232–248 (1968).
Article  Google Scholar 

21.
Jackson, C. J., Rothlisberg, P. C. & Pendrey, R. C. Role of larval distribution and abundance in overall life-history dynamics: a study of the prawn Penaeus semisulcatus in Albatross Bay, Gulf of Carpentaria, Australia. Mar. Ecol. Prog. Ser. 213, 241–252, https://doi.org/10.3354/meps213241 (2001).
ADS  Article  Google Scholar 

22.
Hernroth, L. Sampling and filtration efficiency of 2 commonly used plankton nets – A comparative study of the Nansen net and the UNESCO WP-2 net. J. Plankton Res. 9, 719–728, https://doi.org/10.1093/plankt/9.4.719 (1987).
Article  Google Scholar 

23.
Verheye, H. M., Richardson, A. J., Hutchings, L., Marska, G. & Gianakouras, D. Long-term trends in the abundance and community structure of coastal zooplankton in the southern Benguela system, 1951–1996. S. Afr. J. Mar. Sci. 19, 317–332 (1998).
Article  Google Scholar 

24.
Colton, J. B., Green, J. R., Byron, R. R. & Frisella, J. L. Bongo net retention rates as effected by towing speed and mesh size. Can. J. Fish. Aquat.Sci. 37, 606–623, https://doi.org/10.1139/f80-077 (1980).
Article  Google Scholar 

25.
Pillar, S. C. A comparison of the performance of four zooplankton samplers. S. Afr. J. Mar. Sci. 2, 1–18 (1984).
Article  Google Scholar 

26.
McKinnon, A. D. et al. Zooplankton Growth, Respiration and Grazing on the Australian Margins of the Tropical Indian and Pacific Oceans. Plos One 10, e0140012, https://doi.org/10.1371/journal.pone.0140012 (2015).
CAS  Article  PubMed  PubMed Central  Google Scholar 

27.
Beers, J. R. Determination of zooplankton biomass. Monogr. Oceanogr. Methodol. 4, 35–84 (1976).
Google Scholar 

28.
Postel, L. et al. 83–192 (Elsevier, 2000).

29.
Motoda, S. Plankton sampler for collecting uncontaminated materials from several zones by a single vertical haul. Rapp. P.-V. Réun. Cons. Int. Explor. Mer. 153, 55–58 (1962).
Google Scholar 

30.
Motoda, S. & Osawa, K. Filteration ratio, variance of samples and estimated distance of haul in vertical hauls with Indian Ocean standard net. Inf. Bull. Planktol. Japan 11, 11–24 (1964).
Google Scholar 

31.
Wiebe, P. H., Boyd, S. & Cox, J. L. Relationships between zooplankton displacement volume, wet weight, dry weight, and carbon. Fish. Bull. (Wash. D. C.) 73, 777–796 (1975).
Google Scholar 

32.
Wiebe, P. H. Functional regression equations for zooplankton displacement volume, wet weight, dry weight, and carbon: a correction. U S Fish. Wildl. Serv. Fish. Bull. 86, 833–835 (1988).
Google Scholar 

33.
Moriarty, R. & O’Brien, T. D. Distribution of mesozooplankton biomass in the global ocean. Earth Syst. Sci. Data 5, 45–55, https://doi.org/10.5194/essd-5-45-2013 (2013).
ADS  Article  Google Scholar 

34.
Omori, M. Some factors affecting on dry-weight, organic weight and concentrations of carbon and nitrogen in freshlt prepared and in preserved zooplankton. Int. Rev. Gesamt. Hydrobiol. 63, 261–269, https://doi.org/10.1002/iroh.19780630211 (1978).
Article  Google Scholar 

35.
Balvay, P. G. Equivalence entre quelques paramètres estimatifs de l’abondance du zooplancton total. Swiss J. Hydrol. 49, 75–84, https://doi.org/10.1007/BF02540381 (1987).
Article  Google Scholar 

36.
Tranter, D. J. Comparison of zooplankton biomass determinations by Indian Ocean Standard Net, Juday Net and Clarke-Bumpus sampler. Nature 198, 1179, https://doi.org/10.1038/1981179a0 (1963).
ADS  Article  Google Scholar 

37.
DeVries, D. R. & Stein, R. A. Comparison of three zooplankton samplers: a taxon-specific assessment. J. Plankton Res. 13, 53–59, https://doi.org/10.1093/plankt/13.1.53 (1991).
Article  Google Scholar 

38.
Skjoldal, H. R. et al. Intercomparison of zooplankton (net) sampling systems: Results from the ICES/GLOBEC sea-going workshop. Prog. Oceanogr. 108, 1–42, https://doi.org/10.1016/j.pocean.2012.10.006 (2013).
ADS  Article  Google Scholar 

39.
Eriksen, R. S. et al. Australia’s Long-Term Plankton Observations: The Integrated Marine Observing System National Reference Station Network. Front. Mar. Sci. 6, https://doi.org/10.3389/fmars.2019.00161 (2019).

40.
O’Brien, T. D. COPEPOD: The Global Plankton Database. An overview of the 2014 database contents, processing methods, and access interface. U.S. Dep. Commerce, NOAA Tech. Memo., NMFS-F/ST-37. 29 p. (NOAA, 2014).

41.
O’Brien, T. D. COPEPOD: The Global Plankton Database. An overview of the 2010 database contents, processing methods, and access interface. U.S. Dep. Commerce, NOAA Tech. Memo., NMFS-F/ST-36. 28 p. (NOAA, 2010).

42.
O’Brien, T. D. COPEPOD: The Global Plankton Database. A review of the 2007 database contents and new quality control methodology. U.S. Dep. Commerce, NOAA Tech. Memo., NMFS-F/ST-34. 28 p. (NOAA, 2007).

43.
O’Brien, T. D. COPEPOD: A Global Plankton Database. U.S. Dep. Commerce, NOAA Tech. Memo. NMFS-F/SPO-73. 136 p. (NOAA, 2002).

44.
O’Brien, T. D., Wiebe, P. H. & Falkenhaug, T. ICES Zooplankton Status Report 2010/2011. ICES Cooperative Research Report. 318, 208 (2013).
Google Scholar 

45.
Clark, R. A., Frid, C. L. J. & Batten, S. A critical comparison of two long-term zooplankton time series from the central-west North Sea. J. Plankton Res. 23, 27–39 (2001).
Article  Google Scholar 

46.
John, E. H., Batten, S. D., Harris, R. P. & Hays, G. C. Comparison between zooplankton data collected by the Continuous Plankton Recorder survey in the English Channel and by WP-2 nets at station L4, Plymouth (UK). J. Sea Res. 46, 223–232 (2001).
ADS  Article  Google Scholar 

47.
Pitois, S. G. & Fox, C. J. Long-term changes in zooplankton biomass concentration and mean size over the Northwest European shelf inferred from Continuous Plankton Recorder data. ICES J. Mar. Sci. 63, 785–798, https://doi.org/10.1016/j.iccsjms.2006.03.009 (2006).
CAS  Article  Google Scholar 

48.
Richardson, A. J., John, E. H., Irigoien, X., Harris, R. P. & Hays, G. C. How well does the continuous plankton recorder (CPR) sample zooplankton? A comparison with the Longhurst Hardy Plankton Recorder (LHPR) in the northeast Atlantic. Deep Sea Res. Part I. 51, 1283–1294, https://doi.org/10.1016/j.dsr.2004.04.002 (2004).

49.
Herman, A. W. Design and calibration of a new optical plankton counter capable of sizing small zooplankton. Deep Sea Res. Part I. 39, 395–415, https://doi.org/10.1016/0198-0149(92)90080-D (1992).
ADS  Article  Google Scholar 

50.
Espinasse, B. et al. Conditions for assessing zooplankton abundance with LOPC in coastal waters. Prog. Oceanogr. 163, 260–270, https://doi.org/10.1016/j.pocean.2017.10.012 (2018).
Article  Google Scholar 

51.
Baird, M. E. et al. Biological properties across the Tasman Front off southeast Australia. Deep Sea Res. Part I. 55, 1438–1455, https://doi.org/10.1016/j.dsr.2008.06.011 (2008).
Article  Google Scholar 

52.
Everett, J. D., Baird, M. E. & Suthers, I. M. Three-dimensional structure of a swarm of the salp Thalia democratica within a cold-core eddy off southeast Australia. J. Geophys. Res. 116, https://doi.org/10.1029/2011jc007310 (2011).

53.
Baird, M. E., Everett, J. D. & Suthers, I. M. Analysis of southeast Australian zooplankton observations of 1938–42 using synoptic oceanographic conditions. Deep Sea Res. Part II. 58, 699–711, https://doi.org/10.1016/j.dsr2.2010.06.002 (2011).
ADS  Article  Google Scholar 

54.
Brady, G. S. Report on the Copepoda collected by H.M.S. ‘Challenger’ during the years 1873–76. Rep. Sci. Results, Voy. of H. M. S. “Challenger,” Zool. viii, 1–142 (1883).
Google Scholar 

55.
Knox, G. A. Biology of the Southern Ocean. Second edition (CRC Press/Taylor & Francis Group, 2007).

56.
Scott, A. The Copepoda of the Siboga Expedition, Part I. Free swimming, littoral and semi-parasitic Copepoda Siboga expeditie. 29a, 323 p. (1909).

57.
Farran, G. P. Crustacea, Part X. Copepoda. Natural History Report 8 no. 3. (Trustees Brit. Mus., London, 1929).

58.
Brady, G. S. Copepoda. AAE 1911-1914 Ser C. 5, 1–48 (1918).
Google Scholar 

59.
Russell, F. S. & Colman, J. S. Great Barrier Reef expedition 1928-1929. The zooplankton I. Gear, methods and station list. 2, 2–35 (The Bristish Museum, London, 1931).

60.
Marshall, S. M. The production of microplankton in the Great Barrier Reef region. Sci. Rpts. Gt. Barrier Reef Exped. London 1928–29 2, 111–157 (1933).
Google Scholar 

61.
Jesperson, P. Quantative investigation of the distribution of the macroplankton in different oceanic regions. Rep. “Dana” Exped. 2, 1–44 (1935).
Google Scholar 

62.
Sheard, K. Plankton of the Australian-Antarctic quadrant. Rep. BANZ Antarct. Res. Exped. ser. B 6, 1–120 (1947).
Google Scholar 

63.
Dakin, W. J. & Colefax, A. The marine plankton of the coastal waters of New South Wales. I. The chief planktonic forms and their seasonal distribution. Proc. Linn. Soc. N.S.W. 58, 186–222 (1933).
Google Scholar 

64.
Dakin, W. J. & Colefax, A. The plankton of the Australian coastal waters off New South Wales Part I. Publ. Univ. Syd. 210 (1940).

65.
Thompson, H. Pelagic Tunicates in the plankton of South-eastern Australian Waters, and their place in Oceanographic studies. Bull. CSIR. Aust. 153, 1–51 (1942).
Google Scholar 

66.
Foxton, P. The distribution of the standing crop of zooplankton in the Southern Ocean. Discovery Reports 28, 191–236 (1956).
Google Scholar 

67.
Smith, J. A. et al. A database of marine larval fish assemblages in Australian temperate and subtropical waters. Sci. Data 5, 180207, https://doi.org/10.1038/sdata.2018.207 (2018).
Article  PubMed  PubMed Central  Google Scholar 

68.
Brewer, D. T. et al. Impacts of gold mine waste disposal on a tropical pelagic ecosystem. Mar. Pollut. Bull. 64, 2790–2806, https://doi.org/10.1016/j.marpolbul.2012.09.009 (2012).
CAS  Article  PubMed  Google Scholar 

69.
Dennis, D. et al. Baseline assessment of the status of the tropical marine benthic and pelagic communities adjacent to and distant from the nickel/cobalt refinery site at Basamuk, PNG. Final Report for Ramu NiCo Management Ltd.Commericial-in-confidence (pp. 208. CSIRO Marine and Atmospheric Research, Australia, 2009).
Google Scholar 

70.
Duggan, S., McKinnon, A. D. & Carleton, J. H. Zooplankton in an Australian tropical estuary. Estuar. Coasts 31, 455–467 (2008).
Article  Google Scholar 

71.
McKinnon, A. D., Duggan, S., Holliday, D. & Brinkman, R. Plankton community structure and connectivity in the Kimberley-Browse region of NW Australia. Estuar. Coast Shelf Sci. 153, 156–167, https://doi.org/10.1016/j.ecss.2014.11.006 (2015).
ADS  Article  Google Scholar 

72.
McKinnon, A. D., Duggan, S., Boettger-Schnack, R., Gusmaoa, L. F. M. & O’Leary, R. A. Depth structuring of pelagic copepod biodiversity in waters adjacent to an Eastern Indian Ocean coral reef. J. Nat. Hist. 47, 639–665 (2013).
Article  Google Scholar 

73.
McKinnon, A. D. Planning tools for environmentally sustainable tropical finfish cage culture in Indonesia and northern Australia. Final Report, ACIAR Project FIS 2003/027 (Australian Centre for International Agricultural Research, Canberra, 2009).

74.
Strzelecki, J., Koslow, J. A. & Waite, A. Comparison of mesozooplankton communities from a pair of warm- and cold-core eddies off the coast of Western Australia. Deep Sea Res. Part II. 54, 1103–1112, https://doi.org/10.1016/j.dsr2.2007.02.004 (2007).
ADS  Article  Google Scholar 

75.
McKinnon, A. D., Duggan, S. & De’ath, G. Mesozooplankton dynamics in nearshore waters of the Great Barrier Reef. Estuar. Coast Shelf Sci. Science 63, 497–511, https://doi.org/10.1016/j.ecss.2004.12.011 (2005).
ADS  CAS  Article  Google Scholar 

76.
AIMS. Biological Oceanographic Reconnaissance of the Arafura Sea. https://apps.aims.gov.au/metadata/view/7b5b2d16-052e-054a098-a392-e051a835e64746 (AIMS Data Centre, Australia, 2019).

77.
McKinnon, A. D. & Duggan, S. Summer copepod production in subtropical waters adjacent to Australia’s North West Cape. Mar. Biol. (Berl) 143, 897–907, https://doi.org/10.1007/s00227-003-1153-1 (2003).
Article  Google Scholar 

78.
Taw, N. Studies on the zooplankton and hydrology of south-eastern coastal waters of Tasmania. PhD thesis, Univ. Tasmania (1975).

79.
Terauds, A. In Tasmanian Slope Trophodynamics: Final Report. FRDC Project 91/7 (ed Parslow, J. et al.) 213 (CSIRO Division of Fisheries, Australia, 1995).

80.
Pausina, S. et al. In Moreton Bay Quandamooka & Catchment: Past, present and future. (eds Rothlisberg, P. C. et al.) The Moreton Bay Foundation. Brisbane, Australia., pp. 335–360 (The Moreton Bay Foundation, Australia, 2019).

81.
Lynch, T. P. et al. In Oceans 2008, Vols 1-4 Oceans-Ieee 367–374 (Ieee, 2008).

82.
Armstrong, A. O. et al. Prey Density Threshold and Tidal Influence on Reef Manta Ray Foraging at an Aggregation Site on the Great Barrier Reef. Plos One 11, https://doi.org/10.1371/journal.pone.0153393 (2016).

83.
Furnas, M. J., Mitchell, A. W., Gilmartin, M. & Revelante, N. Phytoplankton biomass and primary production in semi-enclosed reef lagoons of the central Great Barrier Reef, Australia. Coral Reefs 9, 1–10, https://doi.org/10.1007/bf00686716 (1990).
ADS  Article  Google Scholar 

84.
CSIRO. Oceanographical observations in the Indian Ocean in 1959. H.M.A.S. Diamantina cruises Dm1/59 and Dm2/59. Oceanographical Cruise Report 1 (Division of Fisheries and Oceanography, CSIRO Australia, 1962).

85.
CSIRO. Oceanographical observations in the Indian Ocean in 1960. H.M.A.S. Diamantina cruise Dm1/60. Oceanographical Cruise Report 2 (Division of Fisheries and Oceanography, CSIRO Australia, 1962).

86.
CSIRO. Oceanographical observations in the Indian Ocean in 1961. H.M.A.S. Diamantina cruise Dm2/61. Oceanographical Cruise Report 9 (Division of Fisheries and Oceanography, CSIRO Australia, 1963).

87.
CSIRO. Oceanographical observations in the Indian Ocean in 1961. H.M.A.S. Diamantina cruise Dm1/61. Oceanographical Cruise Report 7 (Division of Fisheries and Oceanography, CSIRO Australia 1963).

88.
CSIRO. Oceanographical observations in the Indian Ocean in 1960. H.M.A.S. Diamantina cruise Dm2/60. Oceanographical Cruise Report 3 (Division of Fisheries and Oceanography, CSIRO Australia, 1963).

89.
CSIRO. Oceanographical observations in the Indian Ocean in 1961. H.M.A.S. Diamantina cruise Dm3/61. Oceanographical Cruise Report 11 (Division of Fisheries and Oceanography, CSIRO Australia, 1964).

90.
CSIRO. Oceanographical observations in the Pacific Ocean in 1960. H.M.A.S. Gascoyne cruises G1/60 and G2/60. Oceanographical Cruise Report 5 (Division of Fisheries and Oceanography, CSIRO Australia, 1962).

91.
CSIRO. Oceanographical observations in the Pacific Ocean in 1961. H.M.A.S. Gascoyne cruise G1/61. Oceanographical Cruise Report 8 (Division of Fisheries and Oceanography, CSIRO Australia, 1963).

92.
CSIRO. Oceanographical observations in the Pacific Ocean in 1961. H.M.A.S. Gascoyne cruise G4/61. Oceanographical Cruise Report 12 (Division of Fisheries and Oceanography, CSIRO Australia, 1967).

93.
CSIRO. Oceanographical observations in the Pacific Ocean in 1962. H.M.A.S. Gascoyne cruise G1/62. Oceanographical Cruise Report 13 (Division of Fisheries and Oceanography, CSIRO Australia, 1967).

94.
Crooks, A. D. Coastal hydrological investigations in the New South Wales tuna fishing area, 1960. Oceanographical Station List 53 (Division of Fisheries. CSIRO Australia, 1963).

95.
CSIRO. Investigations by F.R.V. Derwent Hunter on the eastern Australian tuna grounds in 1961. Oceanographical Station List 54 (Division of Fisheries. CSIRO Australia, 1968).

96.
CSIRO. Investigations by F.R.V. Derwent Hunter on the eastern Australian tuna grounds in 1962. Oceanographical Station List 59 (Division of Fisheries, CSIRO Australia, 1968).

97.
CSIRO. Coastal investigations at Port Hacking, New South Wales, 1961. Oceanographical Station List 81 (Division of Fisheries, CSIRO Australia, 1961).

98.
CSIRO. F.R.V. “Derwent Hunter”. Scientific report of Cruises 13/57-16/57. Report of Division of Fisheries and Oceanography 21 (Division of Fisheries and Oceanography, CSIRO Australia, 1959).

99.
CSIRO. F.R.V. Derwent Hunter. Scientific report of Cruises 5/57-8/57. Report of Division of Fisheries and Oceanography 19 (Division of Fisheries and Oceanography, CSIRO Australia, 1958).

100.
CSIRO. F.R.V. Derwent Hunter. Scientific report of Cruises DH 9/57-12/57. Report of Division of Fisheries and Oceanography 20 (Division of Fisheries and Oceanography, CSIRO Australia 1959).

101.
Rothlisberg, P. C. & Jackson, C. J. Temporal & spatial variation of plankton abundance in the Gulf of Carpentaria, Australia 1975–1977. J. Plankton Res. 4, 19–40, https://doi.org/10.1093/plankt/4.1.19 (1982).
Article  Google Scholar 

102.
CSIRO. RV Southern Surveyor dropnet data files: C2012/7164, C2012/7171, C2012/7166. (CSIRO Oceans & Atmosphere DataCentre, Australia, 2012).

103.
Ikeda, T. et al. Biological, chemical and physical observations in inshore waters of the Great Barrier Reef, North Queensland 1975–1978. AIMS Techn. Bull. AIMS-OS-80-1. 56 (AIMS, Australia, 1980).

104.
Tranter, D. J. Zooplankton abundance in Australasian waters. Aust. J. Mar. Freshwater Res. 13, 106–142 (1962).
MathSciNet  Article  Google Scholar 

105.
Young, J. W., Bradford, R. W., Lamb, T. D. & Lyne, V. D. Biomass of zooplankton and micronekton in the southern bluefin tuna fishing grounds off eastern Tasmania, Australia. Mar. Ecol. Prog. Ser. 138, 1–14 (1996).
ADS  Article  Google Scholar 

106.
CSIRO. Australian Equatorial JGOFS Dataset. CSIRO Division of Marine Research, http://www.marine.csiro.au/datacentre/JGOFSweb/inventory/index.htm (2000).

107.
Motoda, S., Kawamura, T. & Taniguchi, A. Differences in productivities between the Great Australian Bight and the Gulf of Carpentaria, in summer. Mar. Biol. (Berl) 46, 93–99, https://doi.org/10.1007/bf00391524 (1978).
CAS  Article  Google Scholar 

108.
Sheard, K. Plankton characteristics at the Cronulla Onshore Station, New South Wales 1943-46. Bull. CSIR. Aust. 246, 1–23 (1949).
Google Scholar 

109.
Thompson, P., Lourey, M., Strzelecki, J., Wild-Allen, K. & McLaughlin, J. In Final report. Southwest Australian Coastal Biogeochemistry (ed J. Keesing) (CSIRO, 2011).

110.
Davis, T. L. O. & Clementson, L. A. Data report on the vertical and horizontal distribution of tuna larvae in the east Indian Ocean, January-February 1987. CSIRO Marine Laboratories Report (CSIRO Australia, 1990).

111.
van Ruth, P. D. Spatial and temporal variation in primary and secondary productivity in the easten Great Australian Bight PhD thesis, Univ. Adelaide (2009).

112.
Stevens, J. D., Hausfeld, H. F. & Davenport, S. R. Observations on the biology, distribution and abundance of Trachurus declivis, Sardinops neopilchardus and Scomber australasicus in the Great Australian Bight. CSIRO Marine Laboratories Report 27 (CSIRO Australia, 1984).

113.
Clementson, L. A., Harris, G. P., Griffiths, F. B. & Rimmer, D. W. Seasonal and inter-annual variability in chemical and biological parameters in Storm Bay, Tasmania. 1. Physics, chemistry and the biomass of components of the food chain. Aust. J. Mar. Freshwater Res. 40, 25–38 (1989).
CAS  Article  Google Scholar 

114.
Swadling, K. M., Eriksen, R. S., Beard, J. M. & Crawford, C. M. Marine currents, nutrients and plankton in the coastal waters of south eastern Tasmania and responses to changing weather patterns. (FRDC Project No 2014/031) 99 (University of Tasmania, Australia, 2017). More

Heterogeneous and affluence-influenced footprints of APAC
The amount of natural resource extractions, environmental emissions, and socio-economic influences associated with satisfying the final demand of average person vary significantly among the eight APAC countries/regions (Fig. 1). The variations correlate with differences in affluence levels (Supplementary Table 2), which is consistent with the findings of previous footprint research focusing on GHG emissions4, and blue water consumption42 across countries worldwide. While the footprints of lower-income countries (e.g., India, Indonesia, RoAP, and China) have also increased with poverty alleviation, yet, most of them are still below the global averages by 2015. In addition, based on the varying environmental footprints (i.e., blue water, energy, and GHG) of the region’s six countries assessed annually from 1995–2015, we find some evidence for the Environmental Kuznets Curve (EKC) hypothesis, i.e., environmental pollution first rises and then falls as economic development proceeds (Supplementary Fig. 1 and Supplementary Fig. 2). Previously, few studies have addressed the socio-economic implications driven by final demand, such as the labor inputs required (employment) and revenues generated (value added)1. Here we find that a country/region’s employment and value-added footprints are both positively correlated with its affluence level.
Fig. 1: The environmental-social-economic footprints of APAC countries/regions in 1995 and 2015.

The six footprint indicators fall into three categories and are presented in a natural resource, b local and global environmental threats, c socio-economic effects. The eight countries/regions are aligned on the y-axis following a descending order by average income per capita in 2015. Footprints are further broken down to the impacts occurred within the country/region, abroad in other APAC countries/regions, and abroad in non-APAC countries/regions. In each subplot, the dashed line indicates the world average level. GHG (CO2, N2O, and CH4) are measured in CO2 equivalent based on the 100-year Global Warming Potential (GWP100). Note, in the plots for value added, we use the insets to show the footprint compositions of the last four countries/region. China refers to Chinese mainland, and we call Taiwan, China as Taiwan for short in the rest.

Full size image

Affluence levels also affected the geo-compositions of the footprint indicators in the APAC region. Richer APAC economies showed high and increasing reliance on outsourcing blue water consumption, PM2.5 emissions, and labor abroad, especially within the APAC region. Specifically, 58, 56, and 52% of the blue water footprints of Japan, South Korea and Taiwan were traced to water consumed in other APAC economies in 2015 (Fig. 1), primarily through the intra-APAC trade of agricultural products, such as sugar cane/sugar beet and paddy rice (Supplementary Fig. 3). In contrast, the final demand of middle- to low-income APAC economies has been largely satisfied by domestic natural and labor resources over the two decades. Moreover, during the same period, the PM2.5 footprint of high-income APAC economies (Australia, Japan, South Korea, and Taiwan) that occurred in other APAC economies increased from ~17 to 40%. Yet, for the less wealthy countries (e.g., China, India, and Indonesia), 75−99% of PM2.5 footprints were indigenous, a considerable fraction of which are attributed to direct household emissions (Supplementary Fig. 3). The consumption of traditional fossil fuels and biomass (e.g., fuelwood and agriculture wastes) for home cooking and heating, and the open-air combustion of biomass, especially in rural areas of India and China has been considered a significant source of PM2.5 emissions43,44,45.
Among all the indicators, we find the value-added footprints of the APAC economies, i.e., the contributions of their final demand to global economic growth (domestic + intra-APAC + non-APAC) experienced the most significant increases. Such increases are especially significant for the middle- and low-income APAC economies. For example, China’s per capita value-added footprint increased by eight times during 1995−2015, followed by India (217%), Indonesia (113%), and RoAP (82%). The rates are much higher than those experienced by the high-income economies (28% on average).
The APAC economies have grown more interdependently linked during the past two decades. Depending on the footprint indicators, the foreign footprint shares traced to APAC countries (Abroad, APAC in Fig. 1) increased by 4–27% while the domestic shares decreased by 3–33%. Given that the compositional changes within the aggregated RoAP cannot be estimated, the variation of RoAP is not included in the ranges. The increased intra-APAC interdependencies are primarily attributable to the strengthened linkages among the six major APAC countries. In contrast, India’s environmental-social-economic footprints remained predominantly domestic (90–99% in 1995 and 82–97% in 2015). Previously, researchers highlighted that resource constraints have already become a bottleneck for India’s social and economic development, such as fresh water scarcity46 and energy deficiency47. Air quality deterioration caused by PM2.5 emissions has also made India under severe health burden48. Our results here confirm that India remains highly dependent on local resource use and labor-intensive production activities to sustain its socio-economic growth. Engaging in international trades has the potential of reducing the depletion of local scare resources (e.g., blue water and energy) through import while adding employment and added-value through exporting goods produced with abundant labor resources and low environmental impacts locally1,30. The trade-environment relationship is primarily rooted in the economic principle of competitive advantages among countries for international trade. Therefore, adopting the strategy of reducing trade barriers rationally, increasing the openness to the outside world to actively promote international economic cooperation may be one solution to alleviating the pressure of domestic resource depletion (i.e., water and energy) and environmental damages (i.e., carbon and air pollutions) during India’s economic growth.
Footprint outsourcing and disparity through intratrade
The linkages and imbalances among the APAC countries/regions, through the virtual flows of environmental-social-economic resources embedded in intra-APAC trade, are highlighted in Fig. 2 (2015) and Supplementary Fig. 4 (1995). The intra-APAC trade patterns observed over the past two decades confirm that developed economies (Japan, Australia, South Korea, and Taiwan) import natural resources and labor from less-developed regions where resources and labors are cheaper (China, India, and RoAP). In 2015, for the 37 billion m3 of net bilateral virtual water trade in the region, nearly 80% was associated with the exports from RoAP and India, while ~50% was driven by the final demand of Japan and South Korea, mainly embedded in a range of water-intensive agricultural crops and products. Yet, as mentioned above, India has already been threatened by serious water crises with low availabilities of safe drinking water49. For the 15 EJ net bilateral energy flow embedded in the 2015 intra-APAC trade, China was the main net supplier to the rest of the region, contributing 43%, followed by South Korea (27%), which possesses strong petrochemical and steel industries. China’s net energy outflows were predominately embedded in energy-intensive products (e.g., coal power, steel, and gasoline), 73% of which ultimately served RoAP’s final demand. Of the 488 Tg net bilateral flows of GHG emissions, 76% was supplied by China, while the final demand of RoAP contributed more than half of it (286 Tg), and the rest is attributed to the final demand of four higher-income economies (Australia, Japan, South Korea, and Taiwan, 123 Tg in total). 81% of the 677 thousand tons of net bilateral PM2.5 flows occurred in China and India. And 38% of the net flows were driven by the final demand of four higher-income economies. The social and economic effects of the intra-APAC trade are more nuanced than those related to natural resources and emissions. The Intra-APAC trade resulted in a net bilateral outpouring of 99 million people, i.e., labor resources, in 2015. RoAP was the largest net labor supplier, contributing 77%, while Australia, Japan, South Korea, and Taiwan had net inflows of labor to satisfy their final demands, equivalent to employing 11.4, 27.6, 9.3, and 5.6 million people from the rest of APAC, respectively. In terms of net flows of value added, the intraregion pattern appears significantly different and almost opposite from those of other indicators. By net outpouring value added to other APAC economies, China turned from the second trade surplus country in 1995 to the largest one within APAC in 2015, followed by the developed economies. Yet, on the per capita basis, South Korea, Australia, Japan and Taiwan achieved the most prominent economic gains through the intra-APAC trade.
Fig. 2: Net environmental-social-economic virtual flows of the intra-APAC trade in 2015.

The top five net flows are shown for each footprint indicator, and fall into three panels a natural resource, b local and global environmental threats, c socio-economic effects. The width of the arrows in each panel represents the magnitude of the net flow within the APAC region. The background colors indicate the specific net footprint (import-export) per capita of each region/country. The negative net footprint indicates net displacements (of resource use, emissions, labors, and economic value added) to other APAC countries/regions. Note, the background color and flows of Taiwan, China could not be shown on the map.

Full size image

Our results further demonstrate that the economic and environmental inequity owing to the intraregional trade worsened along time. From 1995 to 2015, the majority of the environmental externalities (i.e., more than 90% of the virtual water flows) were shifted from higher-income to lower-income countries, with a considerable worsening trend for PM2.5 (Supplementary Table 3). As for the economic gains associated with the intraregional trade, higher-income countries’ share increased from 38 to 59% from 1995 to 2015. At country level, China experienced the most significant transformation in intra-APAC trade, especially in the virtual trade of water and labor. For virtual blue water flows, China turned from the second largest net exporter of the region in 1995 to the third largest net importer in APAC (5.3 billion m3), largely due to the reverse of the import-export relationship with RoAP. Over the same period, the net virtual water export of India surged by 161%, despite the aggravating water stress concerns. China turned from the largest net supplier of employment (Supplementary Fig. 4) to the second net demander in APAC after two decades, while RoAP showed an opposite trend. For example, the labor- and resource-intensive textile and apparel production industries have been shifted from China to other less-developed countries, such as Vietnam and Bangladesh50. Such a role switching stemmed from the rising labor costs in China, which is expected to drive more low-end manufacturers to low-cost foreign economies in the coming years51,52. Overall, based on the consistently calculated virtual flows of multiple indicators, we highlight a growing environmental-socio-economic disparity within the APAC region, owing to the intraregional trade. The undesirable environmental externalities are primarily and increasingly shifted from higher to lower-income economies, while higher-income economies achieved more economic gains.
The role of APAC and intra-APAC trade in globalization
APAC is becoming an important player in the globalization process: as natural resource suppliers and manufacturers for the rest of the world, for managing global environmental emissions, and in the labor and monetary markets (Fig. 3). By 2015, APAC-related share of these categories had surpassed 50% (the red, yellow, and blue parts in Fig. 3). In contrast, the international trade without APAC countries (gray in Fig. 3) shrank for all the indicators. Moreover, for all the footprint indicators, the intra-APAC trade and non-APAC’s outsourcing to APAC (red and yellow parts in Fig. 3, respectively) account for a considerable and increasing fraction of the global trade. The former grew from 17 to 20% on average and the latter from 23 to 27% on average. APAC’s outsourcing to non-APAC countries (blue in Fig. 3) only grew slowly, with an average of 14 and 16% in 1995 and 2015, respectively. Earlier studies have proven that the world’s dominating embedded labor flows originate from developing countries, predominately to satisfy the final demand of wealthier economies23. Here, we further elucidate that non-APAC economies became more dependent on offshoring the resources, emissions and labor-intensive industries to APAC over the investigation period.
Fig. 3: Natural resources, environmental emissions, and socio-economic factors embedded in exports.

a natural resource, b local and global environmental threats, c socio-economic effects are shown in three panels. Note, the red and blue sections correspond to the same components in Fig. 1.

Full size image

Among the six indicators, we find APAC’s share in global value added was the smallest. In particular, APAC’s economic gains from exports (red + yellow) were significantly smaller (27–32%) than APAC’s resources, emissions, and labor embedded in its exports (40–48%). Intuitively, this may be due to the fact that the APAC region is dominated by population from developing countries, thus resource/labor-intensive and low value-added products dominate the region’s export. In comparison, non-APAC countries are relatively more skilled in producing products and with higher value added and lower intensities of resource and emissions4. Also, APAC’s relatively small share in the value-added dimensions of the global supply chains can be a result of the overall low resource efficiency of the region12. As resource extractions and environment emissions become increasingly outsourced to the region, where environmental regulations and efficiency measures are only emerging, more environmental impacts will likely be resulted on the global level, offsetting or even reversing the resource efficiency gains and climate change mitigation efforts achieved in developed countries.
On the positive note, from 1995 to 2015, APAC’s resource and environmental intensities declined substantially, both from the perspective of footprint, i.e., footprint per final expenditure, and from the perspective of trade, i.e., direct impacts/gross trade (see Materials and Methods). Such improvement is most noticeable for labor (Table 1). More specifically, over the past two decades, APAC improved faster than the world averages in blue water, PM2.5, and labor requirements per expenditure of final consumer, which decreased by 30, 39, and 35%, respectively, while lagging behind the global averages in reduction pace of energy and GHG intensity from footprint perspectives. The intra-APAC trade achieved a reduction in the intensities of energy, GHG and labor by 22, 42, and 34%, respectively. Relatively, this was higher than the world averages over the period. However, the PM2.5 embedded in intra-APAC trade per traded values even grew by 7%, whereas the same indicator decreased by 34% at the global scale. This can be explained primarily by the PM2.5 trade magnitude in APAC nearly tripling in 2015, much higher than that of the global average (0.8 time). This indicates that the PM2.5 issue has become graver for the APAC trade after two decades. The ratio of traded value added to total intratrade values measures the economic gains of trade. Value added is a general indicator of economic performance, yet, it could not evaluate all the advantages of trade, which may cause misleading results in some cases. Therefore, more comprehensive indicators need to be incorporated to assess its advantages.
Table 1 Intensity comparison between APAC and the global average from the footprint and trade perspectives.
Full size table

Despite APAC’s improvement, APAC’s resource and environmental intensities remain higher in general than the world averages in 2015. The comparison of intensity changes at the national and regional levels was also shown in Supplementary Fig. 5. This can be attributed to the characteristics of the APAC region and may not be altered soon. Specifically, the socio-economic development is supported by resource- and emission-intensive productions in primary industries, especially those linked with capital development and those to satisfy the export demand. The APAC region has a long way to go to achieve a greener and more sustainable trade pathway.
Our result further reveals that, China played a crucial role in APAC’s transformation to greener trade and greener consumption. Generally, without China’s improvement, all the footprint intensities of APAC would exhibit a further increase from 1995–2015, at 8% on average. By contrast, the value is −20% with China in the picture. For the trade intensities, the decline ranges would be much smaller (−1% on average) without China as opposed to with China (−24% on average). Another noteworthy finding is that when eliminating China from APAC, the intensities of footprint and trade in 1995 would be 20–30% lower, implying that China had turned from the lagger to the leader of efficiency improvement in the APAC region. More

Study area
This study was conducted at the Sierra Morena mountain range in Southwestern Spain (38° 22′ 50.64″ N, 4° 45′ 27.69″ W; Pozoblanco, Córdoba). This area has a continental-Mediterranean climate characterised by cold winters and hot, dry summers. Annual precipitation is 439 mm year−1 and mean annual temperature is 15.2 °C. The studied ecosystem is a “Mediterranean dehesa”, a savannah-type forest managed as a traditional sylvo-pastoral system, with an herbaceous stratum of native pasture and a tree layer of scattered oak trees (Quercus ilex L.). Tree density is 14.5 ± 1.3 trees ha−1, and the herbaceous layer is dominated by a diverse community of annual native species such as: Hordeum murinum (L.), Senecio vulgaris (L.), Bromus madritensis (L.), and Sinapis alba (L.), with a mean species richness of 9.6 ± 0.3 species m−2.
Experimental design
In September 2016, we selected three sites separated by 4.3 ± 1.9 km with similar tree density, topography, slope, orientation, and soil type. At each site, we sampled two areas: a partially shaded site under Q. ilex canopy (under tree canopy) and a nearby open grassland site (open grassland). These two habitat types are characteristic of the dehesa ecosystem in the region. Then, within each area (i.e. habitat type) at each site we established six permanent 4 × 6 m (1.2 m high) fenced plots (N = 36 plots). Each plot was divided into two subplots (N = 72 subplots), one subplot was subjected to a rainfall reduction treatment whereas the other was exposed to ambient rainfall. Rainfall was reduced by placing a 2.5 × 2.5 m rain-exclusion shelter over one of the subplots (Fig. 2), resulting in a 30% reduction in annual rainfall22 as predicted by the IPCC for Mediterranean regions23. In addition, within each subplot we sampled two ca. 1-m2 quadrats (144 quadrats in total, i.e. 2 quadrats × 72 subplots): one was subjected to a warming treatment and the other was not manipulated (i.e. ambient temperature). Warming was achieved by placing a 40 × 50 × 32 cm (ca. 0.65 m2) hexagonal open top chamber of methacrylate material without UV filter (Faberplast, Madrid) at the centre of the quadrat (Fig. 2), increasing air temperature by ca. 2 °C compared to ambient temperature (17.5 ± 0.1 vs. 15.7 ± 0.1 °C, respectively). This matches the temperature increase for the study region by climatic forecasting models (SRES A-2 model by the Intergovernmental Panel on Climate Change23). Ambient air temperature was slightly lower (approx. 0.5 °C) in the shaded habitat than in the open habitat, but the experimental temperature increase was of similar magnitude in both cases. Temperature in each quadrat was measured with data loggers. Based on the above set-up, the experiment followed a randomized split–split plot design replicated across three sites, with habitat type as the whole factor (plots as replicates), rainfall manipulation as the split factor (subplots), and temperature manipulation as the split–split factor (quadrats).
Figure 2

Experimental design. Picture showing open top chambers and rain-exclusion shelters used for our climate manipulations. Photo credit: Ignacio M. Pérez-Ramos.

Full size image

Sampling, measurements and chemical analyses
In April 2017 and April 2018, i.e. 1 and 2 years after establishment of climatic treatments, we identified all plant species and estimated their frequency within each quadrat. Then, in April 2018, after the second survey of species frequencies, we collected 3–4 fully expanded leaves from each plant for chemical analyses which were pooled to obtain a single sample per species and quadrat. Plant sampling was restricted to a central portion of 0.65 m2 in each quadrat (within chamber in the case of warming quadrats) as a function of the surface covered by and sampled within the hexagonal chamber placed in quadrats subjected to warming (see above). Due to lack of plant material, we only sampled 127 out of the 144 quadrats. In total, we collected 436 samples from 28 annual plant species (Table S1). Immediately after collection, we oven-dried leaves for 48 h at 40 °C, ground them with liquid nitrogen, and stored them at room temperature before conducting chemical analyses.
We determined total phenolic content colorimetrically by the Folin-Ciocalteu method24,25. Briefly, we extracted phenolics from 20 mg of plant tissue with 70% methanol in an ultrasonic bath for 15 min, followed by centrifugation. We determined total phenolic content colorimetrically by the Folin–Ciocalteu method in a Biorad 650 microplate reader (Bio-Rad Laboratories, PA, USA) at 740 nm, using tannic acid as standard. We performed three technical replicates of each sample in order to estimate variations due to the experimental procedure, and expressed concentrations based on dry weights (d.w.). We decided to use total phenolics as these could be measured across all species in order to address the goal of testing for community-level effects of climatic stressors on plant secondary chemistry.
Statistical analyses
For each quadrat, we calculated the CWMphenolics using the weighted.mean function in R software26. For this, we used phenolic values at the species level weighed by the abundance (i.e. frequency) of each species. We then performed linear mixed models using PROC MIXED in SAS 9.4 (SAS Institute, Cary, NC)27 to test for the effects of habitat type, rainfall manipulation, temperature manipulation, and their two-way and three-way interactions (all fixed factors) on CWMphenolics. We included site, the site × habitat type interaction, and the site × habitat type × rainfall manipulation interaction as random factors. The two- and three-way interaction interactions tested the main effect of habitat type and rainfall manipulation (respectively) with the appropriate error terms according to a split-split plot design. We log-transformed original values to achieve normality of residuals.
We also assessed whether climatic manipulations influenced the community-level expression of phenolic compounds through changes in plant composition with a PERMANOVA (‘vegan’ package in R) that included the effects of climatic treatment, year, and their interaction. More

Origin of hydrate-bound hydrocarbons
A relationship between C1/(C2 + C3) and C1 δ13C has been applied to identify the sources of hydrocarbons in submarine seeps24. Recently, this diagram was revised based on a large dataset25. As shown in Fig. 4a, hydrate-bound hydrocarbons at the Kedr MV have thermogenic and/or secondary microbial origins, whereas those of other gas hydrate sites (Malenky, Bolshoy, Malyutka, Peschanka P-2, Kukuy K-0, Kukuy K-2 and Goloustnoe; Fig. 1) in Lake Baikal demonstrate microbial or early mature thermogenic origins. The hydrate-bound C1 from all locations except those at the Kedr MV were interpreted to be of microbial origin via methyl-type fermentation23 according to Whiticar’s old diagram26; however, the revised diagram25 suggests early mature thermogenic gases (Fig. 4b). Those of the Kedr MV plot at the boundary of the thermogenic and secondary microbial origin zones. Low C1 and C2 δ13C at the Peschanka P-2 MV indicated that C1 and C2 are of microbial origin27,28, whereas Kedr MV shows high C1 and C2 δ13C indicating their thermogenic origin (Fig. 4c). At other sites, C1 and C2 δ13C suggested that gases are mainly of microbial origin (in terms of C1) with some thermogenic component (13C rich and higher concentration in C2).
Stable isotopes in hydrate-bound C1 at the Kedr-1 and Kedr-2 areas suggested its thermogenic origin. However, it is close to the field of secondary microbial C1 in Fig. 4b, and the data are plotted in the overlap between the fields of thermogenic and secondary microbial in Fig. 4a. Milkov29 mentioned that secondary microbial C1 is characterised by C1-rich dry gas, large C1 δ13C (between − 55‰ and − 35‰) and large CO2 δ13C (more than + 2 ‰). Although hydrate-bound and sediment gases in the Kedr MV were not C1 rich and contained 3%–15% of C2, C1 δ13C was around − 45‰, which agrees with the secondary microbial C1. Because some data of secondary microbial gas are plotted outside the field on the original graph25, we could include the gas data in the category of secondary microbial C1 in Fig. 4b.
Figure 6 shows the relationship between C1 δ13C and CO2 δ13C in the sediment gas obtained using headspace gas method. According to the genetic diagram25, gas hydrate cores are plotted at the zones of the thermogenic and secondary microbial origins, whereas the cores at the peripheral area are primary microbial. The headspace gas data of the hydrate-bearing cores in Fig. 6 seem to be plotted in the field of thermogenic gas (low CO2 δ13C), but the effect of light CO2 produced by methane oxidation in the subsurface layer also decreased CO2 δ13C as shown in Fig. 5. These results suggested that secondary microbial C1 mixes into thermogenic gas. Coal-bearing sediments exist around the Kedr area21,22, and secondary microbial C1 can also form from coal beds30. Hydrate-bound C1 of secondary microbial origin has been only reported at the Alaska North Slope31. This study is another case for it.
Figure 6

A diagram of headspace gases. CO2 δ13C plotted against C1 δ13C, based on the classification of Milkov and Etiope25.

Full size image

Formation process of the sII gas hydrates
As stated before, the crystallographic structure of gas hydrates at the Kedr MV is mainly due to the composition of thermogenic C2 in the volatile hydrocarbons. The concentration of C3, which is one of the sII-forming components, was two to three orders of magnitude smaller than that of C2, because biodegradation occurs and this preferentially reduces C3−5 of n-alkanes19,32, 33. The concentration of n-C4 was smaller than that of i-C4, whereas that of n-C5 was not detected (Table 1). C3 δ13C was around − 10‰, suggesting that light C3 is consumed by microbial activity. Assuming that sediment gas C3+ can be ignored, sediment gas ratio C1/C2 at the study area was 30 ± 17 (mean and standard deviation), and the concentration of C2 was ~ 3%. Such a composition of thermogenic gas is, therefore, considered to be supplied from a deep sediment layer, forming sI gas hydrates composed of mainly C1 and C211,12 in the lake floor sediment.
In the cases where sI gas hydrates plug and block migration pathways, upward fluid flow becomes more focused in other areas16. Once gas supply stops locally, gas hydrates begin to decompose, with the gas dissolving into gas–poor sediment pore water. In the system of C1 and C2, C2 is prone to be encaged in gas hydrate and decreases the equilibrium pressure of mixed-gas hydrate. Therefore, C2-rich gas hydrate forms in parallel with the decomposition of sI gas hydrate. The Colorado School of Mines Hydrate (CSMHYD) program34 showed that C2-rich sII gas hydrate (C2 concentration 17%) forms from mixed gas composed of C1 and C2 (C2 concentration 3%). The C2 concentration of hydrate-bound gas at the Kedr MV was ~ 14%, agreeing fairly well with the results of the CSMHYD program. Such secondary generation of gas hydrates can produce compositions and crystallographic structures that are different from the original crystals. A calorimetric study of synthetic C1 and C2 mixed-gas hydrate revealed that double peaks of heat flow correspond to the dissociation process of C1 and C2 mixed-gas hydrate, suggesting that C2-rich gas hydrate forms simultaneously from dissociated gas and showed that the second heat flow peak correspond to the dissociation of C2-rich gas hydrate18. The PXRD and solid-state 13C nuclear magnetic resonance techniques demonstrated that C2-rich sI gas hydrate forms in the dissociation process of C1 + C2 sII gas hydrate35.
Among twenty hydrate-bound cores in the Kedr area, four cores contained sI only, seven cores had sII only, and seven cores showed sII at the upper layer and sI at the lower layer, as observed at the Kukuy K-2 MV13,16,17. Furthermore, in the cores 2015St1GC15 and 2016St18GC2, gas hydrate structure had sI at the upper and lower layer, and sII at the middle layer. These results suggested that complex gas hydrate layers are composed of sI and sII in subsurface sediments as shown in the schematic illustration in Poort et al.16.
Depth profiles of C2 δ2H of gas hydrate cores from the Kedr MV are shown in Fig. 7. C2 δ2H of hydrate-bound gases varied between − 227‰ and − 206‰, with a grouping around − 210‰. C2 δ2H of sediment gases was also around − 210‰, indicating that C2 δ2H of the original thermogenic gas is − 210‰. As stated above, C2 δ2H of some cores showed low values at their base. Based on the isotopic fractionation of hydrogen in C2 during the formation of sI C2 hydrate36, δ2H of hydrate-bound C2 was 1.1‰ lower than that of residual C2. However, this is too small to explain the wide distribution in C2 δ2H shown in Fig. 7. On the other hand, Matsuda et al.37 reported that isotopic fractionation of hydrogen in C2 is dependent on the crystallographic structure: 1‰–2‰ for sI and ~ 10‰ for sII. Gas hydrates plotting around − 220‰ in C2 δ2H can be explained as a secondary generation of sII from dissociated gas hydrates, of which C2 δ2H was around − 210‰. However, some sII samples showed high C2 δ2H (around − 210‰), whereas some sI samples showed low C2 δ2H (around − 220‰). These results indicated that formation and dissociation processes of gas hydrates produce complicated isotopic profiles in C2 δ2H under non-equilibrium conditions.
Figure 7

Depth profiles of C2 δ2H of hydrate-bound and sediment gases. cmblf, centimetres below lake floor.

Full size image

Characteristics of hydrate-bound gases in sII
C3, i-C4, n-C4 and neo-C5 can be encaged in the larger hexadecahedral cages of sII1. n-C4 and neo-C5 can be encaged using a help gas (e.g. C1) to fill in the smaller dodecahedral cages of sII, because they cannot form pure n-C4 and neo-C5 hydrates, respectively. Figure 8 shows the concentration of C3, i-C4, n-C4, neo-C5 and i-C5 plotted against C2 concentration. The figure illustrates a clear division between sI (3–4%) and sII (14%) C2 concentrations. Data points between C2 concentrations of 5% and 13% were considered to have a mixture of sI and sII. Concentrations of C3, i-C4, n-C4 and neo-C5 had a positive correlation with the concentration of C2, and these concentrations in sII were 1 or 2 orders of magnitude larger than those in sI, suggesting that C3, i-C4, n-C4 and neo-C5 are encaged with C2 in the sII formation process.
Figure 8

Concentration of C3–5 against C2 concentration in the hydrate-bound gases.

Full size image

C3 values of 0.001%–0.01%, ~ 0.0001% of n-C4, and 0.0001%–0.01% of neo-C5 were also detected in sI hydrate-bound gas (Fig. 8), despite these hydrocarbons being unable to be encaged in sI. This can be explained by gases being adsorbed with sediments and gas hydrate crystals, which are then trapped in the grain boundary of polycrystalline gas hydrate crystals, and the gases are encaged if a small amount of sII crystals are present. For example, Uchida et al.38 examined natural gas hydrate retrieved at the Mackenzie Delta (onshore Canada) and detected C3 encaged in sII using Raman spectroscopy, although PXRD results suggested that the sample was sI and the major component of hydrate-bound gas was C1 (more than 99%).
neo-C5 is considered to form from the decomposition of gem-dimethylcycloalkanes derived from the terpenes of terrestrial organic matter39. It is easily enriched by preferential diffusion due to the nearly spherical molecules and its diffusion coefficient, which is higher than that of less branched isomers40. The sII hydrates retrieved at the Kukuy K-2 MV (central Baikal basin) contained 0.026–0.064% of neo-C5 in the volatile hydrocarbons13,14, and those at the Kedr MV had a maximum value of 0.054% of neo-C5 (Supplementary Information Table S1). On the contrary, in the case of natural gas hydrates retrieved at the Joetsu Basin (Japan Sea), neo-C5 was excluded and remained in sediment during the formation of sI gas hydrates from C1-rich gas41. The molecular size of i-C5 is considerably large to be encaged in the large cages of sII. Maximum concentration of i-C5 in the hydrate-bound gases was in several parts per million in both the fields of sI and sII (Fig. 8), indicating that i-C5 is not a hydrate-bound hydrocarbon and adsorbed with gas hydrate crystals and/or trapped in their grain boundary. More

General scheme
We reviewed the literature published until April 2018 on the economic impacts of invasive species. For reasons of feasibility (linguistic skills of the review team, restriction to a reasonable scale of the review), we conducted all searches in the English language assuming that a large body of knowledge (mostly from international peer-reviewed papers and reports) is written in English. The dates of each search process were systematically recorded. We used the following strategy for all repositories (Fig. 1), while also taking into consideration the specificity of their algorithms.
First, a literature search was performed using three online bibliographic sources successively to minimize the risk of omitting relevant materials (Fig. 1, step 1a): ISI Web of Science platform (https://webofknowledge.com/), Google Scholar database (https://scholar.google.com/) and the Google search engine (https://www.google.com/). We carefully composed appropriate search strings that were consensually retained as the most efficient among a set of potential candidates. A decision was taken following preliminary tests based on a handful of relevant articles provided by consulted subject experts on some taxonomic groups (amphibians, reptiles, fishes and ants). Final selection of search strings comprised those considered to have the largest potential to identify key references. Each search string was set to include a combination of two search terms, related to ‘invasive’ and ‘economics’. For both terms, we used a range of synonyms or related words. For example, for ‘invasive’ we used invasi*, invader or exotic; for ‘economics’, we used econom*, cost or monetary. In addition, the search string included exclusion terms to omit mismatches, for example, with studies from the field of medicine that are focused on pathologies or procedures that can be ‘invasive’ for patients. We complemented this search with documents gathered opportunistically (Fig. 1, step 1b). The potentially relevant materials derived from all these sources were combined in a single file and screened for duplicates. Second, retrieved documents were individually assessed at progressive levels (titles, then abstracts, keywords, and finally full text when abstracts were missing; Fig. 1, step 2) based on three criteria. Hence, materials were deemed relevant if (i) they matched with the linguistic competencies of the review team (i.e. written in English, or French where English language was restricted only to the title and/or abstract) for allowing reliable assessment, (ii) they contained at least one cost estimate (studies exclusively providing benefit estimates from direct use or exploitation of invasive species were excluded), and (iii) that this cost estimate is exclusively associated with invasive species (estimates merging non-invasive and invasive species, without the possibility of distinguishing the respective contribution of each group to the overall cost, were excluded). To ensure transparency and validity, each document was checked by two reviewers and in case of a disagreement between assessors, a third reviewer was involved. However, it was often difficult to judge from the topic whether the content of an article was relevant and so consequently many more articles were conservatively kept when final agreement was lacking among assessors.
Finally, relevant materials were scrutinized for data on economic costs (Table 1; Fig. 1, step 3). During this step, additional relevant materials were found as cited by the analysed materials. Obtained cost data were collated in a database and the costs were converted to a common and up to date currency (2017 US$), and then depicted by different descriptors. Categories extracted from relevant materials allow search of the database and data pre-selection to facilitating analysis of costs based on taxonomic groups, geographical areas, impacted sectors, types of costs, or other categories. The reliability of cost estimates and all associated information recorded in the final InvaCost database was systematically checked at least twice, and every ambiguous element was discussed to reach a consensus. We also checked all entries in the database to ensure that there were no obvious duplicate reports (i.e. multiple documents reporting the same cost estimate) or mistakes.
Hereafter, we specifically describe each of the steps made to generate InvaCost.
Literature search
Web of Science
We used the Web of Science (hereafter called WoS) to conduct a search for potentially relevant materials on 7 December 2017 (Fig. 1, step 1a). We applied the following search string: (econom* OR cost OR monetary OR dollar OR euro OR “sterling pound”) AND (invasi* OR alien OR non-indigenous OR nonindigenous OR nonnative OR non-native OR exotic OR introduced OR naturali* OR invader) NOT (cancer* OR cardio* OR surg* OR carcin* OR engineer* OR rotation OR ovar* OR polynom* OR purif* OR respirat* OR “invasive technique” OR carbon OR fuel OR therap* OR vehicle OR cell* OR drug OR fitness OR “operational research” OR banking OR liberalization). The terms were searched in the field code “Topic” which includes title, abstract and keywords, and which also comprises ‘Keywords Plus’ that are generated by WoS through an automatic computer algorithm, based on words and phrases that appear frequently in the titles of article’s bibliographic references and not necessarily in the main text of the article itself. To limit the search to relevant fields of research, we used the function ‘refine’ to exclude subject areas not related to economics and/or invasion biology.
We exported all records (n = 16,875) into an Excel worksheet30 (Table 1) to identify the relevant materials by a two-step procedure. First, we excluded the references identified only based on ‘Keywords Plus’, which were shown to be poor specific descriptors of the content of articles31. We also excluded references identified based on the presence of only a single search term in the topic, as we assumed that words related to both search terms (‘invasive’ and ‘economics’) should be mentioned at least once in the title, abstract and/or keywords of a relevant material. To identify these irrelevant materials within the references collected, we developed a script (see Code Availability) in the R programming language (R v.3.4.3)32. Subsequently, 10,592 references were kept for the next screening step based on the described criteria.
In the second step, the topic of every reference selected was checked manually to ensure potential relevance of its contents. This allowed the elimination of documents incorrectly identified as relevant, such as studies without a true monetary assessment, or those focusing on economic estimates not directly attributable to invasive species only. Finally, 1,333 documents were judged as relevant materials (Table 1) and moved to the final data collation step.
Google scholar
The Google Scholar database is a large source of grey as well as peer-reviewed literature. Nevertheless, we had to modify our approach in order to address inherent limitations of this database as a search tool (see Haddaway et al.33 for a comprehensive analysis). Typically, Google Scholar allows limited Boolean operators (no nesting using parentheses permitted) and search strings are limited to 256 characters. Additionally, only the first 1,000 search results can be viewed and the order in which results are returned is not disclosed. We also wanted to maximize novel information by avoiding too much overlap between the references collected with WoS and those gathered here.
In light of the above, we adapted our search string to generate the most efficient outcome, i.e. sufficiently pertinent to bring the most relevant items to the top of the result list while not unnecessarily large so as to limit the host of non-viewable results. Thus, the following search string was applied on 26 April 2018, using the advanced search facility to search for selected words anywhere in the article (see https://scholar.google.se/intl/en/scholar/help.html#searching for further details): dollars OR euros OR “USD” OR “EUR” OR “NZD” OR “AUD” OR “CAD” OR “GBP” OR “economic cost” OR “economic impact” OR “estimated cost” invasive species. We specified currencies for prioritising materials with monetary data in the top of the resulting list. These currencies were chosen as they were the most often used to express economic costs in the literature collected from the WoS. Nevertheless, any reference evoking economic costs in other currencies was expected to be also captured by some specific combinations of ‘economic’ terms in our search string that we would expect to be mentioned at least once in the full-text of relevant papers. In addition, we included the concomitant presence of ‘invasive’ and ‘species’ terms to restrict the outcomes to papers within the scope of our synthesis. Subsequently, we collected all viewable results (100 pages, n = 992 references of the 668,000 generated), thus going beyond the traditional and arbitrary sample size of first 50–100 results, which is frequently selected in many systematic reviews. We used a web-scraping programme (https://www.webscraper.io/) to extract all the titles’ references returned by the search in an Excel spreadsheet. Because we could not efficiently export the abstract for every reference, we screened them online to assess their potential relevance.
As a result of a search and relevance assessment within Google Scholar, the references, abstracts and specific bibliographic details of 432 documents were added to the sample for further analysis. After excluding duplicates with WoS retrieved references, 310 additional documents were included in the sample as potentially relevant materials (Table 1).
Google
We used the Google search engine to complete the standardised literature search. As when searching with Google Scholar, we took into account specific constraints related to the use of this search engine. Moreover, browsing through Google search results can be overwhelming due to the vast amount of information of highly variable quality. We attempted to implement a search strategy that could allow overcoming these limitations as much as possible. We used the following search string: economic species invasive OR nonnative OR alien OR exotic OR nonindigenous -disease -surgery -fungus -respiratory. We added four exclusion terms (disease, surgery, fungus, respiratory) identified during preliminary tests to restrict the number of irrelevant studies, associated with medical research. We did not use a range of economics-related terms, such as impact or cost, as they returned overly large numbers of mismatches.
The web search was conducted on 8 May 2018 by searching for specified terms within page titles of each document, in order to maximize the likelihood of identifying grey literature. We especially targeted grey literature because searches by the other two platforms mainly led to peer-reviewed publications. We assumed that documents published online by various governmental and non-governmental organisations (NGO), research centres and academic institutes are more likely to contain relevant data than other types of documents such as blogs and catalogues29. Therefore, we restricted our search to the documents located on governmental, academic and NGO webpages to ensure that explicit, traceable and expertise-based information was retrieved. We conducted independent searches for each type of webpage by specifying the type of web extension in the advanced search facility (.gov for governmental,.edu for academic, and.org for organisational webpages).
361 search hits were collected (document name, publishing year and URL of the main website homepage, if available) and stored in the database with the same host of dedicated information (Table 1). If the item analysed was a website homepage, we conducted on-line searches of potentially relevant materials within the website database(s), by filters if available, or by using the search bar with combinations of keywords. Websites that did not contain a database or search bar were searched manually. We then eliminated all duplicates resulting from references being listed on multiple websites, or due to typographical mistakes and/or incomplete records when reporting a reference within different repositories. A total of 119 potentially relevant materials was finally obtained (Table 1).
Targeted collection
Finally, we sourced other potentially relevant materials that did not originate from the above-described processes (Fig. 1, step 1b). On one side, we dedicated specific efforts on gathering cost estimates for particular taxa or areas for which data previously obtained seemed scarce. First, we made sure that some key species were adequately covered; for example, costs associated with invasive mosquito species responsible for much of the burden of mosquito-borne viral diseases worldwide (Aedes aegypti that mainly invaded the intertropical zone from the 15th-17th centuries, and Aedes albopictus for which the global dissemination was more recent34) were searched in a specific way using WoS and PubMed (https://www.ncbi.nlm.nih.gov/pubmed/) repositories (see supplementary file 1 for details on search strings and matching with PRISMA statements). Second, materials were also retrieved following requests to specialists (e.g. Aliens mailing list, https://list.auckland.ac.nz/sympa/info/aliens-l) to bridge gaps identified for Russia and China, two of the five largest countries for which available on-line data were particularly scarce. A typical message first summarized the objectives of our research project and second, requested recipients to provide relevant material and/or suggest further contacts in this regard. On the other side, we also compiled additional materials when establishing the methodology for the project (e.g. when testing different search string combinations at initial stages of the work), from the bibliographic alerts set up by the review team. All 1417 documents obtained from this process were entered in the database, with information on the person providing the document (Table 1;30). Subsequently, 150 documents identified as not previously retrieved were considered relevant for further, full-text screening (Table 1).
Extraction of cost estimates
The Online-only Table 1 comprises all the information of InvaCost that we mention further in this article, using simple quotation marks for ‘Columns’ of the database and italic letters for the different categories within each column. The full-text of each relevant material was scrutinized for any cost estimate that could be incorporated into InvaCost30. The final stage of inclusion/exclusion took place during this data extraction. When the screened documents reported cost estimates by citing sources that were not retrieved by our literature search, whenever possible we assessed the original sources of data in order to better characterize the reported cost. These novel information sources not initially captured by our literature search were then added to the collection list (Table 1). In such cases, we provided information on all documents that were consulted to trace back the original source (‘Previous materials’). In contrast, if no original cost data were found in the cited source, the document was discarded. For all reported costs where the original source was not available or accessible, we emphasized this in a dedicated column (‘Availability’).
Then, we first extracted raw cost data, i.e. how they appear in the material in local currency (‘Raw cost estimate local currency’). When multiple cost estimates were provided for a single instance, we calculated median values (e.g. different cost estimates according to several management scenarios dedicated to the same invasive population) and collated the minimum and maximum estimates provided (columns ‘Min/Max raw cost estimate local currency’). When costs were estimated at different time and/or spatial scales in the same material, we opted to choose – when possible – those estimate(s) that summarise(s) as effectively as possible the figure(s) shown in the study. If such an estimate was not obvious to identify throughout the full-text, we extracted every relevant cost estimate. In these latter cases where several cost estimates were provided in a single study, we also collated the minimum and maximum estimates provided.
Temporal information on the costs were also retrieved: the ‘Period of estimation’ as stated in the material and hence, when possible, the ‘Probable starting/ending year’ of the period of estimation and the ‘Time range’ (year if the estimate is given yearly or for a period up to one year, period if the estimate is given for a period exceeding a year). The ‘Occurrence’ column gives the status of the cost estimate as potentially ongoing (if the cost can be expected to continue beyond the period of estimation) or one-time (if the cost was deemed as unlikely to continue). For cost estimates provided without a clear indication on the timeframe considered, or covering periods shorter than a year, we considered them with a year ‘Time Range’ and a one-time ‘Occurrence’ to avoid the risk of overestimating the duration of collated costs. The ‘Raw cost estimate’– with complementary information on the ‘Time range’, ‘Period of estimation’ and ‘Occurrence’ – can be used to estimate total costs over a given period of time. We then transformed the raw cost estimates to cost estimates per year (‘Cost estimate per year’) by dividing the raw costs with a period ‘Time Range’ by the duration of the ‘Period of estimation’ (obtained from the difference between the ‘Probable ending year’ and ‘Probable starting year’). The raw costs with a year ‘Time Range’ were reported as they are, because they are already considered at the scale of a year.
Description of cost estimates in InvaCost
Each of the cost estimates recorded was characterized by a number of information, including (a) the reference from which the cost was extracted, (b) the taxonomy of the associated species, (c) the spatial and temporal coverage of the study, (d) the typology of each cost estimate and (e) the evaluation of the reliability of the estimation method(s). For most of the variables considered in InvaCost, a non-negligible part of the cost estimates was not attributable to a single existing category due to the lack of precise information provided by the authors or because they simultaneously belong to multiple categories. In such cases, we respectively reported them as either Diverse/Unspecified or as slash-separated lists of categories (e.g. Artiodactyla/Carnivora for the ‘Order’).
Details about the nature of the information retrieved as well as the choices made to characterize each cost are synthesized in Online-only Table 1:
(a) We provided bibliographic information on each reference (e.g. ‘Reference title’, ‘Authors’, ‘Publication year’). Others specific details (e.g. abstract, journal, download link) are given in a dedicated file30 with which the columns ‘Repository’ and ‘Reference ID’ of InvaCost allow correspondence of information.
(b) We normalised and harmonised all taxonomic information on the invasive species (‘Kingdom’ to ‘Species’ level) using the GBIF.org Backbone Taxonomy35. At this stage, spelling and other taxonomic errors were corrected. While each cost extracted was generally associated with a single invasive alien species, in some cases the data was related to multiple species without the possibility of disentangling species-specific costs. In this case, we mentioned either all species concerned if explicitly indicated by the author(s), or Diverse/Unspecified if not.
(c) We dedicated seven columns to describing the impacted area according to its environment (terrestrial and/or aquatic habitats), the temporal extent as mentioned earlier (e.g. ‘Period of estimation’, ‘Time range’) and the spatial coverage from the ‘Geographic region’ (e.g., Central America, South America, Oceania-Pacific Islands) – rather than the official continent for better accuracy – down to the exact site (‘Location’) when possible. Each area was related to its country of attachment, leading to some mismatches between the ‘Geographic region’ and ‘Official country’ columns due to the existence of countries with non-contiguous overseas territories. For instance, costs found from invaders in La Réunion (a French oversea department) were attributed to Africa as ‘Geographic region’ and France as ‘Country’, while France obviously belongs to European continent.
(d) We characterised the typology of each cost mainly based on the following descriptors. The ‘Implementation’ at the moment of the cost evaluation states whether the reported cost was observed (i.e. cost actually incurred by an invasive species within its invasive distribution area) or potential (i.e. not incurred but expected cost for an invasive species beyond its actual distribution area and/or predicted over time within or beyond its actual distribution area). The ‘Acquisition method’ provides information on how the cost data was obtained, i.e. report/estimation directly obtained or derived (using inference methods) from field-based information, or extrapolation relying on computational modelling. The ‘Impacted sectors’ indicates which activity, societal or market sectors were related to the cost estimate (see Table 2 for details). The ‘Type of cost’ ranges from the economic damages and losses incurred by an invasion (e.g. value of crop losses, damage repair) to different levels of means dedicated to the management of biological invaders (e.g. control, eradication, prevention).
Table 2 The different market and/or activity sectors mentioned in InvaCost.
Full size table

(e) Lastly, we evaluated the level of ‘Reliability’ of the methodology reported by the authors to provide cost estimates (Fig. 2). Prejudging the relevance of each cost estimate is not straightforward and could suffer from a high level of subjectivity. Here, we rather aimed to evaluate in the most objective manner whether the approach used for cost estimation was documented and traceable. Hence, materials that could not be accessed for full-text investigation were conservatively considered as of low reliability. Alternatively, each cost estimate recorded from any accessible material was qualitatively assessed as of high or low reliability following a procedure depending on the ‘Type of material’ analysed (peer-reviewed article or grey material; Fig. 2). Peer-reviewed articles and official documents (e.g. institutional or governmental reports) are likely validated by experts before publication. We assumed therefore that all cost estimates collected from these materials may likely be of high reliability. Conversely, for grey materials other than official reports, the attribution to one or other of these categories (high vs low reliability) was based on specific analysis of each cost estimate. We checked whether the method estimation was fully described, independently of its comprehensiveness, i.e. if the original sources or potential assumptions were properly documented or justified, and/or the calculation methodology was explicitly demonstrated. Here, we opted for a conservative strategy that might be not optimal, as depending mostly on the nature of the publication.
Fig. 2

Decision tree approach for assessing the reliability of the method used for estimating each cost. The colour of the boxes indicates which decision was taken: green when material was deemed as of high reliability, red when material was deemed as of low reliability, blue when taking any decision needs further investigation. The intended purpose of this process was not to evaluate the quality, relevance or realism of the studies performed for providing cost estimates, but rather to assess if the methodology (i) has been reviewed and validated by peers or experts prior any publication, or (ii) if not, whether this methodology was clearly stated and demonstrated.

Full size image

Beyond the factual elements included in the descriptors from (a) to (c), those presented in (d) and (e) (to which we can add the descriptor ‘Spatial scale’) are the result of a conceptual and analytical framework created based on our own experience. This experience was gained when collecting and getting acquainted with the diversity and complexity of situations one can find behind the “economic costs” linked to biological invasions, as well as the strategies used for estimating them. We think that the different subcategories identified therein (e.g. observed vs potential costs within the descriptor ‘Implementation’) should not be aggregated to limit potential confusions in future analysis. Also, we acknowledge that the possible sub-categories of these descriptors might be improved and adapted according to the scope of future analyses made using InvaCost. We are convinced that the descriptors thus defined and categorised may strongly help in this perspective.
Standardisation of cost data
Using definitions, data and indicators provided by the World Bank Open Data and the Organisation for Economic Cooperation and Development (OECD), we expressed all retrieved costs (raw costs and costs per year) in US dollars (US$) for the year 201730 using a multi-step procedure. We provided here two ways for standardising cost estimates according to the conversion factor: one based on the market exchange rate (local currency unit per US$, calculated as an annual average), and another based on the Purchasing Power Parity (PPP, local currency unit per US$, calculated as an annual average) that is the rate of currency conversion that standardises the purchasing power of different currencies by eliminating the differences in price levels between countries. Opting for one strategy or the other for further investigation or discussion is beyond the scope of this paper and will befall on the author(s) of future analyses made using InvaCost.
We first converted the cost estimates from local currencies to US$, by dividing the cost estimate with the official market exchange rate (https://data.worldbank.org/indicator/PA.NUS.FCRF?end=2017&start=1960) corresponding to the year of the cost estimation (‘Applicable year’, that is the year of the ‘Currency’ value, but not necessarily the year of the cost occurrence). The cost obtained in US$ of that year was then converted in 2017 US$ using an inflation factor that takes into account the evolution of the value of the US$ since the year of cost estimation. The inflation factor was computed by dividing the Consumer Price Index (CPI, which is a measure of the average change over time in the prices paid by consumers for a market basket of consumer goods and services; https://data.worldbank.org/indicator/FP.CPI.TOTL?end=2017&start=1960) of 2017 by the CPI of the year of the cost estimation.
As an alternative, we also converted costs to 2017 US$ value based on PPP instead of the classical market exchange rates in the initial conversion step. PPP values were primarily collected from data provided by the World Bank (https://data.worldbank.org/indicator/PA.NUS.PPP?end=2017&start=1990), or by the OECD (https://data.oecd.org/conversion/purchasing-power-parities-ppp.htm) when information was not retrievable through the World Bank database. For this purpose, we had to deal with published costs that were expressed in currency that was different from the country where the costs were estimated (e.g. published cost in African countries expressed in US or Canadian $). Thus, prior to using PPP as a conversion index, we had to perform a retro-conversion by multiplying the original cost estimate by the official market exchange rate (local currency unit per currency unit used). For PPP-based standardisation, it was not possible to perform the process for all cost estimates as PPP data do not exist for all countries and/or specific periods (we mentioned NA in the database when such information was missing).
In summary, we used the following formula to convert and standardise each cost estimate:

$${C}_{e}=left({{boldsymbol{M}}}_{{boldsymbol{V}}}/{{boldsymbol{C}}}_{{boldsymbol{F}}}right),times ,{{boldsymbol{I}}}_{{boldsymbol{F}}}$$

with Ce = Converted cost estimate (to 2017 US dollars based on exchange rate or Purchase Power Parity), MV = Cost estimate (either the ‘Raw cost estimate local currency’ extracted from analysed paper or the ‘Cost per year local currency’ transformed by us), CF = Conversion factor (either the official market exchange rate or the purchasing power parity, in US dollars), IF = Inflation factor since the year of cost estimation, calculated as CPI2017/CPIy with CPI corresponding to the Consumer Price Index and y corresponding to the year of the cost estimation (‘Applicable year’).
We thus provided four columns with the raw cost estimates or the cost estimates per year, expressed in 2017 USD based on the exchange rate or PPP.
Data summary
InvaCost currently contains 2419 cost estimates (1215 from peer-reviewed articles, 1204 from grey materials), collected from 849 references, of which 1769 estimates were deemed as of high reliability. In total, twenty currencies are reported in our database, the majority being US dollars, n = 1348 cost estimates. Not all cost estimates were successfully converted to 2017 US$ as (i) conversion data from official sources are available only since 1960 (cost estimates range from 1945 to 2017 in InvaCost) or simply not found for some years and countries, and/or (ii) cost data are sometimes simultaneously associated with several countries, constraining the PPP-based standardisations. Hence, respectively 2416 and 2126 estimates were successfully converted using market exchange rates and PPPs. Cost estimates are either direct reports/estimations (n = 2127) or values gathered from extrapolative computations (n = 292). At a taxonomic level, these estimates are associated with 343 species belonging to six kingdoms (Animalia, Bacteria, Chromista, Fungi, Plantae, Riboviria). InvaCost has global coverage (90 countries) and includes continental, insular and overseas territories. Data are associated with terrestrial as well as aquatic (freshwater, brackish and marine) environments. Costs were estimated at different spatial scales (continental (n = 35), country (n = 1111), global (n = 17), intercontinental (n = 9), regional (n = 67), site (n = 836), unit (n = 329)). The Table 3 summarises quantitative data and information reported in InvaCost for each geographic region considered (see also Supplementary file 2).
Table 3 Quantitative summary of information recorded in InvaCost according to the ‘Geographic region’ of the cost estimates.
Full size table

Possible applications
InvaCost is expected to help bridge the gap between a growing scientific understanding of invasion impacts and still inadequate management actions. This work is thus in line with the aims of a panel of decisions recently adopted by the Convention on Biological Diversity (Decision XIII/13, https://www.cbd.int/doc/decisions/cop-13/cop-13-dec-13-en.pdf) advocating the incorporation of invasion science knowledge into management planning. In addition to offer unique opportunities for future research, InvaCost will provide a strong quantitative and evidence-based support for impacts of invasive species reported in other databases such as the Global Register of Introduced and Invasive Species (GRIIS)20, helping refine information in this database. Also, invasive populations recorded in InvaCost but data deficient in the GRIIS should be ultimately classified in that database.
Additionally, InvaCost could be considered as another data-based component, adding novel and significant information on invader impacts categorised by the Socio-Economic Impact Classification of Alien Taxa (SEICAT)36. The latter is a classification system, applicable across a broad range of taxa and spatial scales, providing a consistent procedure for translating the broad range of measures and types of impacts into ranked levels of socio-economic impacts, assigning alien taxa on the basis of the best available evidence of their documented deleterious impacts. Quantitative support provided by InvaCost will strongly contribute to impact classification. Ultimately, integrating data from these diverse sources could allow a complete description of the overall impacts of biological invasions at regional and global scales.
Caveats and directions for further database improvement
Rather than claiming exhaustiveness of data collated, we highlight that InvaCost should be considered as the most current, standardised, accurate and globally representative repository of various economic losses and expenditures documented for the largest possible set of invaders. We are aware that our database can be improved in at least three ways.
First, InvaCost mostly does not include publications and reports not yet available in electronic format and/or using non-English language, leaving open the possibility of increasing data comprehensiveness and limiting potential biases. Indeed, local reports as well as research results from some countries (e.g., China, Russia) are likely to be published in non-English language37. Again, accessing grey literature is challenging as it is not systematically digitalised and/or included in well-curated bibliographic databases29. We strongly encourage future users of InvaCost to help gathering this currently unreachable information when possible. Furthermore, some mistakes might have occurred despite our best efforts when constructing InvaCost. In this regard, we advocate for regular public updates of InvaCost in order to improve it both quantitatively (by adding currently inaccessible or missed information) and qualitatively (if errors are identified).
Second, as the distribution and impacts of invaders are inherently dynamic for a number of reasons38, InvaCost should further consider the status of the species recorded for their economic impacts in order to improve both the relevance and the usefulness of the database. As an illustration, InvaCost likely includes invasive populations currently extirpated from particular areas after successful eradication campaign(s) as well as those still established but for which impacts are locally reduced as a result of management efforts. Attempting to obtain and integrate such information into InvaCost was beyond the scope of this work. Nonetheless, it should be reciprocally beneficial to establish connections between InvaCost and other databases such as the GRIIS that provides a harmonised, open source, multi-taxon database including verified information on the continued presence of introduced and invasive species for most countries20. In light of such additional information, the value of InvaCost will be its application for policy purposes, such as identification of exotic invaders that are currently associated with economic losses in particular areas. Also, crossing information between databases may allow the refinement of the descriptor ‘Spatial scale’ we propose here.
Third, we would recommend, for a future updated version of InvaCost that would require screening back all the materials, to improve the ‘Acquisition method’, ‘Implementation’ and ‘Reliability’ descriptors, to pay attention to the specificity of “avoided costs” and to create a new descriptor for ‘non-market values’. We detail these possibilities below.
Improving descriptors
An improved version of the ‘Acquisition method’ could lead to a subdivision of the extrapolation category into spatial, temporal and spatio-temporal extrapolation. This would allow simultaneous refinement from the currently binary ‘Implementation’ descriptor (observed vs potential) into several levels of certainty regarding the incurred cost (e.g. taking into consideration the temporality (past/current or predicted) of the onset of the cost and of the status of the invasive species in the study area). The next step for deeming the ‘Reliability’ of the cost estimates recorded in InvaCost would consist of assessing the repeatability of the methodology used, by adapting the approach previously developed by Bradshaw et al.14. The latter evidenced that assuming the reproducibility of published methods should not rely only on the nature of the materials and recognized the qualitative nature of the procedure, although applying this approach to InvaCost was constrained by the large sample size and high diversity in our database (Bradshaw et al.’s study focused on a single taxonomic class). Also, because InvaCost involves several collaborators and potential future contributors, consistent and objective criteria should be further defined to cope with the large array of materials, methods and situations encountered.
Avoided costs
Introducing certain actions against biological invasions leads to avoided costs. Such avoided costs are sometimes evaluated, for instance to examine the relevance of different potential actions or to assess the effectiveness of an action that was taken. However, avoided costs cover a great variety of situations and require a careful consideration for future analysis, even if they do not have to be analysed separately from the other economic costs gathered in InvaCost. For instance, in the case of hypothetical actions, avoided costs can be considered as minimum estimates of the “real” costs (if they are unknown). However, in the case of completed or planned actions, the reported data should be the original costs (if known) minus the avoided costs, because the latter do no longer exist. Some avoided costs are probably already included in InvaCost but they are likely underestimated because keywords such as “savings” or “benefits” were not included in the search strings. Also, even if they are sometimes mentioned as “benefits” in the literature, care should be taken not to confuse these avoided-costs with the benefits incurred by direct use or exploitation of invasive species. The latter have been ignored in InvaCost since they were relatively few (and beyond of the scope of this database), but might constitute a twin project.
A new ‘Non-market values’ descriptor
The means dedicated to preventing or managing an invasion (e.g. manual removal of invasive plants) and certain economic losses and damage due to an invasion (e.g. the value of crop losses or the repair costs of damaged infrastructures due to an invasive insect) are observable on markets. However, some costs are not observable on markets but can be translated in monetary terms using several valuation methods – for instance, the willingness to pay for the conservation of a native species that is impacted by an invasive species is considered as the value given by a group of people to preserving the native species (i.e. the value that would be lost if this native species was impacted). We recognize the importance of informing the public about “non-market values”, as giving an economic value to ecosystems or biodiversity can be a way of recognising and taking them into account in public decision-making processes39, but attention should be paid to the issues linked to their assessment40,41. Among others, the different methods for assessing non-market values do not necessarily capture the same aspects of the values, so the resulting estimates might be different. Moreover, the very principle of giving a value to “benefit from nature” through economic valuation is not necessarily acknowledged by the entirety of scientific and civil communities39,42. For future analysis, the ‘non-market values’ should not be systematically aggregated with the other economic costs gathered in InvaCost. It is to note that while some non-market values are probably already included in InvaCost within the losses and damage ‘Type of cost’, the loss of non-market values is probably largely underestimated in the database because they were not the primary focus of InvaCost and therefore the related keywords were not included in the search strings.
These possible ways of improvement call for completion and/or refinement of existing entries as well as integration of newly published or acquired data by future contributors in InvaCost, with the aim to consolidate its long-term relevance (cf. Usage Note paragraph). More

Tree-ring data
We used tree-ring records from over 210,000 trees of 110 species, distributed globally in habitats ranging from the tropics to the Arctic region over more than 70,000 sites (Supplementary Fig. 1, Supplementary Table 1). The largest publicly available data source from which we used data is the International Tree-Ring Data Bank (ITRDB, https://www.ncdc.noaa.gov/data-access/paleoclimatology-data/datasets/tree-ring). These were complemented with other datasets to maximize the number of records for each species and to fill in spatial gaps. A particularly large tree-ring dataset used in our analyses is the National Forestry Inventory data from the Ministère des Forêts de la Faune et des Parcs from Quebec, Canada52,53 (hereafter NFI-Quebec). This data consists of a complete set of ring-width data from 156,711 trees from 79.381 sites across the province of Quebec. Field tree-ring data were collected according to specific standard protocols52,53, which consisted of selection of up to nine trees in each plot, with 3–5 trees ( >91 mm diameter at breast height, DBH) randomly selected, 1–2 selected from the largest trees, and 1–2 from trees closest to the mean tree diameter of the plot52,53. From selected trees one core per tree was collected. Tropical tree-ring data were compiled from the ITRDB and from unpublished and published records26,27,54,55. For species with larger sample sizes, we distinguished between tree-ring data from trees that died naturally before the moment of sampling and trees that were alive at the moment of sampling, allowing us to test the assumption that living tree ages can be used to estimate trees natural lifespans (see section “Trade-off estimates and assessment of possible artefacts”). Part of these dead tree data were obtained from the ITRDB by selecting tree-ring records of which the last measured ring width was dated to before AD1900. We assumed that most of these trees must have been dead at the time of sampling as no records were collected before 1900. In addition, we compiled published dead tree data from refs. 21,56, and used subfossil tree-ring data from refs. 57,58. Supplementary Table 1 provides an overview of the datasets, and full details of each dataset are available online as supplemental info.
Various data controls and selection procedures were used to assure high confidence in our dataset. Where possible we tried to identify duplicate records, i.e., multiple cores taken from the same individual tree. This is only a problem for ITRDB, but not for the NFI-Quebec dataset where only one core per tree existed, or for datasets from co-authors. We thus merged ITRDB records that had identical ID′s except for the last character of their ID (e.g. 01a, 01b, or ID1-1, ID1-2, etc). From the ITRDB, we only used species for which we could obtain data from a minimum of 3 different sites with at least 20 records each, and only selected species that had a minimum total of 100 separate ring width series. We excluded those sites from the ITRDB that showed relative even age structures, and are thus unlikely to represent old-growth populations that provide robust estimates of trees’ maximum lifespans. To this end, we calculated for each ITRDB site the coefficient of variation in tree ages (CVAge = StandDevAge/ MeanAge × 100) and excluded sites with a CVAge lower than 10%. A large subset of ITRDB-data from 46 species has previously been inspected for data quality by co-author S. Voelker59. In this subset of data, each ring width series was manually re-aligned by cambial age (i.e., ring number from pith), providing more reliable estimates of tree ages. For the datasets that were not acquired from the ITRDB, we used slightly different criteria. From NFI-Quebec, we used all available sites, excluding those that were classified with evidence of recent management (commercial thinning or clear-cutting) and where fire or insect disturbances destroyed more than 25% of the forest cover.
For estimates of species-level early growth rates and lifespans (cf. Fig. 1a), we included only species with a minimum of 30 records, as lower sample size is unlikely to provide good approximations of tree life spans. This resulted in the inclusion of 110 species, with a median sample size of 305 trees and 12 sites per species. To assess within-species relationships between early growth and tree lifespan, we included 82 species. As a general rule, we included only species with more than a total of 150 trees and from at least 3 sites. About half of our species had more than 300 tree records (see Supplementary Information).
To assess what minimum sample size is needed to get a representative estimate of the true maximum age of a species or a site, and to evaluate how sample size affected estimates of trade-offs between early growth and longevity, we randomly resampled 500 times varying sample sizes—from 25 to 600 trees—from a subset of 11,752 Picea mariana trees from NFI-Quebec sites located north of 50.7°N. Comparison of the maximum ages of these random subsets of trees with the true observed maximum ages shows that a sample size of 100 trees results in 99.4% of the cases in maximum age estimates larger than the 95th percentile of the original dataset, and in 67.2% of the cases in ages larger than the 99th percentile original age (Supplementary Fig. 2a, b). As more than 70% of the species had at least 100 trees we thus assumed that for most species, the estimates of their lifespan were close to true lifespans. We used this same approach to assess how sample size affected the estimate of the trade-offs (i.e., estimation of the negative exponential decay constant; see next section for details). This analysis showed that sample size of 300 trees (corresponding to median sample sizes for trade-off analysis), leads to mean errors in the estimated slope of 12% (Supplementary Fig. 2c). Thus, for most species we achieve relatively accurate estimates of the trade-off strength. Low sample sizes for some species will nevertheless result in small errors of the mean slope, but we expect that positive and negative errors will cancel out against each other. Indeed, we do not observe a specific bias towards over- or under-estimation for low sample sizes, as the mean exponential decay constant for a simulated sample size of 150 trees is very similar to that observed (i.e., −0.409 versus −0.399).
Trade-off estimates and assessment of possible artefacts
The strength of the trade-offs between growth and tree lifespan was assessed for each species using a 95th quantile regression between mean early growth rate and the natural logarithm of age using the QUANTREG package in R60, as

$$begin{array}{*{20}{c}} {{{log}}(A_{(95{mathrm{th}},{mathrm{quantile}})}) = a + b cdot overline {{mathrm{RW}}} } \ {{mathrm{or}}} \ {{mathrm{Lifespan}} approx A_{left( {95{mathrm{th}},{mathrm{quantile}}} right)} = {mathrm{exp}},left( {a + b cdot overline {{mathrm{RW}}} } right)} end{array}$$
(1)

where A is age of the tree, (overline {{mathrm{RW}}}) is the mean ring width over the first 10 years. The constant b describes the negative exponential decay constant (i.e., exponential rate of decrease of tree lifespan with increasing early growth rate). This quantile regression fit results in similar estimates of the maximum ages of trees as the 95th percentile ages in binned early growth rate categories (see Supplementary Fig. 3a–c). Note that in contrast the maximum diameter does not vary strongly between slow and fast-growing trees (Supplementary Fig. 3d). We chose a relative short period, the first ten years, for estimating early growth as our study included some relative short-lived species. Previous studies have found similar results when using longer periods (50 years)56, and we expect no substantial difference using different early growth periods as tree growth is usually strongly auto-correlated in time61.
To assess trade-off strengths within species, we calculated the mean decay constant (b, Eq. 2) for each species using relative age, A/max(A), and relative mean early ring width, (overline {{mathrm{RW}}})/max((overline {{mathrm{RW}}})). Maximum of A and (overline {{mathrm{RW}}}) are species level maxima. The mean slope calculation across all species was weighted by the cube root of the sample size to account for the large differences between species in sample size, and confidence of the trade-off estimates.
While these relationships suggest true trade-offs, they may also be affected (or even driven by) the approaches or analytical methods used here. In particular, we here evaluate the effect of the following four possible artefacts on our results; (1) the use of living trees to estimate tree lifespans, (2) effect of recent growth increases on early growth-age relationship, (3) effects of pith offsets and wood decay on early growth-age relationship, (4) sampling artefacts, such as disproportionate selection of large trees.
(1)
Use of living trees: our analysis includes mostly trees that were sampled when still alive, and may thus not be representative for the true lifespan trees may achieve. To assess to which degree use of living trees may affect our results, we analyse and compare the trade-off strengths of trees that died before 1900 to living trees for 12 species with sufficient data availability (minimum of 150 dead and 150 living trees). As the slopes of dead and living trees do not differ significantly (Supplementary Fig. 5f, g, paired t-test exponential decay coefficient, t = −0.1095, p = 0.915, n = 12), we conclude that 95th quantile regressions on living trees can be used to approximate tree lifespan.

(2)
Effect of recent growth increases: recent growth stimulation of trees due e.g. to CO2 fertilisation, warming in higher latitudes, and/or nitrogen deposition, may result in observation of a trade-off. This is because recent increases in growth will lead to higher early growth rates for young trees compared to old trees, resulting in a negative relationship between early growth and tree age. The comparison of trade-off strength of dead versus living trees provides strong evidence that this effect does not drive the trade-off. In addition to this, we use a data driven forest simulation (see section “Examining the effects of growth stimulation on forest dynamics”) to assess how growth increases affect estimations of the trade-off strength. In this simulation, we used the actual tree-ring data to simulate realistic growth increases of Picea mariana tree-ring trajectories in response to high latitude warming. By sampling from these trajectories at the end of the growth increase period (i.e., year 350), and in a period without any recent growth increases (i.e., year 600, see Fig. 3e), we establish that growth increases result in only a small over-estimation of the trade-off, decreasing the exponential decay coefficient from −0.37 to −0.44 (see Supplementary Fig. 9c). Thus, it is unlikely that recent growth stimulation is the cause for the negative relation between early growth and tree lifespan.

(3)
Pith offset: tree-ring data, especially those acquired from ITRDB, may miss the innermost sections due to incomplete cores, decayed centres, or imperfect increment borer alignment. Missing rings will result in underestimation of tree ages and inaccurate early growth rates estimation and could thus affect the estimated relationship between early growth and lifespan. However, ring widths in most species decrease with tree age and size17, and even trees showing constant wood production with age, will show decreasing ring width because of geometry. Thus early growth in these samples will underestimate true growth rates and would most likely weaken the observed trade-off, rather than strengthening it. A comparison of species present in both the NFI-Quebec and the ITRDB datasets confirms this. The NFI-Quebec dataset was less affected by pith offset problems, as the trees were carefully screened and trees with substantial differences between cumulative ring widths and field diameters were excluded. Yet, we find that slopes were more negative for NFI-Quebec compared to ITRDB (mean b of −0.25 for Quebec vs. −0.10 for ITRDB, two-sided paired t.test, t = 2.49, p = 0.047, n = 7 species) and pith offsets thus do not explain the relationship. This comparison also shows that estimates of the strength of the trade-offs between early growth and longevity inferred from ITRDB data are probably conservative, as the Quebec data can be considered to be of higher quality, and were collected according to standard protocols. In contrast, data from the ITRDB may contain incomplete series and were collected for unknown purposes, and these issues probably weaken trade-offs in the ITRDB.

(4)
Sampling biases: one potential bias in our dataset may arise due to the tendency of tree-ring studies to sample predominantly large trees in the field (i.e., big tree selection bias62,63,64). This may result in a negative relationship between early growth and tree age, as young slow-growing trees tend to be underrepresented in the tree-ring sample (i.e., have not reached the field minimum size threshold yet), compared to fast-growing young trees that are much larger, and therefore more likely to be sampled. This effect would reduce the number of trees with slow early growth and young ages in the tree-ring sample (i.e., trees in the lower left-hand corner of the early growth-lifespan graphs, cf. Supplementary Fig. 6a), and results in overestimation of the 95th percentile age estimates for slow-growing trees. Our approach to estimate to which degree this bias affected our estimates of growth-lifespan trade-offs was as follows. We first used the tree-ring NFI-Quebec data of Picea mariana, combined with plot data from Quebec to reconstruct a new artificial tree-ring dataset with a size frequency distribution identical to the population size frequency distribution for this species in Quebec (Supplementary Fig. 6b). For each tree of Picea mariana sampled for their tree-rings we know the early growth rate and age, and also the complete diameter- and age-trajectory up to the year of sampling. From these data, we resampled for each size class (in bin widths of 2 cm) the same number of trees as that observed in the field. By doing this we filled in the lacking data of trees smaller than 91 mm, and created a new artificial tree-ring dataset that had an identical size structure to that observed in the field. We know the mean growth rate over the first ten years and the age at which each individual tree reached the diameter of their respective size class, and could thus reconstruct the early growth rate versus tree age graphs for the full population, including the smaller size classes which were missing from our original tree-ring sample. We then compared the early growth -lifespan relationship for the complete population to that of the trees larger than 91 mm, mimicking the NFI-Quebec field collection protocol. This shows that the exponential decrease is marginally larger (b = −0.505 compared to −0.470 for trees >91 mm) and that the intercept is lower (159 years compared to 220 for trees >91 mm) when sampling all trees compared to only trees with diameters >91 mm (Supplementary Fig. 6). Hence, the use of a minimum size threshold (91 mm) in the NFI from Quebec results in a slight underestimation of the trade-off (by ~7%) for the Picea mariana dataset. We also resampled from this artificial dataset the 10% largest trees, to mimic a hypothetical standard tree-ring sampling scenario that only samples the largest trees. Such a sampling scenario resulted in a decay constant of −0.432, thus again causing a small underestimation of the true trade-off. This simulation proves that the trade-off is not a result of a sampling bias.

Possible environmental drivers of the trade-offs
We evaluated whether the observed trade-off between early growth and tree lifespans could be caused by covariance of growth and lifespan with climate, soil or competition. Temperature variation for example reduces tree growth and lifespan in various species (cf. Fig. 2a, b). To this end, we calculated site-level mean early growth rates and the maximum tree age for a set of species covering different geographic regions (North America, Europe and Quebec). For Quebec, we combined multiple nearby locations to obtain a minimum of 30 trees per site, as sample sizes were low for each location. Site-level mean annual temperature and precipitation was obtained from WorldClim65. We then assessed for nine different species the effect of temperature and precipitation on site-level mean early life growth and maximum tree age using major axis regression from the package smart-366. These analyses confirm that early life growth is positively related to temperature for all nine species studied, and that lifespan decreases significantly with temperature for seven out of nine species (see Fig. 2a, b). Using linear mixed effect models with species as random factor (nlme-package-R67), we find that across all nine species, mean early life growth increases on average by 0.11 mm for each degree temperature increase, while lifespan decreases by 13 years for each degree temperature increase. Precipitation has no significant effect on early life growth or tree lifespan.
To disentangle whether lifespan decreases are a direct effect of temperature increases, or due to increases in early life growth, mixed effect models were run for all nine species that simultaneously included temperature and mean early life growth rate as explanatory variables for variation in tree lifespan. To account for species differences in growth and age, we used relative mean early ring width ((overline {RW})/max((overline {RW}))), and relative maximum age (A/max(A)), and used species as a random factor with random intercepts for both early life growth and temperature. This analysis shows that mean early life growth is a stronger predictor of tree lifespans than temperature (t-value early growth = −7.2, p  , D_0} end{array}} right..$$
(4)

Here b = 0.025, k = 3 × 10−11 and D0 = 91 mm, which is the minimum sampling diameter used for the NFI-Quebec. The rationale to set the mortality below D0 to zero in this model is that the trees sampled for tree rings did all survive to this diameter (as only trees with diameters >91 mm were cored). Thus, only by setting mortality to zero for trees , A_0} end{array}} right.,$$
(5)

with a = 0.021 and b = 0.0000015. Analogous to the diameter-dependent mortality model, we set mortality to zero for trees younger than 74 years (A0), which is the age at which the average Picea mariana tree reaches 91 mm in diameter (D0).
We next created a 600 year sequence of annually seeded, 1250 member, tree cohorts. Each member of a cohort is a randomly selected tree diameter growth trajectory derived from the tree-ring cores, extended in time to an age of 500 years. Short growth trajectories were extended by using the mean growth of the 10 oldest trees that had a similar ring width in the first ten years of growth. To this end all early mean ring width was grouped into six equal early ring width classes.
All trees of this so-constructed tree cohort set live, by construction, exactly 600 years. Realistic age structures were realised by sequentially (year-on-year and tree-by-tree) assigning death to trees where a random number generator identified those individuals  smaller than μ(D) or μ(A). To test the realism of this procedure, we compared the predicted and observed tree age versus early ring width relationship—or i.e. the growth-longevity trade-off. The slope of the relationship for the diameter-dependent mortality model μ(D) is very close to observed (Supplementary Fig. 9a), and is thus a realistic representation of the observed mortality process and justifies its use to examine the effect of a growth stimulation on standing stocks. In contrast, we find that the age-dependent mortality model μ(A) does not result in a significant trade-off between early growth and tree lifespan (Supplementary Fig. 9b), providing an ideal comparison for models that fail to incorporate the observed trade-offs.
To mimic growth stimulation, we boosted growth of trajectories from year t0 = 300 year onwards of the 600 year sequence of cohorts, while exposing the trajectories over the entire 600 year period to the mortality algorithm just described. We stimulated growth rate, (overline {{mathrm{RW}}}) (mm year−1) from year t0 = 300 year onwards according to

$$overline {{mathrm{RW}}_{stim}} left( t right) = overline {{mathrm{RW}}} left( t right) cdot exp left( {lambda left( A right) cdot delta Tleft( t right)} right),$$
(6)

and

$$delta Tleft( t right) = left{ {begin{array}{*{20}{c}} {frac{{overline {dT} }}{{dt}} cdot left( {t – t_0} right),} & {t – t_0 ,150 yrs) as younger trees are more sensitive to temperature increases than older trees (Supplementary Fig. 9d), and use an exponential model of the form

$$overline {{mathrm{RW}}} left( {A,,T} right) = a cdot {mathrm{exp}}left( {lambda left( A right) cdot delta T} right),$$
(8)

where a is a constant, λ(A) is the exponential increase rate for age class [A, A+ δA] and δT = (T − 0) (°C). We then used the exponential growth increase with temperature for each age band, λ(A), to estimate the relationship between tree age and λ(A) (Supplementary Fig. 9e). In all cases the age modulation of the stimulus is

$$lambda left( A right) = left{ {begin{array}{*{20}{c}} {0.0000132 cdot A^2-0.00291 cdot A + 0.229,} & {A , < ,135} \ {0.07,} & {A , ge ,135} end{array}} right..$$ (9) Finally, we compared the effect of growth stimulation on forests dynamics against simulation without growth increases (baseline) for the two different mortality models. We also performed one simulation where we multiplied the full growth series by 2, as a representation of the effect of growth stimulation on a faster-growing species. Finally, we evaluated for each simulation the mean ring width growth, the stem mortality rate, the age of the largest (75th percentile) trees that died and the change in the total basal area stock over time for the full population. For the stem mortality and the basal area stocks, we calculate and present the change in dynamics of the growth stimulation scenario relative to the baseline scenario without growth stimulation. All analyses and simulations were performed using R-studio, version 0.99.90370. Maps in SI figures were produced using the ggmap function from the R-package ‘ggplot2’71 More