Philippa Kaur

Results of chemical analysisThe results of the study showed that the concentrations of Cd and Pb among all analyzed fruit samples (n = 242) were below the associated LOQs in only 87 and 96 samples, respectively. Similarly, in vegetable samples (n = 128) we found that Cd and Pb concentrations were below the LOQ in 31 and 69 samples, respectively. The levels of the Cd and Pb in the analyzed food samples were compared and contrasted with the maximum levels in foodstuffs regulated by legal acts: Commission Regulation (EU) No 488/2014 of 12 May 2014 amending Regulation (EC) No 1881/2006 as regards maximum levels of cadmium in foodstuffs and Commission Regulation (EU) 2015/1005 of 25 June 2015 amending Regulation (EC) No 1881/2006 as regards maximum levels of lead in certain foodstuffs3,4. It was found that in 12 food samples, the Cd content exceeded the maximum acceptable level. Among the fruit samples, this result was observed in: frozen raspberries (n = 1; 122% of maximum level) and frozen strawberries (n = 1; 114% of maximum level). In the case of vegetables, this result was observed in: fresh beetroots (n = 2; 203% and 670% of maximum level), frozen carrot (n = 1; 113% of maximum level), fresh celery (n = 4; 130%, 150%, 345%, 356% of maximum level) and processed tomatoes (n = 3; 102%, 112%, 134% of maximum level). The maximum permissible Pb level was exceeded in 3 analyzed food samples: fresh beetroot (n = 1; 135% of maximum level), frozen carrot (n = 1; 117% of maximum level) and 1 sample of frozen tomatoes in which the Pb concentration was up to 1074% of the acceptable limit (Table 5).Table 5 The number and type of food samples in which the maximum level of Cd or Pb has been exceeded.Full size tableTables 6 and 7 present the mean and SD, as well as the minimum and maximum values for the Cd and Pb contents in each of the analyzed fruits (Table 6) and vegetables (Table 7). Heavy metals concentrations were reported in mg/kg f.m. (fresh mass) in the fresh, frozen and processed products, while the content of Cd and Pb in dried products were presented in mg/kg d.w. (dry weight). Lack of a value in the tables means that the Cd or Pb value was below the LOQ for that particular sample.Table 6 The mean value, standard deviation, minimum and maximum values ​​of Cd and Pb concentrations in particular types of fruit samples.Full size tableTable 7 The mean value, standard deviation, minimum and maximum values of Cd and Pb concentrations in particular types of vegetable samples.Full size tableThe analysis of Cd and Pb contents in all food products is necessary due to the possibility of assessing the health risks associated with consumption of contaminated ready-to-eat different types of food. A review of the scientific literature showed that the issue of food contamination with heavy metals is discussed by several researchers. However, they mostly include only fresh fruits and vegetables. Additionally, there is a little data concerning the level of heavy metals contamination of vegetables and fruits cultivated in other European countries in the available literature. Consequently, the results presented in this paper may form the basis for further research on the scale of food contamination with heavy metals such as Pb and Cd.Among fruits such as apples, pears, raspberries and strawberries, the highest average values of both Cd and Pb were observed in dried products (Cd: 0.023, 0.015, 0.116, 0.131 mg/kg d.w., respectively; Pb: 0.127, 0.036, 0.111, 0.161 mg/kg d.w., respectively). In cranberry samples, the highest levels of Cd were determined in fresh fruits (0.008 mg/kg f.m.), while Pb—in processed products (0.01 mg/kg f.m.). In the case of grape samples, the same average Cd concentration was recorded in both dried and fresh products (0.001 mg/kg), while the highest Pb content was observed in processed products (0.07 mg/kg f.m.). In most fruit samples the lowest average Cd concentrations were determined in processed products (grapes, pears, raspberries and strawberries—0.0004, 0.0008, 0.009, 0.003 mg/kg f.m., respectively), while Pb—in fresh fruits (cranberries, grapes, pears—0.004, 0.005, 0.008 mg/kg f.m.) or processed (raspberries and strawberries—0.011 and 0.006 mg/kg f.m.). In apple samples, the same average Pb value was recorded in both fresh fruit and processed products (0.009 mg/kg f.m.).The content of Cd and Pb in fruits, in the results available in the literature, is very diverse. The demonstrated average Cd content in apples (0.001 mg/kg f.m.) is lower compared to studies from other regions of the world, including Great Britain (0.002 mg/kg f.m.)23. The amounts of Cd in raspberries and strawberries tested in Poland were higher compared to those investigated by Norton et al. (2015) (0.002 mg/kg f.m. vs 0.011 mg/kg f.m. and 0.002 mg/kg f.m. vs 0.018 mg/kg f.m.)23. Additionally, in samples collected in Turkey and Serbia, the Cd content in the analyzed products was below the LOQ24,25.Our results of Pb values in fruit samples are similar to those reported by some researchers and the range of values presented for this element in other analyses were very wide. However, as in the case of Cd content in apples purchased in Poland, Pb concentrations in these fruits (0.009 mg/kg f.m.) were also lower than other studies—minimum of 200%23. The average Pb content in grapes (0.009 mg/kg f.m.) was comparable to that obtained by Bağdatlıoğlu et al. (2010) (0.006 mg/kg f.m.)24. The results of author’s research regarding the content of Pb in raspberries (0.012 mg/kg f.m.) exceeded 2.5 times those published by Norton et al. (2015)23. Pb concentrations in strawberries (0.009 mg/kg f.m.) compared to other studies are in their lower range (0.010 mg/kg–0.027 mg/kg f.m.)23,24.The highest average concentrations of Cd were determined in fresh vegetables (beetroot and celery—0.235 and 0.152 mg/kg f.m., respectively) and dried—carrots and tomatoes (0.2 and 0.103 mg/kg d.w.), while Pb—in frozen vegetables (beetroots and tomatoes—0.173 and 0.294 mg/kg f.m.), as well as dried (carrots and celery—0.206 and 0.259 mg/kg d.w.). For most samples, the lowest average Cd and Pb levels were observed in processed products (beetroots, carrots, celery). Exceptions were samples of tomatoes—the lowest average Cd and Pb concentration values were observed in fresh foodstuffs (0.003 and 0.016 mg/kg f.m., respectively).Analyses conducted by other scientists indicate lower average Cd content in fresh beetroots (0.018–0.09 mg/kg f.m.)23,26 and higher by almost 600% in the case of Pb (0.58 mg/kg f.m.)26 compared to our research (Cd—0.235 mg/kg f.m.; Pb—0.095 mg/kg f.m.). Only the British study has shown lower Pb content (0.033 mg/kg f.m.)23. Our results—concentration of Cd (0.041 mg/kg f.m.) and Pb (0.027 mg/kg f.m.) in fresh carrot samples were similar to those obtained by other authors from the same territory in Poland, but also those from Great Britain, China or Brazil—Cd values ranged from 0.014 mg/kg f.m. to 0.03 mg/kg f.m., while Pb from 0.023 mg/kg f.m. to 0.971 mg/kg f.m.23,26,27,28. In the scientific literature we found only individual articles regarding celery heavy metal contamination. Guerra et al. (2012) showed 3 times lower Cd content in this vegetable—0.05 mg/kg f.m.26. The concentration of Pb in Brazilian research indicates higher content (0.47 mg/kg f.m.) than those obtained in this study (0.031 mg/kg f.m.)26. Tomatoes are the most frequently analyzed products, probably due to the easiness and simplicity of processing. Our analysis showed relatively low concentration of Cd and Pb in fresh tomatoes (Cd—0.003 mg/kg f.m.; Pb—0.016 mg/kg f.m.). In the most available scientific data Cd levels were in the range of 0.028 mg/kg f.m. to 0.033 mg/kg f.m., and Pb from 0.078 mg/kg f.m. to 0.18 mg/kg f.m.26,28. Only Norton et al. (2015) and Bagdatlioglu et al. (2010) noted lower or equal Cd and Pb values in the corresponding product23,24.Massadeh et al. (2018) in Jordan determined Pb and Cd of various canned fruits and canned vegetables including canned juice (pineapple), canned tomato sauce, canned whole carrots and canned green beans. They showed metal concentration levels in the samples were in the range of 0.50–0.60 mg/kg f.m. for Cd and 2.6–3.0 mg/kg f.m. for Pb29. These results significantly exceed the values shown in present study, as well as the results presented by Domagała-Świątkiewicz and Gąstoł (2012) in the analysis of vegetable juices (beetroot, carrot, celery)30.The high contamination found in vegetables might be closely related to the pollutants in irrigation water, farm soil, fertilizers and also industrial and low pollution household emissions. Differences in levels of contamination between fruits and vegetables may result from the specificity of the geographical area from which they are collected, their diverse capacity to accumulate heavy metals, as well as the way they are processed. It should be pointed out that in polluted environments (soil, water, and air), the presence of toxic metals in elevated concentrations is not uncommon. Due to the structure of consumption of various groups of food products both in Poland and other countries, a significant risk of exposure to heavy metals is associated with the consumption of fruits and vegetables, which are one of the main elements of the diet. Unfortunately, complete elimination of elements such as Cd or Pb from these products is impossible, and the technological processes used in food production can only remove a small part of the impurities from selected products or even contribute to their increased contamination. Thus, there is a need for regular monitoring of heavy metals on every kind of foodstuff, not only in fresh products, in order to estimate the health risk from heavy metals in the human food chain.Statistical analysisANOVAFor the purpose of ANOVA carried out to detect significant differences in the heavy metal concentrations of the four types of food (fresh, dried, frozen, and processed), samples with concentration value below the LOQ were removed from the analysis. In the case of Cd concentration, the value of F statistic was 11.15 for fruits and 4.049 for vegetables, leading to significant results with p-values below 0.001 and 0.01 respectively. For the of Pb concentration, the ANOVA results were even more extreme with F values of 56.59 for fruits and 7.13 for vegetables with associated p-values being below 0.001 in both cases. These results show that there is strong evidence to believe that mean Cd and Pb contents in the four types of fruits and vegetables are not equal (Table 8).Table 8 Analysis of variance (ANOVA) for variates in four groups.Full size tableOutlier analysisThe boxplots depicted in Fig. 1 were used to illustrate the outlier analysis for Cd and Pb. Each plot shows the median of the observations along with the lower quartile (Q1) and the upper quartile (Q3). The highest and the lowest observations are shown by the whiskers. From Fig. 1a, there appears to be two outliers in the dried fruits with values 0.277 and 0.210. From Fig. 1b, there seems to be six outliers in the fresh vegetables with values of 0.203, 0.670, 0.260, 0.690, 0.300 and 0.712. In Fig. 1c, we see two outliers in the processed fruits with values of 0.127 and 0.047. Finally, Fig. 1d shows that there is one one outlier in the frozen vegetable category with the value of 0.537.Figure 1Outlier analysis in case: Cd concentration in: (a) fruits, (b) vegetables, and Pb concentration in: (c) fruits, (d) vegetables.Full size imageOutliers associated with high Cd and Pb values in fruit and vegetable samples may be the result of sample contamination during technological processes or vegetables/fruits cultivation in a polluted agricultural area.Post-hoc multiple comparisonSince the ANOA results indicated significant differences among the mean concentrations of Cd and Pb both in fruits and vegetables, to further detect the specific different means, the Tukey HSD test22 was applied. The results are presented in Fig. 2. For the Cd concentration, comparison of all pairs of means indicated that the content of Cd in dried fruits is significantly different from mean concentrations of other types of food namely fresh, frozen, and processed fruits, see Fig. 2a. In the case of vegetables, the mean Cd contents of fresh and processed vegetables are different, see Fig. 2b, although mean Cd content of frozen and fresh vegetables are also significantly different if a significance level of 10% is used. Upon analyzing the mean concentrations of Pb in fruits, we found that the mean content of dried fruits was significantly different from the other three types, namely fresh, frozen and processed, see Fig. 2c. For the Pb concentrations in vegetables, a highly significant difference was detected between the means of processed and dried vegetables. In addition, mean Pb concentrations of fresh versus dried and processed versus frozen vegetables were significantly different, see Fig. 2d.Figure 2Post-hoc Multiple Comparison Tukey-Test of Cd and Pb in all samples of fruits and vegetables; differences in Cd mean concentration of: (a) fruits, (b) vegetables; differences in Pb mean concentration of: (c) fruits, (d) vegetables.Full size image More

1.Stearns, S.C. The Evolution of Life Histories. https://doi.org/10.1046/j.1420-9101.1993.6020304.x (Oxford University Press/Wiley, 1992).2.Dingle, H. The evolution of life histories. in Population Biology (Wöhrmann K., Jain S.K. eds.). https://doi.org/10.1007/978-3-642-74474-7_9 (Springer, 1990).3.Yang, L., Walck, J. L. & El-Kassaby, Yousry, A. Roles of the environment in plant life-history tradeoffs. in Advances in Seed Biology (Jimenez-Lopez, J.C. ed.) 674. https://doi.org/10.5772/intechopen.70312.4.Hamrick, J. L. & Godt, M. J. W. Effects of life history traits on genetic diversity in plant species. Philos. Trans. R. Soc. Lond.. Ser. B Biol. Sci. 351, 1291–1298 (1996).5.Hobbs, R. J., Higgs, E. & Harris, J. A. Novel ecosystems: Implications for conservation and restoration. Trends Ecol. Evol. 24, 599–605 (2006).Article 

Google Scholar 
6.Prozherina, N., Freiwald, V., Rousi, M. & Oksanen, E. Interactive effect of springtime frost and elevated ozone on early growth, foliar injuries and leaf structure of birch (Betula pendula). New Phytol. 159, 623–636 (2003).CAS 
Article 

Google Scholar 
7.Zvereva, Elena, L., Roitto, M. & Kozlov, M. V. Growth and reproduction of vascular plants in polluted environments: A synthesis of existing knowledge. Environ. Rev. 18, 355–367 (2010).8.Niinemets, Ü. Responses of forest trees to single and multiple environmental stresses from seedlings to mature plants: Past stress history, stress interactions, tolerance and acclimation. For. Ecol. Manag. 260, 1623–1639 (2010).Article 

Google Scholar 
9.Possen, B. J. H. M. et al. Adaptability of birch (Betula pendula Roth) and aspen (Populus tremula L.) genotypes to different soil moisture conditions. For. Ecol. Manag. 262, 1387–1399 (2011).10.Řehounková, K. & Prach, K. Life-history traits and habitat preferences of colonizing plant species in long-term spontaneous succession in abandoned gravel–sand pits. Basic Appl. Ecol. 11, 45–53 (2010).Article 

Google Scholar 
11.Franiel, I. The Biology and Ecology of Betula pendula Roth on Post-Industrial Waste Dumping Grounds: The Variability Range of Life History Traits (Silesia University, 2012).
Google Scholar 
12.Kompała-Ba̧ba, A. & Ba̧ba, W. Participation of grasses (Poaceae) in the communities, which developed on iron smelter affected lands in the Silesian Uplands. Fragm. Florist. Geobot. Pol. 20, 267–284 (2013).13.Kompała-Bąba, A. & Bąba, W. The spontaneous succession in a sand-pit—The role of life history traits and species habitat preferences. Polish J. Ecol. 61, 13–22 (2013).
Google Scholar 
14.Darbah, J. N. T. et al. Impacts of elevated atmospheric CO2 and O3 on paper birch (Betula papyrifera): Reproductive fitness. Sci. World J. 7, 240–246 (2007).Article 

Google Scholar 
15.Alvarez-Cansino, L., Zunzunegui, M., Díaz Barradas, M. C. & Esquivias, M. P. Gender-specific costs of reproduction on vegetative growth and physiological performance in the dioecious shrub Corema album. Ann. Bot. 106, 989–998 (2010).16.Zvereva, E. L., Roitto, M. & Kozlov, M. V. Growth and reproduction of vascular plants in polluted environments: A synthesis of existing knowledge. Environ. Rev. 18, 355–367 (2010).CAS 
Article 

Google Scholar 
17.Franiel, I. & Babczyńska, A. The growth and reproductive effort of Betula pendula Roth in a heavy-metals polluted area. Polish J. Environ. Stud. 20, 1097–1101 (2011).CAS 

Google Scholar 
18.Lancaster, L. T., Morrison, G. & Fitt, R. N. Life history trade-offs, the intensity of competition, and coexistence in novel and evolving communities under climate change. Philos. Trans. R. Soc. B Biol. Sci. 372, (2017).19.Obeso, J. R. Costs of reproduction in ilex aquifolium: Effects at tree, branch and leaf levels. J. Ecol. 85, 159–166 (1997).Article 

Google Scholar 
20.Obeso, J. R. The costs of reproduction in plants. New Phytol. 155, 321–348 (2002).Article 

Google Scholar 
21.Cipollini, M. L. & Stiles, E. W. Costs of reproduction in Nyssa sylvatica: Sexual dimorphism in reproductive frequency and nutrient flux. Oecologia 86, 585–593 (1991).ADS 
Article 

Google Scholar 
22.Seidling, W., Starfinger, U. & Stöcklin, J. Plant population ecology. Prog. Bot. 55, 345–370 (1994).Article 

Google Scholar 
23.Possen, B. J. H. M. et al. Variation in 13 leaf morphological and physiological traits within a silver birch ( Betula pendula ) stand and their relation to growth. Can. J. For. Res. 44, 657–665 (2014).Article 

Google Scholar 
24.Körner, C. Limitation and stress—Always or never?. J. Veg. Sci. 14, 141–143 (2003).
Google Scholar 
25.Giuliani, C., Lazzaro, L., Mariotti Lippi, M., Calamassi, R. & Foggi, B. Temperature-related effects on the germination capacity of black locust (Robinia pseudoacacia L., Fabaceae) seeds. Folia Geobot. 50, 275–282 (2015).26.Wolkovich, E. M. et al. Warming experiments underpredict plant phenological responses to climate change. Nature https://doi.org/10.1038/nature11014 (2012).Article 
PubMed 

Google Scholar 
27.Koski, V. & Rousi, M. A review of the promises and constraints of breeding silver birch (Betula pendula Roth) in Finland. For. Int. J. For. Res. 78, 187–198 (2005).
Google Scholar 
28.Marguí, E., Queralt, I., Carvalho, M. L. & Hidalgo, M. Assessment of metal availability to vegetation (Betula pendula) in Pb–Zn ore concentrate residues with different features. Environ. Pollut. 145, 179–184 (2007).Article 

Google Scholar 
29.Hynynen, J. et al. Silviculture of birch (Betula pendula Roth and Betula pubescens Ehrh.) in northern Europe. Forestry 83, 103–119 (2010).30.Frouz, J. et al. Development of canopy cover and woody vegetation biomass on reclaimed and unreclaimed post-mining sites. Ecol. Eng. 84, 233–239 (2015).Article 

Google Scholar 
31.Řehounková, K., Lencová, K. & Prach, K. Spontaneous establishment of woodland during succession in a variety of central European disturbed sites. Ecol. Eng. 111, 94–99 (2018).Article 

Google Scholar 
32.Franiel, I. & Więski, K. Leaf features of silver birch (Betula pendula Roth). Variability within and between two populations (uncontaminated vs Pb-contaminated and Zn-contaminated site). Trees 19, 81–88 (2005).33.Kozlov, M. V. & Zvereva, E. L. Industrial barrens: Extreme habitats created by non-ferrous metallurgy. Rev. Environ. Sci. Bio/Technol. 6, 231–259 (2007).CAS 
Article 

Google Scholar 
34.Zvereva, E. L. & Kozlov, M. V. Growth and reproduction of dwarf shrubs, Vaccinium myrtillus and V. vitis-idaea, in a severely polluted area. Basic Appl. Ecol. https://doi.org/10.1016/j.baae.2004.11.003 (2005).35.Kozlov, M. V. Pollution resistance of mountain birch, Betula pubescens subsp. czerepanovii, near the copper–nickel smelter: Natural selection or phenotypic acclimation? Chemosphere 59, 189–197 (2005).36.Kondracki, J. Geografia fizyczna Polski. (Physical Geography of Poland). (PWN (in Polish), 2001).37.Pełka-Gościniak, J. Environmental aspects of relief transformation (Silesian Upland, Southern Poland). Environ. Socio-Econ. Stud. 2, 13–20.38.Tomusiak, R. et al. Age tables for silver birch (Betula pendula Roth) trees for early succession stands on abandoned agricultural lands. Sylwan 158, 579–589 (2014).
Google Scholar 
39.Szymkiewicz, B. Tablice Zasobności i Przyrostu Drzewostanów Ważniejszych Gatunków Drzew Leśnych. (Państwowe Wydawnictwo Rolnicze i Leśne, 2001).40.Dubois, H., Verkasalo, E. & Claessens, H. Potential of Birch (Betula pendula Roth and B. pubescens Ehrh.) for forestry and forest-based industry sector within the changing climatic and socio-economic context of Western Europe. Forests 11, 336 (2020).41.Ostrowska, A., Gawliński, S. & Szczubiałka, Z. Metody Analizy i Oceny Właściwości Gleb i Roślin. Katalog (Methods for Analysis and Assessment of Soil and Plant Properties) Catalog. (Instytut Ochrony Środowiska, 1991).42.Baskin, C. C. & Baskin, J. M. Seeds : Ecology, Biogeography, and Evolution of Dormancy and Germination. (Elsevier, 2014).43.Soltani, E., Ghaderi-Far, F., Baskin, C. C. & Baskin, J. M. Problems with using mean germination time to calculate rate of seed germination. Aust. J. Bot. 63, 631–635 (2015).Article 

Google Scholar 
44.Flores, P., Poggi, D., García, S. M., Catraro, M. & Gariglio, N. Effects of pre-stratification storage conditions on black walnut seed post-stratification germination capacity. Int. J. Fruit Sci. 17, 29–40 (2017).Article 

Google Scholar 
45.Ranal, M. A., de Santana, D. G., Ferreira, W. R. & Mendes-Rodrigues, C. Calculating germination measurements and organizing spreadsheets. Rev. Bras. Botân. 32, 849–855 (2009).Article 

Google Scholar 
46.Matthews, S. & Khajeh Hosseini, M. Mean germination time as an indicator of emergence performance in soil of seed lots of maize (Zea mays). Seed Sci. Technol. 34, 339–347 (2006).47.Ter Braak, C. J. F. & Šmilauer, P. Canoco Reference Manual and User’s Guide: Software for Ordination, Version 5.0. (Microcomputer Power, 2012).48.Inc., D. Dell Statistica (Data Analysis Software System, Version 13). (2016).49.Team, R. D. C. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2008).50.Nicotra, A. B. et al. Plant phenotypic plasticity in a changing climate. Trends Plant Sci. 15, 684–692 (2010).CAS 
Article 

Google Scholar 
51.Gray, S. B. & Brady, S. M. Plant developmental responses to climate change. Dev. Biol. 419, 64–77 (2016).CAS 
Article 

Google Scholar 
52.Dahl, Å. E. & Fredrikson, M. The timetable for development of maternal tissues sets the stage for male genomic selection in Betula pendula (Betilaceae). Am. J. Bot. 83, 895–902 (1996).Article 

Google Scholar 
53.Kozlov, M. V. & Zvereva, E. L. Reproduction of mountain birch along a strong pollution gradient near Monchegorsk, Northwestern Russia. Environ. Pollut. 132, 443–451 (2004).CAS 
Article 

Google Scholar 
54.Sultan, S. E. Phenotypic plasticity in plants: A case study in ecological development. Evol. Dev. 5, 25–33 (2003).Article 

Google Scholar 
55.Zvereva, E. L. & Kozlov, M. V. Effects of pollution-induced habitat disturbance on the response of willows to simulated herbivory. J. Ecol. 89, 21–30 (2001).Article 

Google Scholar 
56.Eränen, J. K. Local adaptation of mountain birch to heavy metals in subarctic industrial barrens. For. Snow Landsc. Res. 80, 161–167 (2006).
Google Scholar 
57.Neuvonen, S., Nyyssonen, T., Ranta, H. & Kiilunen, S. Simulated acid rain and the reproduction of mountain birch [Betula pubescens ssp. tortuosa (Ledeb.) Nyman]: A cautionary tale. New Phytol. 118, 111–117 (1991).58.Cuinica, L. G., Abreu, I., Gomes, C. R. & Esteves da Silva, J. C. G. Exposure of Betula pendula Roth pollen to atmospheric pollutants CO, O3 and SO2. Grana 52, 299–304 (2013).59.Pasonen, H.-L., Pulkkinen, P. & Kärkkäinen, K. Genotype-environment interactions in pollen competitive ability in an anemophilous tree, Betula pendula Roth. Theor. Appl. Genet. 105, 465–473 (2002).Article 

Google Scholar 
60.Sarvas, R. On the flowering of birch and the quality of seed crop. Commun. Inst. For. Fenn. 1–38 (1952).61.Midmore, E., McCartan, S., Jinks, R. & Cahalan, C. Using thermal time models to predict germination of five provenances of silver birch (Betula pendula Roth) in southern England. Silva Fenn. 49, 1–12 (2015).62.Bojarczuk, K. et al. Effect of polluted soil and fertilisation on growth and physiology of Silver Birch (Betula pendula Roth.) seedlings. Polish J. Environ. Stud. 11, 483–492 (2002).63.Wierzbicka, M. & Rostański, A. Microevolutionary changes in ecotypes of calamine waste heap vegetation near Olkusz, Poland: A review. Acta Bot. Cracoviensia 44, 7–10 (2002).
Google Scholar 
64.Mahdi, T. & Whittaker, J. B. Do birch trees (Betula pendula) grow better if foraged by wood ants?. J. Anim. Ecol. 62, 101–116 (1993).Article 

Google Scholar 
65.Barantal, S. et al. Contrasting effects of tree species and genetic diversity on the leaf—Miner communities associated with silver birch. Oecologia 189, 687–697 (2019).ADS 
Article 

Google Scholar 
66.Santamour, Frank, S. & Greene, A. European hornet damage to ash and birch trees. J. Arboric. 12, 273–279 (1986).67.Klingeman, B., Oliver, J. & F., H. Who’s doin’ all that chewin’? The European hornet. Tenn. Green Times 2, 34–36 (2001).68.Ylioja, T., Roininen, H., Heinonen, J. & Rousi, M. Susceptibility of Betula pendula clones to Phytobia betulae, a dipteran miner of birch stems. Can. J. For. Res. 30, 1824–1829 (2000).Article 

Google Scholar 
69.Franiel, I. Development of Betula pendula seedlings growing on the Silesia Steelworks dumping grounds in Katowice. Acta Agrophys. 51, 51–55 (2001).
Google Scholar  More

1.Díaz et al. Summary for Policymakers of the Global Assessment.pdf.2.Kehoe, L. et al. Biodiversity at risk under future cropland expansion and intensification. Nat. Ecol. Evolut. 1, 1129–1135 (2017).Article 

Google Scholar 
3.Hendershot, J. N. et al. Intensive farming drives long-term shifts in avian community composition. Nature 579, 393–396 (2020).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
4.Bommarco, R., Kleijn, D. & Potts, S. G. Ecological intensification: Harnessing ecosystem services for food security. Trends Ecol. Evol. 28, 230–238 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
5.Michael, D. R., Wood, J. T., O’Loughlin, T. & Lindenmayer, D. B. Influence of land sharing and land sparing strategies on patterns of vegetation and terrestrial vertebrate richness and occurrence in Australian endangered eucalypt woodlands. Agr. Ecosyst. Environ. 227, 24–32 (2016).Article 

Google Scholar 
6.Tittonell, P. Ecological intensification of agriculture—Sustainable by nature. Curr. Opin. Environ. Sustain. 8, 53–61 (2014).Article 

Google Scholar 
7.Willer, E. H., Schlatter, B., Trávní, J., Kemper, L. & Lernoud, J. The World of Organic Agriculture Statistics and Emerging Trends 2020. 337.8.Reganold, J. P. & Wachter, J. M. Organic agriculture in the twenty-first century. Nat. Plants 2 (2016).9.Connor, D. J. Organic agriculture cannot feed the world. Field Crop Res. 106, 187–190 (2008).Article 

Google Scholar 
10.Seufert, V. & Ramankutty, N. Many shades of gray—The context-dependent performance of organic agriculture. Sci. Adv. 3, e1602638 (2017).11.Smith, O. M. et al. Landscape context affects the sustainability of organic farming systems. Proc. Natl. Acad. Sci. 117, 2870–2878 (2020).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
12.Bengtsson, J., Ahnström, J. & Weibull, A.-C. The effects of organic agriculture on biodiversity and abundance: A meta-analysis: Organic agriculture, biodiversity and abundance. J. Appl. Ecol. 42, 261–269 (2005).Article 

Google Scholar 
13.Tuck, S. L. et al. Land-use intensity and the effects of organic farming on biodiversity: A hierarchical meta-analysis. J. Appl. Ecol. 51, 746–755 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
14.Lichtenberg, E. M. et al. A global synthesis of the effects of diversified farming systems on arthropod diversity within fields and across agricultural landscapes. Glob. Change Biol. 23, 4946–4957 (2017).ADS 
Article 

Google Scholar 
15.Lori, M., Symnaczik, S., Mäder, P., De Deyn, G. & Gattinger, A. Organic farming enhances soil microbial abundance and activity—A meta-analysis and meta-regression. PLOS ONE 12, e0180442 (2017).16.Kleijn, D., Rundlöf, M., Scheper, J., Smith, H. G. & Tscharntke, T. Does conservation on farmland contribute to halting the biodiversity decline?. Trends Ecol. Evol. 26, 474–481 (2011).PubMed 
Article 
PubMed Central 

Google Scholar 
17.Birkhofer, K., Ekroos, J., Corlett, E. B. & Smith, H. G. Winners and losers of organic cereal farming in animal communities across Central and Northern Europe. Biol. Cons. 175, 25–33 (2014).Article 

Google Scholar 
18.Mackie, K. A., Müller, T., Zikeli, S. & Kandeler, E. Long-term copper application in an organic vineyard modifies spatial distribution of soil micro-organisms. Soil Biol. Biochem. 65, 245–253 (2013).CAS 
Article 

Google Scholar 
19.Buchholz, J. et al. Soil biota in vineyards are more influenced by plants and soil quality than by tillage intensity or the surrounding landscape. Sci. Rep. 7 (2017).20.Hole, D. G. et al. Does organic farming benefit biodiversity?. Biol. Cons. 122, 113–130 (2005).Article 

Google Scholar 
21.Power, A. G. Ecosystem services and agriculture: Tradeoffs and synergies. Philos. Trans. R. Soc. B Biol. Sci. 365, 2959–2971 (2010).Article 

Google Scholar 
22.Peigné, J. et al. Earthworm populations under different tillage systems in organic farming. Soil Tillage Res. 104, 207–214 (2009).Article 

Google Scholar 
23.Biondi, A., Desneux, N., Siscaro, G. & Zappalà, L. Using organic-certified rather than synthetic pesticides may not be safer for biological control agents: Selectivity and side effects of 14 pesticides on the predator Orius laevigatus. Chemosphere 87, 803–812 (2012).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
24.Mehrabi, Z., Seufert, V., Ramankutty, N. The conventional versus alternative agricultural divide: A response to Garibaldi et al. Trends Ecol. Evolut. 32, 720–721 (2017).25.Tscharntke, T., Klein, A. M., Kruess, A., Steffan-Dewenter, I. & Thies, C. Landscape perspectives on agricultural intensification and biodiversity – ecosystem service management. Ecol. Lett. 8, 857–874 (2005).Article 

Google Scholar 
26.Gámez-Virués, S. et al. Landscape simplification filters species traits and drives biotic homogenization. Nat. Commun. 6 (2015).27.Holzschuh, A., Steffan-Dewenter, I. & Tscharntke, T. Agricultural landscapes with organic crops support higher pollinator diversity. Oikos 117, 354–361 (2008).Article 

Google Scholar 
28.Muneret, L., Auriol, A., Thiéry, D. & Rusch, A. Organic farming at local and landscape scales fosters biological pest control in vineyards. Ecol. Appl. 29, e01818 (2019).29.Gabriel, D. et al. Scale matters: The impact of organic farming on biodiversity at different spatial scales: Scale matters in organic farming. Ecol. Lett. 13, 858–869 (2010).PubMed 
Article 
PubMed Central 

Google Scholar 
30.Agreste. Pratiques Phytosanitaires en Viticulture. Campagne 2016. (2020)31.Agreste. La Viticulture Bio en Nouvelle-Aquitaine: Un Dynamisme à Tous les Stades de la Filière. (2020).32.Gruber, S. & Claupein, W. Effect of tillage intensity on weed infestation in organic farming. Soil Tillage Res. 105, 104–111 (2009).Article 

Google Scholar 
33.Pfingstmann, A. et al. Contrasting effects of tillage and landscape structure on spiders and springtails in vineyards. Sustainability 11, 2095 (2019).Article 

Google Scholar 
34.Dittmer, S. & Schrader, S. Longterm effects of soil compaction and tillage on Collembola and straw decomposition in arable soil. Pedobiologia 44, 527–538 (2000).Article 

Google Scholar 
35.Kolb, S., Uzman, D., Leyer, I., Reineke, A. & Entling, M. H. Differential effects of semi-natural habitats and organic management on spiders in viticultural landscapes. Agric. Ecosyst. Environ. 287, 106695 (2020).36.Birkhofer, K. et al. Relationships between multiple biodiversity components and ecosystem services along a landscape complexity gradient. Biol. Cons. 218, 247–253 (2018).Article 

Google Scholar 
37.Kratschmer, S. et al. Tillage intensity or landscape features: What matters most for wild bee diversity in vineyards?. Agric. Ecosyst. Environ. 266, 142–152 (2018).Article 

Google Scholar 
38.Ullmann, K. S., Meisner, M. H. & Williams, N. M. Impact of tillage on the crop pollinating, ground-nesting bee, Peponapis pruinosa in California. Agric. Ecosyst. Environ. 232, 240–246 (2016).Article 

Google Scholar 
39.Jiang, X., Wright, A. L., Wang, X. & Liang, F. Tillage-induced changes in fungal and bacterial biomass associated with soil aggregates: A long-term field study in a subtropical rice soil in China. Appl. Soil. Ecol. 48, 168–173 (2011).Article 

Google Scholar 
40.Zuber, S. M. & Villamil, M. B. Meta-analysis approach to assess effect of tillage on microbial biomass and enzyme activities. Soil Biol. Biochem. 97, 176–187 (2016).CAS 
Article 

Google Scholar 
41.Luff, M. L. The biology of the ground beetle Harpalus rufipes in a strawberry field in Northumberland. Ann. Appl. Biol. 94, 153–164 (1980).Article 

Google Scholar 
42.Shearin, A. F., Reberg-Horton, S. C. & Gallandt, E. R. Direct effects of tillage on the activity density of ground beetle (Coleoptera: Carabidae) weed seed predators. Environ. Entomol. 36, 1140–1146 (2007).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
43.Rundlöf, M. et al. Seed coating with a neonicotinoid insecticide negatively affects wild bees. Nature 521, 77–80 (2015).ADS 
PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
44.Martin, E. A. et al. The interplay of landscape composition and configuration: new pathways to manage functional biodiversity and agroecosystem services across Europe. Ecol. Lett. 22, 1083–1094 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
45.Goded, S., Ekroos, J., Azcárate, J. G., Guitián, J. A. & Smith, H. G. Effects of organic farming on plant and butterfly functional diversity in mosaic landscapes. Agric. Ecosyst. Environ. 284, 106600 (2019).46.Rusch, A., Valantin-Morison, M., Sarthou, J.-P. & Roger-Estrade, J. Multi-scale effects of landscape complexity and crop management on pollen beetle parasitism rate. Landsc. Ecol. 26, 473–486 (2011).Article 

Google Scholar 
47.Tamburini, G., De Simone, S., Sigura, M., Boscutti, F. & Marini, L. Conservation tillage mitigates the negative effect of landscape simplification on biological control. J. Appl. Ecol. 53, 233–241 (2016).Article 

Google Scholar 
48.Le Féon, V. et al. Intensification of agriculture, landscape composition and wild bee communities: A large scale study in four European countries. Agric. Ecosyst. Environ. 137, 143–150 (2010).Article 

Google Scholar 
49.Sousa, J. P. et al. Changes in Collembola richness and diversity along a gradient of land-use intensity: A pan European study. Pedobiologia 50, 147–156 (2006).Article 

Google Scholar 
50.Vanbergen, A. J. et al. Scale-specific correlations between habitat heterogeneity and soil fauna diversity along a landscape structure gradient. Oecologia 153, 713–725 (2007).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
51.Lehmitz, R., Russell, D., Hohberg, K., Christian, A. & Xylander, W. E. R. Active dispersal of oribatid mites into young soils. Appl. Soil. Ecol. 55, 10–19 (2012).Article 

Google Scholar 
52.Concepción, E. D., Díaz, M. & Baquero, R. A. Effects of landscape complexity on the ecological effectiveness of agri-environment schemes. Landsc. Ecol. 23, 135–148 (2008).Article 

Google Scholar 
53.Tscharntke, T. et al. Landscape moderation of biodiversity patterns and processes – Eight hypotheses. Biol. Rev. 87, 661–685 (2012).PubMed 
Article 
PubMed Central 

Google Scholar 
54.Desneux, N., Decourtye, A. & Delpuech, J.-M. The sublethal effects of pesticides on beneficial arthropods. Annu. Rev. Entomol. 52, 81–106 (2007).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
55.Naveed, M. et al. Simultaneous loss of soil biodiversity and functions along a copper contamination gradient: When soil goes to sleep. Soil Sci. Soc. Am. J. 78, 1239–1250 (2014).ADS 
Article 
CAS 

Google Scholar 
56.Eijsackers, H., Beneke, P., Maboeta, M., Louw, J. P. E. & Reinecke, A. J. The implications of copper fungicide usage in vineyards for earthworm activity and resulting sustainable soil quality. Ecotoxicol. Environ. Saf. 62, 99–111 (2005).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
57.Le Provost, G. et al. Land-use history impacts functional diversity across multiple trophic groups. Proc. Natl. Acad. Sci. 117, 1573–1579 (2020).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
58.Muneret, L. et al. Organic farming expansion drives natural enemy abundance but not diversity in vineyard-dominated landscapes. Ecol. Evol. https://doi.org/10.1002/ece3.5810 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
59.Lechenet, M., Dessaint, F., Py, G., Makowski, D. & Munier-Jolain, N. Reducing pesticide use while preserving crop productivity and profitability on arable farms. Nat. Plants 3 (2017).60.Le Féon, V. et al. Solitary bee abundance and species richness in dynamic agricultural landscapes. Agric. Ecosyst. Environ. 166, 94–101 (2013).Article 

Google Scholar 
61.McCravy, K. & Ruholl, J. Bee (Hymenoptera: Apoidea) diversity and sampling methodology in a midwestern USA deciduous forest. Insects 8, 81 (2017).PubMed Central 
Article 

Google Scholar 
62.Bano, R. & Roy, S. Extraction of Soil Microarthropods: A Low Cost Berlese-Tullgren Funnels Extractor. 4.63.Grueber, C. E., Nakagawa, S., Laws, R. J. & Jamieson, I. G. Multimodel inference in ecology and evolution: Challenges and solutions: Multimodel inference. J. Evol. Biol. 24, 699–711 (2011).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
64.Burnham, K. P. & Anderson, D. R. Multimodel inference: Understanding AIC and BIC in model selection. Sociol. Methods Res. 33, 261–304 (2004).MathSciNet 
Article 

Google Scholar 
65.Bartoń, K. MuMIn: Multi-Model Inference. R Package Version 1.43.17. https://CRAN.R-project.org/package=MuMIn (2020).66.R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/ (2020).67.Hartig, F. DHARMa: Residual Diagnostics for Hierarchical (Multi-Level/Mixed) Regression Models (2020). R Package Version 0.3.3.0. https://CRAN.R-project.org/package=DHARMa68.Bates, D., Maechler, M., Bolker, B., Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01 (2015). More

1.Mutreja, A. et al. Evidence for several waves of global transmission in the seventh cholera pandemic. Nature 477, 462–465 (2011).CAS 
Article 

Google Scholar 
2.Ali, M., Nelson, A. R., Lopez, A. L. & Sack, D. A. Updated global burden of cholera in endemic countries. PLoS Negl. Trop. Dis. 9, e0003832 (2015).Article 

Google Scholar 
3.Bompangue, D. N. et al. Lakes as source of cholera outbreaks, Democratic Republic of Congo. Emerg. Infect. Dis. 14, 798–800 (2008).Article 

Google Scholar 
4.Weill, F. X. et al. Genomic history of the seventh pandemic of cholera in Africa. Science 358, 785–789 (2017).CAS 
Article 

Google Scholar 
5.Hounmanou, Y. M. G. et al. Genomic insights into Vibrio cholerae O1 responsible for cholera epidemics in Tanzania between 1993 and 2017. PLoS Neglect Trop. D. 13, e0007934 (2019).Article 

Google Scholar 
6.Colwell, R. R. Global climate and infectious disease: the cholera paradigm. Science 274, 2025–2031 (1996).CAS 
Article 

Google Scholar 
7.Singleton, F., Attwell, R., Jangi, M. & Colwell, R. Influence of salinity and organic nutrient concentration on survival and growth of Vibrio choleraein aquatic microcosms. Appl. Environ. Microbiol. 43, 1080–1085 (1982).CAS 
Article 

Google Scholar 
8.Kirschner, A. K. T. et al. Rapid growth of planktonic Vibrio cholerae non-O1/non-O139 strains in a large alkaline lake in Austria: Dependence on temperature and dissolved organic carbon quality. Appl. Environ. Microb. 74, 2004–2015 (2008).CAS 
Article 

Google Scholar 
9.Reid, P. C. et al. The continuous plankton recorder: concepts and history, from Plankton Indicator to undulating recorders. Prog. Oceanogr. 57, 117–173 (2003).Article 

Google Scholar 
10.Vezzulli, L. et al. Climate influence on Vibrio and associated human diseases during the past half-century in the coastal North Atlantic. Proc. Natl. Acad. Sci. USA. 113, E5062–E5071 (2016).CAS 
Article 

Google Scholar 
11.Huq, A. et al. Detection, isolation, and identification of Vibrio cholerae from the environment. Curr. Protoc. Microbiol. https://doi.org/10.1002/9780471729259.mc06a05s26 (2012).12.Thompson, F. L. et al. Phylogeny and molecular identification of vibrios on the basis of multilocus sequence analysis. Appl. Environ. Microb. 71, 5107–5115 (2005).CAS 
Article 

Google Scholar 
13.Vezzulli, L. et al. GbpA as a novel qPCR target for the species-specific detection of Vibrio cholerae O1, O139, Non-O1/Non-O139 in environmental, stool, and historical continuous plankton recorder samples. s. PLoS ONE 10, e0123983 (2015).Article 

Google Scholar 
14.Alam, M. et al. Toxigenic Vibrio cholerae in the aquatic environment of Mathbaria, Bangladesh. Appl. Environ. Microb. 72, 2849–2855 (2006).CAS 
Article 

Google Scholar 
15.Senoh, M. et al. Isolation of viable but nonculturable Vibrio cholerae O1 from environmental water samples in Kolkata, India, in a culturable state. Microbiol. Open 3, 239–246 (2014).CAS 
Article 

Google Scholar 
16.Vezzulli, L. et al. Whole-genome enrichment provides deep insights into vibrio cholerae metagenome from an African river. Microb. Ecol. 73, 734–738 (2017).CAS 
Article 

Google Scholar 
17.Kaboré, S. et al. Occurrence of Vibrio cholerae in water reservoirs of Burkina Faso. Res Microbiol 169, 1–10 (2018).Article 

Google Scholar 
18.Bwire, G. et al. Environmental surveillance of Vibrio cholerae O1/O139 in the five African Great Lakes and other major surface water sources in Uganda. Front. Microbiol. 9, 1560 (2018).Article 

Google Scholar 
19.Vezzulli, L., Baker-Austin, C., Kirschner, A., Pruzzo, C. & Martinez-Urtaza, J. Global emergence of environmental non-O1/O139 Vibrio cholerae infections linked with climate change: a neglected research field? Environ. Microbiol. 22, 4342–4355 (2020).CAS 
Article 

Google Scholar  More

The data needed for estimating the SR relationship consist of spawning biomass (S) and recruitment (R) observed over time. A lognormal distribution is frequently used as the distribution of errors for SR relationships13. We therefore assume that the residuals from a regression model having (r=log (R)) as a response variable and the logarithm of the latent SR relationship as the mean will have a normal distribution. In addition, we assume that the latent SR relationship is likely to be contaminated by some outliers given that fish populations often suffer from nonnegligible contamination, such as sporadic strong cohorts5.Figure 3Parameter estimates of the density-independent parameter (a), density-dependent parameter (b), and autocorrelation ((rho)) for the simulation using the HS SR function with autocorrelation (true (rho = 0.8)) in the residuals.Full size imageFigure 4Application of the robust SR model to fish population data from Japan. (Top) Estimates of ((b-min (S))/(max (S)-min (S))) using the LS and RSR methods. (Bottom) Examples of fitted SR curves using the LS (black line) and RSR (red line) methods (left, walleye pollock in the Sea of Japan; right, round herring in the Tsushima warm current).Full size imageA robust regression approachSuppose that the logarithm of recruitment ((r_t = log (R_t), (t = 1, ldots , T))) has the following autocorrelated normal distribution,$$begin{aligned} r_t = f(S_t|{varvec{theta }})+varepsilon _t, end{aligned}$$
(1)
where (varepsilon _t) is a scaled autoregressive error of order one, that is, (sqrt{lambda _t}(varepsilon _t-rho sqrt{lambda _{t-1}} varepsilon _{t-1})= e_t) with a gaussian noise (e_t) of mean zero and variance (sigma ^2), (S_t) is the spawning biomass, (f(S_t|{varvec{theta }})) is the logarithm of a density-dependent population growth model (spawner-recruitment (SR) curve), ({varvec{theta }}) is the parameter (vector) of the SR curve, (rho) is the autocorrelation, and (sigma ^2) is the base variance of the normal distribution. (lambda _t , (in (0,1])) is the weight for a datum in year t. Rearranging the equation for (varepsilon _t), we have (varepsilon _t sim N(rho sqrt{lambda _{t-1}} varepsilon _{t-1}, sigma ^2/lambda _t)) (Appendix A). We define (lambda _t) to be related to the magnitude of the residual (varepsilon _t),$$begin{aligned} lambda _t = exp left( – phi varepsilon _t^2 right) , end{aligned}$$where (phi , ( >0)) is the parameter that adjusts the influence of outliers. Given that the base variance (sigma ^2) is divided by (lambda _t), the variance is inflated when the difference between the datum and the SR curve is large. The model is equivalent to the AR(1) model when (lambda _t equiv 1) (i.e., (phi =0)) for any t. (sqrt{lambda _t}) is interpreted as the probability of the datum being generated from an uncontaminated normal distribution. When changing the (phi) parameter with (rho =0), the shapes of the probability density function and its derivative are similar to the Tukey’s biweight (also called bisquare) function14, which is close to the gaussian function near zero but decays swiftly as the datum becomes farther from zero (Fig. 1).By solving the equation at equilibrium, the mean deviance residual at (t=1) is zero and the variance at (t=1) is given by ({text{var}}({varepsilon_{1}} ){ = }{sigma ^{2}} {{/}}left[ {lambda _{1}} left( {1} – {rho ^{2}} {tilde{lambda }} right) right]), where ({tilde{lambda }}) is calculated by substituting the sample mean of (lambda _t), (tilde{lambda } = (1/T) sum _{t=1}^T lambda _t) (Appendix B). Incorporating the initial status, the log-likelihood function to be maximized is given by$$begin{aligned} log (L) = sum _{t=1}^T log left( N(r_t|f(S_t|{varvec{theta }})+delta _t, nu _t sigma ^2 lambda _t^{-1}) right) , end{aligned}$$
(2)
where (delta _{t} = 0) and (nu _{t} = (1-rho ^2 tilde{lambda })^{-1}) if (t = 1), and (delta _{t} = rho sqrt{lambda _{t-1}} varepsilon _{t-1}) and (nu _{t} = 1) if (t > 1). Because (varepsilon _{t-1}) increases and (lambda _{t-1}) decreases when there is an outlier at (t-1), the multiplication of (rho) and (sqrt{lambda _{t-1}}) mitigates the influence of an extreme outlier on autocorrelation and contributes to the restoration of the original autocorrelation.We need to estimate the parameters (sigma), (rho), and (phi) in addition to the SR relationship parameters ({varvec{theta }}). The parameter (phi) determines the mixing proportion of contamination and governs the predictive ability of the model. We use time series cross-validation15, which is also called retrospective forecasting16 (RF), to stably determine the value of (phi). First we delete the last datum. Then we use the SR relationship estimated from the data excluding the last datum to forecast recruitment and calculate its error assuming that the deleted recruitment for the last year is true. Next, we delete the two last data, forecast the second-to-last recruitment, and calculate the error assuming that the deleted second-to-last year’s recruitment is true. After the procedure is repeated on a rolling basis, the (phi) parameter having the smallest average error is finally selected. The optimum (phi) is determined by minimizing the following RF error:$$begin{aligned} RF_R = exp left( frac{1}{P} sum _{t=1}^P log left[ left( r_{T-(t-1)} -hat{r}_{T-(t-1)}^{1:(T-t)} right) ^2 right] right) . end{aligned}$$
(3)
This is the geometric mean of predicted errors, which stabilizes the performance of retrospective forecasting. (r_{T-(t-1)}) is the logarithm of observed recruitment in year (T-(t-1)) and (hat{r}_{T-(t-1)}^{1:(T-t)}) is the predicted value estimated using the data from years 1 to (T-t), which is given by$$begin{aligned} hat{r}_{T-(t-1)}^{1:(T-t)} = f(S_{T-(t-1)}|hat{varvec{theta }})+hat{rho } sqrt{hat{lambda }_{T-t}} hat{varepsilon }_{T-t}, end{aligned}$$where (t = 1, ldots , P). We adopt (P=10) for stable estimation in this paper, though we commonly take 5 as the minimum P17.All subsequent analyses are performed using R18 and its package TMB19 (Template Model Builder).SimulationWe generate the simulated data ((left{ (R_t, S_t) ; t = 1, ldots , T right})) with some outliers and autocorrelated errors and test the performance of our robust SR (RSR) method in comparison with the LS and LAD methods. LAD was chosen because it is a typical robust method and is generally superior to the least median squares method used in Chen & Paloheimo (1995)11. The average recruitment data are generated from the Hockey–Stick (HS) SR function12, (f(S_t|{varvec{theta }}) = log left( a min (S_t, b) right)), where ({varvec{theta }} = (a, b) = (1.2, 500)). Stochastic normal errors are added to the log recruitment data with or without autocorrelation. When there is an autocorrelation in the residuals of log recruitment, the autocorrelation is set to (rho = 0.8). To examine the effect of outliers, we add the outliers that occur at the expected frequency of twice per 10 years ((p=0.2)) to the residuals of log recruitment. The patterns of outlier occurrence are threefold: evenly occurring positive and negative outliers ((q=0.5)), all positive outliers ((q=1.0)), and all negative outliers ((q=0.0)) (see Appendix C for the definition of q). We then have eight types of simulated data (no outliers, positive and negative outliers, all positive outliers, and all negative outliers for autocorrelation in the normal residual (rho = 0) and (rho =0.8), respectively). The simulations are replicated 1,000 times for each of the eight types. The length of each SR data time series (T) is set to 30 years which is typical for SR time series data9,12. The performance of the methods is evaluated by two indicators that represent long-term and short-term predictive abilities ((hat{R}_0 – R_0)/R_0) and ((hat{R}_{T+1} – R_{T+1})/R_{T+1}), respectively, where the former is the asymptotic maximum recruitment ((R_0 = ab) for the HS SR function) and the latter is recruitment in the ensuing year (T+1), which is given by (R_{T+1} = exp (f(S_{T+1}|{varvec{theta }}) + rho omega _{T} + eta _{T+1})), where (omega _T) and (eta _{T+1}) are independent gaussian noises (Appendix C). Note that the true recruitment at (T+1) does not include any outliers. The mathematical details of the simulation are given in Appendix C. Autocorrelation is always estimated such that (rho) is set to zero when an estimate of (rho) is equal to or less than zero because a negative autocorrelation is usually impractical20. The parameter (log (phi )) in RSR is chosen from the grid values from (-3.0) to 3.0 in increments of 0.5. The best (phi) is a minimizer of the RF error (RF_R) (Eq. 3).For sensitivity tests, we conduct the following additional simulations: (S1) same as the above base case scenario (S0) except that (a = 1.8); (S2) same as S0 except that (p = 0.1) (the expected frequency of outliers is once every 10 years) in place of (p=0.2); (S3) same as S0 except that (p = 0.3) (the expected frequency of outliers is three times every 10 years) in place of (p=0.2); (S4) same as S0 except that (f(S_t|{varvec{theta }})) is the logarithm of the Beverton–Holt function; (S5) same as S0 except that (f(S_t|{varvec{theta }})) is the logarithm of the Ricker function; S6) same as S0 except for the spawner-abundance dependent p, in which the expected frequency of outliers is higher for lower spawner abundances than for higher spawner abundances.Finally, we calculate biological reference points related to maximum sustainable yield (MSY), i.e., fishing rate at MSY ((F_{rm {msy}})) and spawning biomass at MSY ((S_{rm {msy}})), for each scenario and evaluate their relative biases. To calculate (F_{rm {msy}}) and (S_{rm {msy}}), we require additional information on survival and growth as well as an assumption about population dynamics. For simplicity, we use the delay-difference model as the population dynamics model5. The mathematical details are given in Appendix D.Real data analysisIchinokawa, Okamura & Kurota (2017) fitted the SR curves to fish population data from Japan which comprise 26 SR datasets (Appendix E), demonstrating that some populations showed strong density dependence but others had weak or low density dependence. We fit the HS SR curves to the same 26 SR datasets used in Ichinokawa, Okamura & Kurota (2017). Because Ichinokawa, Okamura & Kurota (2017) used LS as the fitting method, we use LS and RSR to compare the density-independent parameter (log (hat{a})), standardized density-dependent parameter (( hat{b}-min (S) )/( max (S) – min (S) )), autocorrelation in the residuals (hat{rho }), and predictability (hat{RF}_R) in the HS SR curves. More

1.Climate Change: Vital Signs of the Planet (NASA, 2021); https://climate.nasa.gov/vital-signs/carbon-dioxide2.Kimball, B. A. Crop responses to elevated CO2 and interactions with H2O, N, and temperature. Curr. Opin. Plant Biol. 31, 36–43 (2016).CAS 
Article 

Google Scholar 
3.Leakey, A. D. B. et al. Elevated CO2 effects on plant carbon, nitrogen, and water relations: six important lessons from FACE. J. Exp. Bot. 60, 2859–2876 (2009).CAS 
Article 

Google Scholar 
4.Lee, T. D., Barrott, S. H. & Reich, P. B. Photosynthetic responses of 13 grassland species across 11 years of free-air CO2 enrichment is modest, consistent and independent of N supply. Glob. Change Biol. 17, 2893–2904 (2011).Article 

Google Scholar 
5.Ainsworth, E. A. & Long, S. P. 30 years of free-air carbon dioxide enrichment (FACE): what have we learned about future crop productivity and its potential for adaptation? Glob. Change Biol. 27, 27–49 (2021).Article 

Google Scholar 
6.Dusenge, M. E., Duarte, A. G. & Way, D. A. Plant carbon metabolism and climate change: elevated CO2 and temperature impacts on photosynthesis, photorespiration and respiration. New Phytol. 221, 32–49 (2019).CAS 
Article 

Google Scholar 
7.Sardans, J. et al. Ecometabolomics for a better understanding of plant responses and acclimation to abiotic factors linked to global change. Metabolites 10, 239 (2020).CAS 
Article 

Google Scholar 
8.Poorter, H. & Navas, M.-L. Plant growth and competition at elevated CO2: on winners, losers and functional groups. New Phytol. 157, 175–198 (2003).Article 

Google Scholar 
9.Parmesan, C. & Hanley, M. E. Plants and climate change: complexities and surprises. Ann. Bot. 116, 849–864 (2015).Article 

Google Scholar 
10.Ode, P. J., Johnson, S. N. & Moore, B. D. Atmospheric change and induced plant secondary metabolites—are we reshaping the building blocks of multi-trophic interactions? Curr. Opin. Insect Sci. 5, 57–65 (2014).Article 

Google Scholar 
11.Robinson, E. A., Ryan, G. D. & Newman, J. A. A meta-analytical review of the effects of elevated CO2 on plant–arthropod interactions highlights the importance of interacting environmental and biological variables. New Phytol. 194, 321–336 (2012).CAS 
Article 

Google Scholar 
12.Willeit, M., Ganopolski, A., Calov, R. & Brovkin, V. Mid-Pleistocene transition in glacial cycles explained by declining CO2 and regolith removal. Sci. Adv. 5, eaav7337 (2019).CAS 
Article 

Google Scholar 
13.Busch, F. A. & Sage, R. F. The sensitivity of photosynthesis to O2 and CO2 concentration identifies strong Rubisco control above the thermal optimum. New Phytol. 213, 1036–1051 (2017).CAS 
Article 

Google Scholar 
14.Drake, B. G., Gonzàlez-Meler, M. A. & Long, S. P. More efficient plants: a consequence of rising atmospheric CO2? Annu. Rev. Plant Physiol. Plant Mol. Biol. 48, 609–639 (1997).CAS 
Article 

Google Scholar 
15.Ziska, L. H., Sicher, R. C., George, K. & Mohan, J. E. Rising atmospheric carbon dioxide and potential impacts on the growth and toxicity of poison ivy (Toxicodendron radicans). Weed Sci. 55, 288–292 (2007).CAS 
Article 

Google Scholar 
16.Ziska, L. H. & Caulfield, F. A. Rising CO2 and pollen production of common ragweed (Ambrosia artemisiifolia L.), a known allergy-inducing species: implications for public health. Funct. Plant Biol. 27, 893 (2000).Article 

Google Scholar 
17.Ziska, L. H., Panicker, S. & Wojno, H. L. Recent and projected increases in atmospheric carbon dioxide and the potential impacts on growth and alkaloid production in wild poppy (Papaver setigerum DC.). Clim. Change 91, 395 (2008).CAS 
Article 

Google Scholar 
18.Del Fabbro, C. & Prati, D. The relative importance of immediate allelopathy and allelopathic legacy in invasive plant species. Basic Appl. Ecol. 16, 28–35 (2015).Article 

Google Scholar 
19.Ni, G. et al. Exploring the novel weapons hypothesis with invasive plant species in China. Allelopath. J. 29, 199–214 (2012).
Google Scholar 
20.Peñuelas, J. et al. Higher allocation to low cost chemical defenses in Iinvasive species of Hawaii. J. Chem. Ecol. 36, 1255–1270 (2010).Article 

Google Scholar 
21.Bajwa, A. A., McClay, A. & Adkins, S. W. in Parthenium Weed: Biology, Ecology and Management (eds Adkins, S., Shabbir, A. et al.) 7–39 (CABI, 2019).22.Adkins, S. & Shabbir, A. Biology, ecology and management of the invasive parthenium weed (Parthenium hysterophorus L.): management of parthenium weed. Pest Manag. Sci. 70, 1023–1029 (2014).CAS 
Article 

Google Scholar 
23.Niranjan, A. et al. Identification and quantification of heterologous compounds parthenin and organic acids in Parthenium hysterophorus L. using HPLC-PDA-MS-MS. Anal. Lett. 46, 48–59 (2013).Article 

Google Scholar 
24.Belz, R. G., van der Laan, M., Reinhardt, C. F. & Hurle, K. Soil degradation of parthenin—does it contradict the role of allelopathy in the invasive weed Parthenium hysterophorus L.? J. Chem. Ecol. 35, 1137–1150 (2009).CAS 
Article 

Google Scholar 
25.Hanif, Z., Adkins, S. W., Prentis, P. J., Navie, S. C. & O’Donnell, C. J. Characterization of the reproductive behaviour and invasive potential of parthenium weed in Australia. Pak. J. Weed Sci. Res. 18, 767–774 (2012).
Google Scholar 
26.Bajwa, A. A., Chauhan, B. S. & Adkins, S. Morphological, physiological and biochemical responses of two Australian biotypes of Parthenium hysterophorus to different soil moisture regimes. Environ. Sci. Pollut. Res. 24, 16186–16194 (2017).CAS 
Article 

Google Scholar 
27.Nguyen, T., Bajwa, A. A., Navie, S., O’Donnell, C. & Adkins, S. Parthenium weed (Parthenium hysterophorus L.) and climate change: the effect of CO2 concentration, temperature, and water deficit on growth and reproduction of two biotypes. Environ. Sci. Pollut. Res. 24, 10727–10739 (2017).Article 

Google Scholar 
28.Chadwick, M., Trewin, H., Gawthrop, F. & Wagstaff, C. Sesquiterpenoids lactones: benefits to plants and people. Int. J. Mol. Sci. 14, 12780–12805 (2013).Article 

Google Scholar 
29.Ojija, F., Arnold, S. E. J. & Treydte, A. C. Impacts of alien invasive Parthenium hysterophorus on flower visitation by insects to co-flowering plants. Arthropod Plant Interact. 13, 719–734 (2019).Article 

Google Scholar 
30.Bajwa, A. A., Chauhan, B. S. & Adkins, S. W. Germination ecology of two Australian biotypes of ragweed parthenium (Parthenium hysterophorus) relates to their invasiveness. Weed Sci. 66, 62–70 (2018).Article 

Google Scholar 
31.Bajwa, A. A. et al. Toxic potential and metabolic profiling of two Australian biotypes of the invasive plant parthenium weed (Parthenium hysterophorus L.). Toxins 12, 447 (2020).CAS 
Article 

Google Scholar 
32.Grime, J. P. Plant Strategies, Vegetation Processes, and Ecosystem Properties (Wiley, 2001).33.Grime, J. P. in Plant Evolutionary Biology (eds Gottlieb, L. D. & Jain, S. K.) 371–393 (Springer, 1988).34.Craine, J. M. Reconciling plant strategy theories of Grime and Tilman. J. Ecol. 93, 1041–1052 (2005).Article 

Google Scholar 
35.Bae, J. et al. Effect of elevated atmospheric carbon dioxide on the allelopathic potential of common ragweed. J. Ecol. Environ. 43, 21 (2019).Article 

Google Scholar 
36.Wang, R.-L. et al. Responses of Mikania micrantha, an invasive weed to elevated CO2: induction of β-caryophyllene synthase, changes in emission capability and allelopathic potential of β-caryophyllene. J. Chem. Ecol. 36, 1076–1082 (2010).CAS 
Article 

Google Scholar 
37.Robinson, J. M. Photosynthetic carbon metabolism in leaves and isolated chloroplasts from spinach plants grown under short and intermediate photosynthetic periods. Plant Physiol. 75, 397–409 (1984).CAS 
Article 

Google Scholar 
38.R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2019).39.Filion, M., Dutilleul, P. & Potvin, C. Optimum experimental design for Free-Air Carbon dioxide Enrichment (FACE) studies. Glob. Change Biol. 6, 843–854 (2000).Article 

Google Scholar 
40.Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. https://www.jstatsoft.org/article/view/v067i01 (2015).41.Searle, S. R., Speed, F. M. & Milliken, G. A. Population marginal means in the linear model: an alternative to least squares means. Am. Stat. 34, 216–221 (1980).
Google Scholar  More

Our findings illustrate both the range and extent of wildlife exploitation in Wuhan markets, prior to new trading bans linked to the COVID-19 outbreak, along with the poor conditions under which these animals were kept prior to sale. Circumstantially, the absence of pangolins (and bats, not typically eaten in Central China; media footage generally depicts Indonesia) from our comprehensive survey data corroborates that pangolins are unlikely implicated as spill-over hosts in the COVID-19 outbreak. This is unsurprising because live pangolin trading has largely ceased in China13.We should therefore not be complacent, because the original source of COVID-19 does not seem to have been established. This is doubly important because false attribution can lead to extreme and irresponsible animal persecution. For instance, civets were killed en masse following the SARS-CoV outbreak5, and any unwarranted vilification or persecution of pangolins and bats in relation to COVID-19 would risk undermining otherwise very successful efforts to better protect and conserve wildlife in China.Regarding our insights into broader IWT issues in Wuhan, the animals sold were relatively expensive, representing luxury food items, not cheap bushmeat (Table 1). We thus make an ethical distinction here between the subsistence consumption of bush meat in poorer nations, versus the sort of cachet attached to wild animal consumption in parts of the developed world, notably China14, but also Japan15. While c. 30% of mammals were clearly wild-caught, indicated by trapping and shooting wounds, the captive breeding of other species is commonplace in China. Raccoon dog fur farming is legal in China; however, due to a drop in fur prices, raccoon dogs are now frequently sold off in live animal markets, augmented by wild-caught individuals. Similarly, all American mink (Neovison vison) originated from fur farms—noting that SARS-CoV-2 has been reported in mink farms in Europe and North America16, 17. In contrast, the captive breeding and sale of Siberian weasels (Mustela sibirica), is totally illegal in China, yet they are easy to breed, and sold openly, without attracting law enforcement. Indeed, prior to COVID-19 reforms, although enforcement officers from the Wuhan Forestry Bureau issued permits to market vendors, they were broadly disinterested in what species were sold. Furthermore, although animals were required to have an origin certificate and be quarantined to ensure they did not exhibit overt disease symptoms, no clear policy was enforced on these conditions. This is important because the species that were traded are capable of hosting a wide range of infectious zoonotic diseases or disease-baring parasites (Supplementary Table S1), aside from COVID-19. These range from potentially lethal viruses, for example, rabies, SFTS, H5N1, to common bacterial infections that, nevertheless, represent a risk to human health (e.g., Streptococcus). Indeed, globally, wildlife is thought to be the source of at least 70% of all emerging diseases18.Legislative reform is also vital to clarify unequivocally which species are considered ‘wild’ and cannot be traded legally and safely. Another problem, as encountered by the WHO report is that, retrospectively, it proved difficult to ascertain which species were on sale, even to the genus level, relying solely on the responsible market authority’s official sales records and disclosures1. As we19, 20, and others21, have proposed previously, China’s LFSSP and LESS must be updated to apply proper binomials, and to align with recent taxonomic revisions; for instance, cobra snakes (Nada atra) can be farmed legally for food with permits, but wild caught species, such as water snakes and wolf snakes were also sold in Wuhan, labelled simply as ‘snakes’. Such an application of clear species names would allow for more effective prosecutions19. Furthermore, the WHO reports that market authorities claimed all live and frozen animals sold in the Huanan market were acquired from farms officially licensed for breeding and quarantine, and as such no illegal wildlife trade was identified1. In reality, however, because China has no regulatory authority regulating animal trading conducted by small-scale vendors or individuals it is impossible to make this determination1, 21. Similar discrepancies concerning species identification and origins afflict investigations around the world22.Another important animal trade that requires attention, outside of exploitation as food, is the supply of pets, like the squirrels and crested myna birds sold in Wuhan’s market. Our previous research found annual trade volumes equivalent to c. 17,000 parrots and c. 160,000 turtles (many turtles being invasive if escaping to the wild) sold online as pets via Taobao.com between 2016–2017, in contravention of China’s WACL and/or the Animal Epidemic Prevention Law23,24,25. While not currently the vector of any major viral epidemics, it would be naive to imagine that unconventional pets do not still also pose a serious concern for public health26. This potential for disease is likely exacerbated by poor sanitary and welfare conditions (Fig. 2). More

Sample collection and incubationThree replicates of soil samples were collected from the top 10 cm in of plant-free patches in four ecosystems along the C. Hart Merriam elevation gradient in Northern Arizona25 beginning at high desert grassland (1760 m), and followed at higher elevations by piñon-pine juniper woodland (2020 m), ponderosa pine forest (2344 m), and mixed conifer forest (2620 m). Soils were air-dried for 24 h at room temperature, homogenized, and passed through a 2 mm sieve before being stored at 4 °C for another 24 h. Soil incubations were performed on soils with mass of 20 g of dry soil for measurements of CO2 and microbial biomass carbon (MBC), while 2 g of dry soil aliquots were incubated separately (but under equivalent conditions) for quantitative stable isotope probing (qSIP). We applied three treatments to these soils through the addition of water (up to 70% water-holding capacity): water alone (control), with glucose (C treatment; 1000 µg C g−1 dry soil), or with glucose and nitrogen (C + N treatment; [NH4]2SO4 at 100 µg N g−1 dry soil). All samples for qSIP were incubated with 18O-enriched water (97 atom%) and matching controls necessary to calculate the change in 18O enrichment across the microbial community. We applied water at natural abundance (i.e., no 18O-enriched water) to the larger soil samples prepared for measurement of carbon flux. All soils were incubated in the dark for one week. Following incubation, soils were frozen at −80 °C for 1 week prior to DNA extraction.Soil, CO2, and microbial biomass measurementsWe analyzed headspace gas of soils for CO2 concentration and δ13CO2 three times during the week-long incubation using a LI-Cor 6262 (LI-Cor Biosciences Inc. Lincoln, NE, USA) and a Picarro G2201 (Picarro Inc., Sunnyvale, CA, USA), respectively. Prior to incubation we analyzed soil MBC using the chloroform-fumigation extraction method on 10 g of soil. One sub-sample was immediately extracted with 25 ml of a 0.05 M K2SO4 solution, while a second sub-sample was first fumigated with chloroform (for 5 days), after which it was similarly extracted. Following K2SO4 addition, we agitated soils for 1 h, filtered the extract through a Whatman #3 filter paper, and dried the filtered solution (60 °C, 4 days). Salts with extracted C were ground and analyzed for total C using an elemental analyzer coupled to a mass spectrometer. MBC was calculated as the difference between the fumigated and immediately extracted samples’ soil C using an extraction efficiency of 0.45 (as per Liu et al.26).Quantitative stable isotope probingWe performed DNA extraction and 16S amplicon sequencing on 18O-incubated qSIP soils11,12,13. The procedures targeted the V4 region of the 16S gene as specified by the Earth Microbiome Project (EMP, http://www.earthmicrobiome.org) standard protocols27,28. We used PowerSoil DNA extraction kits following manufacture instructions to isolate DNA from soil (MoBio laboratories, Carlsbad, CA, USA). We quantified extracted DNA using the Qubit dsDNA High-Sensitivity assay kit and a Qubit 2.0 Fluorometer (Invitrogen, Eugene, OR, USA). To quantify the degree of 18O isotope incorporation into bacterial DNA, we performed density fractionation and sequenced 15–18 fractions separately following methods modified from the canonical publication7. We added 1 µg of DNA to 2.6 mL of saturated CsCl solution in combination with a gradient buffer (200 mM Tris, 200 mM KCL, 2 mM EDTA) in a 3.3 mL OptiSeal ultracentrifuge tube (Beckman Coulter, Fullerton, CA, USA). The solution was centrifuged to produce a gradient of increasingly labeled (heavier) DNA in an Optima Max bench top ultracentrifuge (Beckman Coulter, Brea, CA, USA) with a Beckman TLN-100 rotor (127,000 × g for 72 h) at 18 °C. We separated each sample from the continuous gradient into approximately 20 fractions (150 µL) using a modified fraction recovery system (Beckman Coulter). We then measured the density of each separate fraction with a Reichart AR200 digital refractometer (Reichert Analytical Instruments, Depew, NY, USA) and retained fractions with densities between 1.640 and 1.735 g cm−3. We cleaned and purified DNA in these fractions using isopropanol precipitation, quantified DNA using the Quant-IT PicoGreen dsDNA assay (Invitrogen) and a BioTek Synergy HT plate reader (BioTek Instruments Inc., Winooski, VT, USA), and quantified bacterial 16S gene copies using qPCR (primers: Supplementary Table 1) in triplicate. We used 8 µL reactions consisting of 0.2 mM of each primer, 0.01 U µL−1 Phusion HotStart II Polymerase (Thermo Fisher Scientific, Waltham, MA), 1× Phusion HF buffer (Thermo Fisher Scientific), 3.0 mM MgCl2, 6% glycerol, and 200 µL of dNTPs. We amplified DNA using a Bio-Rad CFX384 Touch real-time PCR detection system (Bio-Rad, Hercules, CA, USA) with the following cycling conditions: 95 °C at 1 min and 44 cycles of 95 °C (30 s), 64.5 °C (30 s), and 72 °C (1 min).We sequenced the 16S V4 region (primers: EMP standard 515F—806R; Supplementary Table 1) on an Illumina MiSeq (Illumina, Inc., San Diego, CA, USA). Sequences were amplified using the same reaction mix as qPCR amplification but cycling at 95 °C for 2 min followed by 15 cycles of 95 °C (30 s), 55 °C (30 s), and 60 °C (4 min). In addition to post-incubation soils, we extracted, amplified, and sequenced DNA of the bacterial community at the start of the incubation.Sequence processing and qSIP analysisThe raw sequence data of forward and reverse reads (FASTQ) were processed within the QIIME 2 environment (release 2018.6)29,30, denoising sequences with the available DADA2 pipeline31. We clustered the remaining sequences into amplicon sequence variants or ASVs (at 100% sequence identity) against the SILVA 132 database32 using an open-reference Naïve Bayes feature classifier33. We removed global singletons and doubleton ASVs, non-bacterial lineages, and samples with less than 4000 sequence reads. Removal of global singletons and doubletons resulted in the removal of 2241 unique ASVs from the feature table yielding 115,647 out of 117,888 (a retention of 98% of all ASVs) as well as the loss of 4018 sequences leaving 37,765,678 (a retention >99% of all sequences). We combined taxonomic information and ASV sequence counts with per-fraction qPCR and density measurements using the phyloseq package (version 1.24.2), in R (version 3.5.1)34. Because high-throughput sequencing produces relativized measures of abundance, we converted ASV sequencing abundances in each fraction to the number of 16S rRNA gene copies per g dry soil based on the known amount of dry soil added and the amount of DNA in each soil sample. All data and analytical code have been made publicly accessible35.To perform qSIP analysis and calculate per-capita growth rates of each ASV, we used our in-house qsip package (https://github.com/bramstone/qsip) based on previously published research7,10. Because rare and infrequent taxa are more likely to be lost in samples with poor sequencing depth with their absences affecting DNA density changes, we invoked a presence or absence-based filtering criteria on ASVs prior to calculation of per-capita growth rates. Within each ecosystem, we kept only ASVs that appeared in two of the three replicates of a treatment (18O, C, and C + N) and at that appeared in at least five of the fractions within each of those two replicates. ASVs filtered out of one treatment were allowed to appear in another if they met the frequency threshold.For all remaining ASVs (1081 representing less than 1% of all ASVs but 58% of all sequence reads), we calculated per-capita gross growth (i.e., cell division) rates observed in each replicate using an exponential growth model10. We applied these per-capita rates to the number of 16S rRNA gene copies to estimate the production of new 16S rRNA gene copies of each ASV per g dry soil per week using the following equation:$$frac{{rm{d}}{N}_{{rm{i}}}}{{{rm{d}}t}}={N}_{{rm{i,t}}}-{N}_{{rm{i,t}}}{e}^{-{g}_{{rm{i}}}t},$$
(1)
Where Ni,t is the number of 16S rRNA gene copies of taxon i at time t (here after 7 days) and gi represents the per-capita growth rate (calculated as a daily rate). See Supplementary Fig. 3 for results on the production of 16S gene copies.Calculation of 16S rRNA gene copy numbers and cell massIn parallel to taxonomic assignment, we compared quality-filtered 16S sequences against a database of 12,415 complete prokaryote genomes obtained from GenBank. From these genomes, we extracted data on 16S rRNA gene copy number, total genome size, and 16S gene sequence. We used BLAST to find matches against this database to the ASVs generated from QIIME 2 to make per-taxon assignments of 16S rRNA gene copy number and total genome size13. For ASVs that did not find an exact match, we assigned 16S rRNA gene copy number values and genome sizes based on the median values observed in the most specific possible taxonomic rank. We estimated the mass of individual cells for each population using published allometric scaling relationships between genome length and cellular mass from West and Brown:36$${{{log }}}_{10}({M}_{{rm{i}}})=frac{{{{log }}}_{10}left({G}_{{rm{i}}}right)-9.4}{0.24},$$
(2)
where Mi indicates cellular mass (g) and Gi indicates genome length (bp) for taxon i. We obtained this relationship by digitizing Fig. 436 using DataThief III and re-fitting the trend line in log–log space. We estimated that 20% of the cellular mass was carbon37. To validate this approach, cellular mass estimates and initial 16S copy number measurements were used to estimate population-level biomass C values which were summed and compared to initial community-level MBC. We found that these values overestimated initial MBC by an order of magnitude. As such, cellular carbon mass was divided by 10 in our final calculations. We applied cellular mass and 16S copy number estimates to the production of 16S copies to estimate the production of biomass carbon for each taxon during the incubation period (t):$${P}_{{rm{i}}}=frac{{rm{d}}{N}_{{rm{i}}}/{{rm{d}}t}}{C_{{rm{i}}}}cdot {M}_{{rm{i}}}cdot 0.2,$$
(3)
where Pi indicates production of biomass carbon (µg C g dry soil−1 week−1) and Ci indicates 16S copy number per cell for taxon i. The 0.2 coefficient represents an estimate that 20% of cellular mass is composed of carbon.Efficiency and respiration modelingWe estimated rates of respiration using qSIP-informed growth rates and community-level carbon use efficiency (CUE). CUE estimates were based on the incorporation of 18O-water into DNA as a measure of gross biomass production38,39 and measured CO2 in headspace gas from soil incubations. We calculated the production of 18O-labeled biomass carbon (18P) at the community-level for each sample by summing the products of per-taxon 18O enrichment (excess atom fraction, EAF) and relative abundance:$${, }^{18}{P}=mathop{sum }limits_{i=1}^{n}({,}^{18}{{{rm{EAF}}}}_{{rm{i}}}cdot {y}_{{rm{i}}})cdot {rm{DN}}{rm{A}}_{0}cdot fleft({{rm{MB}}}{rm{C}}_{0} sim {rm{DN}}{rm{A}}_{0}right),$$
(4)
where 18P indicates the gross production of 18O-labeled microbial biomass carbon per gram of dry soil per week, 18EAFi indicates the enrichment of DNA of taxon i and yi indicates its relative abundance, DNA0 indicates the concentration of DNA per gram of dry soil prior to incubation, and MBC0 indicates the microbial biomass carbon per gram of dry soil prior to incubation. Here, the MBC0 ~ DNA0 function indicates the linear relationship between MBC and DNA concentration. We used the output from Eq. 4 to calculate community CUE for each sample:$${{rm{CUE}}}=frac{{,}^{18}{{P}}}{(!{,}^{18}P+R)},$$
(5)
where R indicates the total CO2 respired per gram dry soil per week.We used the community CUE values from each sample (Eq. 5) to constrain/as upper and lower limits our estimates of per-taxon CUE. For a group of three replicates from a given ecosystem and treatment, we used the minimum and maximum observed community-level CUE values as the acceptable range of per-taxon CUE values. These constraints were used to control the shape of the function of per-taxon CUE and growth rate, though functions were modeled both with and without constraints (i.e., per-taxon CUE values were bounded only by 0 and 0.7). The range of community-level CUE values for each treatment were 0.18–0.53 for control soils, 0.04–0.13 for carbon amended soils and 0.03–0.08 for carbon and nitrogen amended soils and did not vary much between ecosystems. As a result of uncertainty in the literature about the relationship between growth rate and CUE14, several different relationships were postulated to model per-taxon CUE as a function of per-taxon growth rate: linear increase, linear decrease, exponential decrease, unimodal with peak CUE at growth rate of 0.5, and unimodal with peak CUE at a growth rate of 0.05 (the median of all per-taxon growth rates in the data). Comparisons between functions were made by calculating AIC values from per-taxon respiration, summed, and regressing against measured respiration values. Likewise, for each function, we tested how well per-taxon CUE estimates reconstructed community-level CUE by weighting the CUE value of each taxon by its relative abundance, summing, and regressing against community-level CUE. To select the best per-taxon CUE function, AIC values from both scaling efforts were combined. To make AIC values comparable, all respiration and CUE terms were z-transformed prior to regression scaling. To reflect our priority of estimating per-taxon respiration, AIC values from the respiration scaling regression models were multiplied by two and summed with AIC values from CUE scaling such that AICTotal = 2(AICResp) + AICCUE. Across these comparisons, the best estimate of per-taxon CUE was the unimodal function of growth rate, constrained by community-level CUE and peaking at growth rates of 0.5 (Table 1), such that:$${{rm{CUE}}}_{{rm{i}}}=-4({{rm{CUE}}}_{{rm{E}}{rm{:}}{rm{T}}{rm{:}}{{rm{range}}}})cdot {left({g}_{{rm{i}}}-0.5right)}^{2}+({{rm{CUE}}}_{{rm{E}}{rm{:}}{rm{T}}{rm{:}}{max }}),$$
(6)
where CUEi indicates per-taxon CUE, CUEE:T:max indicates the maximum CUE values observed for a group of replicates within a given ecosystem and treatment (E:T). With this function, higher per-capita growth rate values were parameterized to produce higher CUE values initially and then decrease reflecting a growth-CUE tradeoff14, here bound by the difference in maximum and minimum CUE values. We applied per-taxon CUE estimates from Eq. 6 to per-taxon growth rates to yield estimates of per-taxon respiration:$${r}_{{rm{i}}}={r}_{{rm{g,i}}}+{r}_{{rm{m,i}}}=left(frac{{g}_{{rm{i}}}}{{{rm{CUE}}}_{{rm{i}}}}-{g}_{{rm{i}}}right)+left(frac{{g}_{{rm{i}}}}{{{rm{CUE}}}_{{rm{i}}}}-{g}_{{rm{i}}}right)cdot beta,$$
(7)
where ri indicates per-capita respiration for taxon i, rg,i indicates growth-related respiration, rm,i indicates maintenance-related respiration, and β is a constant of 0.01 that represents the maintenance requirements as a proportion of total energy use40. We used these values of per-taxon, per-capita respiration rates to estimate per-taxon respiration per gram of dry soil per week:$${R}_{{rm{i}}}={P}_{{rm{i}}}cdot {r}_{{{rm{g,i}}}}+{P}_{{rm{i}}}cdot {r}_{{{rm{m,i}}}},$$
(8)
where Ri indicates respiration of CO2–C (µg C g dry soil−1 week−1) for taxon i.In addition to per-taxon respiration estimates based on 18O enrichment, we used another model for comparison. Here, respiration was calculated based on 16S abundance alone:$${R}_{{rm{i}}}={N}_{{rm{i}}}cdot f(R sim N+0),$$
(9)
where Ni indicates final 16S abundance for taxon i, R indicates microbial respiration of CO2-C (µg C g dry soil−1 week−1) and N indicates total 16S abundance at the end of the incubation. Here, the R ~ N function indicates the linear relationship, with an intercept of 0, between CO2 respiration and 16S gene concentration across all samples.Diversity, compositional, and statistical analysisFor patterns of evenness in bacterial carbon use and relative abundance, we used Pielou’s evenness which is the quotient of Shannon’s diversity and the observed richness. For each sample, we applied Pielou’s evenness to bacterial abundances as well as bacterial carbon use (relativized to sum to one, in both cases).We created a linear mixed model to test the relationship between the carbon use (the sum of biomass production and respiration) and relative abundance of bacterial genera from the dominant phyla, which accounted for >90% of all C flux. Here, we averaged carbon use and relative abundance for all replicates in a given ecosystem and treatment. We used the lme4 R package (version 1.1-20)41 and obtained p-values using the Satterthwaite method in the lmerTest R package (version 3.1-0)42. To limit pseudo-replication, we accounted for differences in carbon use across ecosystems and due to bacterial Genus by implementing random intercepts. We selected for the optimal random and fixed components by dropping individual terms and comparing models with likelihood ratio tests, disregarding models that failed to converge. Our final model fit was:$${{{log }}}_{10}({C}_{{rm{i}}}) sim {{{log }}}_{10}left({y}_{{rm{i}}}right)ast T+left(1|Eright)+(1|{{rm{Genus}}}),$$
(10)
where Ci indicates the relativized carbon use for taxon i (averaged across all three replicates in a given ecosystem and treatment), yi indicates the relative abundance of taxon i (averaged across all three replicates), T indicates soil treatment, and E indicates ecosystem.For differences in composition, we created species abundance tables using the traditional abundances, as well as measures of carbon use (growth and maintenance respiration) of each ASV in each sample. To account for differences in absolute abundances and flux rates between sites, we relativized all abundance tables. We summarized compositional differences using Bray–Curtis dissimilarities then identified multivariate centroids for all replicates in a site by treatment group. We tested the effect of site and nutrient amendment on the resulting group centroids using PERMANOVA tests implemented with the adonis function in the vegan package (version 2.5-3)43. We related compositional shifts in relative abundance to those in relativized growth and maintenance using Mantel tests with the mantel function in vegan.To test for changes in the type of soil C preferred by microbial genera (either 13C-labeled glucose or 12C soil carbon) in response to nitrogen addition, we used Levene’s test with the car package (version 3.0-10)44. Specifically, we analyzed the relationship between 13C use and 12C use (both relativized) on bacterial genera across all replicates and in C and C + N treatments using a linear model. We then extracted model residuals and tested whether variance was significantly different across treatments by focusing on the interaction between individual replicates and treatment. This produced a significance test describing treatment-level differences in 13C–12C use.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More