More stories

  • in

    Resistance status of lepidopteran soybean pests following large-scale use of MON 87701 × MON 89788 soybean in Brazil

    1.Macrae, T. C. et al. Laboratory and field evaluations of transgenic soybean exhibiting high-dose expression of a synthetic Bacillus thuringiensis cry1A gene for control of Lepidoptera. J. Econ. Entomol. 98, 577–587 (2005).Article 

    Google Scholar 
    2.Miklos, J. A. et al. Characterization of soybean exhibiting high expression of a synthetic Bacillus thuringiensis cry1A transgene that confers a high degree of resistance to lepidopteran pests. Crop Sci. 47, 148–157 (2007).CAS 
    Article 

    Google Scholar 
    3.Spark. BIP soja 2021. Spark Smarter Decisions, Valinhos, SP. http://spark-ie.com.br (2021).4.Bernardi, O. et al. Assessment of the high dose concept and level of control provided by MON 87701 × MON 89788 soybean against Anticarsia gemmatalis and Pseudoplusia includens (Lepidoptera: Noctuidae) in Brazil. Pest Manag. Sci. 68, 1083–1091 (2012).CAS 
    Article 

    Google Scholar 
    5.Bernardi, O. et al. High levels of biological activity of Cry1Ac protein expressed on MON 87701 × MON 89788 soybean against Heliothis virescens (Lepidoptera: Noctuidae). Pest Manag. Sci. 70, 588–594 (2014).CAS 
    Article 

    Google Scholar 
    6.Dourado, P. M. et al. High susceptibility to Cry1Ac and low resistance allele frequency reduce the risk of resistance of Helicoverpa armigera to Bt soybean in Brazil. PLoS ONE 11, e0161388 (2016).Article 

    Google Scholar 
    7.Horikoshi, R. J. et al. Large-scale assessment of lepidopteran soybean pests and efficacy of Cry1Ac soybean in Brazil. Sci. Rep. 11, 15956 (2021).ADS 
    CAS 
    Article 

    Google Scholar 
    8.Gould, F. Sustainability of transgenic insecticidal cultivars: integrating pest genetics and ecology. Annu. Rev. Entomol. 43, 701–726 (1998).CAS 
    Article 

    Google Scholar 
    9.Tabashnik, B. E., Van Rensburg, J. B. J. & Carrière, Y. Field-evolved insect resistance to Bt crops: Definition, theory and data. J. Econ. Entomol. 102, 2011–2025 (2009).CAS 
    Article 

    Google Scholar 
    10.Yano, S. A. et al. High susceptibility and low resistance allele frequency of Chrysodeixis includens (Lepidoptera: Noctuidae) field populations to Cry1Ac in Brazil. Pest Manag. Sci. 72, 1578–1584 (2015).Article 

    Google Scholar 
    11.Tabashnik, B. E. & Carrière, Y. Surge in insect resistance to transgenic crops and prospects for sustainability. Nat. Biotechnol. 35, 926–935 (2017).CAS 
    Article 

    Google Scholar 
    12.CTNBio– Comissão Técnica Nacional de Biossegurança. Technical Opinion No. 2542/2010-Commercial release of genetically modified insect-resistant and herbicide-tolerant soy containing genetically modified events MON 87701 and MON 89788. http://ctnbio.mctic.gov.br/en/liberacao-comercial#/liberacao-comercial/consultar-processo (2010).13.Van Rensburg, J. B. J. First report of field resistance by stem borer, Busseola fusca (Fuller) to Bt-transgenic maize. S. Afr. J. Plant Soil 24, 147–151 (2007).Article 

    Google Scholar 
    14.Storer, N. P. et al. Discovery and characterization of field resistance to Bt maize: Spodoptera frugiperda (Lepidoptera: Noctuidae) in Puerto Rico. J. Econ. Entomol. 103, 1031–1038 (2010).Article 

    Google Scholar 
    15.Dhurua, S. & Gujar, G. T. Field-evolved resistance to Bt toxin Cry1Ac in the pink bollworm, Pectinophora gossypiella (Saunders) (Lepidoptera: Gelechiidae), from India. Pest Manag. Sci. 67, 898–903 (2011).CAS 
    Article 

    Google Scholar 
    16.Gassmann, A. J., Petzold-Maxwell, J. L., Keweshan, R. S. & Dunbar, M. W. Field-evolved resistance to Bt maize by western corn rootworm. PLoS ONE 6, e22629 (2011).ADS 
    CAS 
    Article 

    Google Scholar 
    17.Farias, J. R. et al. Field-evolved resistance to Cry1F maize by Spodoptera frugiperda (Lepidoptera: Noctuidae) in Brazil. Crop Prot. 64, 150–158 (2014).ADS 
    Article 

    Google Scholar 
    18.Bernardi, D. et al. Cross-resistance between Cry1 proteins in fall armyworm (Spodoptera frugiperda) may affect the durability of current pyramided Bt maize hybrids in Brazil. PLoS ONE 10, e0140130 (2015).Article 

    Google Scholar 
    19.Mohan, K. S., Ravi, K. C., Suresh, P. J., Sumerford, D. & Head, G. P. Field resistance to the Bacillus thuringiensis protein Cry1Ac expressed in Bollgard hybrid cotton in pink bollworm, Pectinophora gossypiella (Saunders), populations in India. Pest Manag. Sci. 72, 738–746 (2016).CAS 
    Article 

    Google Scholar 
    20.Omoto, C. O. et al. Field-evolved resistance to Cry1Ab maize by Spodoptera frugiperda in Brazil. Pest Manag. Sci. 72, 1727–1736 (2016).CAS 
    Article 

    Google Scholar 
    21.Naik, V. C., Kumbhare, S., Kranthi, S., Satija, U. & Kranthi, K. R. Field-evolved resistance of pink bollworm, Pectinophora gossypiella (Saunders) (Lepidoptera: Gelechiidae), to transgenic Bacillus thuringiensis (Bt) cotton expressing crystal 1Ac (Cry1Ac) and Cry2Ab in India. Pest Manag. Sci. 74, 2544–2554 (2018).CAS 
    Article 

    Google Scholar 
    22.Sanchez, N. E. & Pereyra, P. C. Life tables of the soybean looper Rachiplusia nu (Lepidoptera: Noctuidae) in the laboratory. Rev. Soc. Entomol. Arg. 54, 89–96 (1995).
    Google Scholar 
    23.Sánchez, N. E., Pereyra, P. C. & Gentile, M. V. Population parameters of Epinotia aporema (Lepidoptera: Tortricidae) on soybean. Rev. Soc. Entomol. Arg. 56, 151–153 (1997).
    Google Scholar 
    24.Barrionuevo, M. J., Murúa, M. G., Goane, L., Meagher, R. & Navarro, F. Life table studies of Rachiplusia nu (Guenée) and Chrysodeixis (= Pseudoplusia) includens (Walker) (Lepidoptera: Noctuidae) on artificial diet. Fla. Entomol. 95, 944–951 (2012).Article 

    Google Scholar 
    25.Ferreira, B. S. C. Sampling Epinotia aporema on Soybean in Sampling Methods in Soybean Entomology (eds. Kogan, M. & Herzog, D. C) 374–381 (Springer 1980).26.Specht, A. et al. Biotic potential and life tables of Chrysodeixis includens (Lepidoptera: Noctuidae), Rachiplusia nu, and Trichoplusia ni on soybean and forage turnip. J. Insect Sci. 19, 1–8 (2019).Article 

    Google Scholar 
    27.Santos, S. R. D., Specht, A., Carneiro, E. & Casagrande, M. M. The influence of agricultural occupation and climate on the spatial distribution of Plusiinae (Lepidoptera: Noctuidae) on a latitudinal gradient in Brazil. Rev. Bras. Entomol. 65, e20200103 (2021).Article 

    Google Scholar 
    28.Greene, G. L., Leppla, N. C. & Dickerson, W. A. Velvetbean caterpillar: A rearing procedure and artificial medium. J. Econ. Entomol. 69, 487–488 (1976).Article 

    Google Scholar 
    29.Andow, D. A. & Alstad, D. N. F2 screen for rare resistance alleles. J. Econ. Entomol. 91, 572–578 (1998).Article 

    Google Scholar 
    30.Andow, D. A. & Alstad, D. N. Credibility interval for rare resistance allele frequencies. J. Econ. Entomol. 92, 755–758 (1999).Article 

    Google Scholar 
    31.R Core Team. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2020). https://www.R-project.org/.32.Herzog, D. C. Sampling soybean looper on soybean. In Sampling Methods in Soybean Entomology (eds. Koogan, M. & Herzog, D. C.) 141–168 (Springer, 1980).33.Navarrro, F. R., Saini, E. D. & Leiva, P. D. Clave pictórica de polillas de interés agrícola. INTA, EEA PERGAMINO, Pergamino, Buenos Aires, Argentina (2009).34.Gilligan, T. M. & Passoa, S. C. LepIntercept–An identification resource for intercepted Lepidoptera larvae. (Identification Technology Program (ITP), 2014). http://idtools.org/id/leps/lepintercept/key.html (2014).35.Chen, D. et al. Bacillus thuringiensis chimeric proteins Cry1A. 2 and Cry1B. 2 to control soybean lepidopteran pests: New domain combinations enhance insecticidal spectrum of activity and novel receptor contributions. PLoS ONE 16, e0249150 (2021).CAS 
    Article 

    Google Scholar 
    36.SAS Institute Inc., SAS/STAT: 9.1. Statistical Analysis System: getting started with the SAS learning. SAS Institute, Cary, NC (2002).37.Farias, J. R. B., Nepomuceno, A. L., Neumaier, N. Ecofisiologia da soja. Londrina: Embrapa-Centro Nacional de Pesquisa de Soja, Circular Técnica, 48 (2007).38.Sosa-Gómez, D. R. et al. Manual de Identificação de Insetos e Outros Invertebrados da Cultura da Soja. (Embrapa Soja-Documentos (INFOTECA-E), 2014).39.Kaster, M. & Farias, J. R. B. Regionalização dos Testes de Valor de Cultivo e Uso e da indicação de Cultivares de Soja-terceira Aproximação (Embrapa Soja-Documentos (INFOTECA-E), 2012).40.IBGE. Divisão Regional do Brasil em Mesorregiões e Microrregiões Geográficas. v. 1. (IBGE, 1990).41.GraphPad Prism. GraphPad Prism version 8.1.2 for Windows, GraphPad Software, San Diego. www.graphpad.com.42.Huang, F., Parker, R., Leonard, R., Yong, Y. & Liu, J. Frequency of resistance alleles to Bacillus thuringiensis-corn in Texas populations of the sugarcane borer, Diatraea saccharalis (F.) (Lepidoptera: Crambidae). Crop Prot. 28, 174–180 (2009).CAS 
    Article 

    Google Scholar 
    43.Rodrigues-Silva, N. et al. Negative cross-resistance between structurally different Bacillus thuringiensis toxins may favor resistance management of soybean looper in transgenic Bt cultivars. Sci. Rep. 9, 1–9 (2019).CAS 
    Article 

    Google Scholar 
    44.Yano, S. A. C. et al. Tolerância de Anticarsia gemmatalis Hübner, Pseudoplusia includens (Walker) e Rachiplusia nu (Guenée) à proteína Cry1Ac. In: Embrapa Soja-Artigo em anais de congresso (ALICE). (Embrapa, 2012).45.Perini, C. R. et al. Genetic structure of two Plusiinae species suggests recent expansion of Chrysodeixis includens in the American continent. Agric. For. Entomol. 23, 250–260 (2020).Article 

    Google Scholar 
    46.Bueno, A. F. et al. Challenges for adoption of integrated pest management (IPM): The soybean example. Neotrop. Entomol. 50, 5–20 (2021).CAS 
    Article 

    Google Scholar 
    47.Carrière, Y. et al. Governing evolution: A socio-ecological comparison of resistance management for insecticidal transgenic Bt crops among four countries. Ambio 49, 1–16 (2020).Article 

    Google Scholar 
    48.Heinrichs, E. A. & Silva, R. F. P. Estudo de níveis de população de Anticarsia gemmatalis Hubner, 1818 e Plusia sp. em soja no Rio Grande do Sul. Agron. Sulriograndense 11, 29–35 (1975).
    Google Scholar 
    49.Guedes, J. V. C. Lagartas da soja: das lições do passado ao manejo do futuro. Rev. Plantio Direto 144, 10–22 (2015).
    Google Scholar 
    50.Silva, C. S. et al. Population expansion and genomic adaptation to agricultural environments of the soybean looper, Chrysodeixis includens. Evol. Appl. 13, 2071–2085 (2020).Article 

    Google Scholar 
    51.Bacalhau, F. B. et al. Performance of genetically modified soybean expressing the Cry1A. 105, Cry2Ab2, and Cry1Ac proteins against key Lepidopteran pests in Brazil. J. Econ. Entomol. 113, 2883–2889 (2020).CAS 
    Article 

    Google Scholar  More

  • in

    Crop production and nitrogen use in European cropland and grassland 1961–2019

    Method and data overviewThe main objective of this study was to assemble a consistent and detailed dataset of annual cropland N budgets in Europe 1961–2019. To this end, we collated a range of different datasets. The FAOSTAT database is the primary data source, providing the longest and most complete data on crop and livestock production with global coverage. However, limiting our scope to the former EU28 allows us to leverage the Eurostat database which in certain respects has more comprehensive coverage. The FAOSTAT and Eurostat databases are both based on national statistics, but have several differences in scope and statistical nomenclature. Therefore, as this paper will demonstrate, the two databases in combination provide a richer picture than either one separately. In addition to these international datasets, we have compiled data from several national statistical databases and yearbooks as well as other literature sources. The main international databases used in this study are listed in Table 1. Remaining data sources are referred to in the text.Table 1 Main international databases used.Full size tableThe focus of this study is cropland, defined as land under arable crops (including temporary grassland) and permanent crops (e.g., fruit orchards). Cropland does not include permanent grasslands, defined by both FAOSTAT and Eurostat as grassland lasting for more than five years22,23. In some cases, however, the results obtained in this work have obvious relevance for permanent grassland as well. For example, our method to estimate synthetic N fertilizer input to cropland also produces an estimate of input to permanent grassland.As explained above, this study covers a territory which today comprises 26 countries. The territory is covered in the period 1961–2019 with the exceptions of Croatia, Estonia, Latvia, Lithuania, and Slovenia, which all gained independence in 1991/1992 and due to data limitations are covered in the period 1992–2019. We follow the nomenclature and regional divisions of the FAOSTAT database. For the period 1961–1999, the FAOSTAT database reports Belgium-Luxembourg as one unit, but since 2000 it reports Belgium and Luxembourg separately. Data are reported for Czechoslovakia 1961–1992, and separately for Czechia and Slovakia since 1993. For the whole period 1961–2019, data are reported for Germany as one unit, the territory known as Germany since the 1990 reunification. In some cases, the FAOSTAT regional nomenclature differs from other data sources—notably, the Eurostat database covers Belgium and Luxembourg separately since the start in 1955, and only West Germany until 1989—and to reconcile such differences we have collected additional data and calculated area-weighted or quantity-weighted averages as appropriate to the extent possible. To facilitate further analyses, we also provide merged 1961–2019 results for Belgium and Luxembourg and former Czechoslovakia. Throughout the paper, we loosely refer to these geographical areas as “countries”.The remainder of this Method section describes and motivates the data and methods used to estimate the following quantities:

    Harvested areas and harvested crop N in 17 categories of arable and permanent crops.

    Symbiotic N fixation in pulses and forage legumes on cropland.

    Areas of cropland, permanent grassland, and cropland in use.

    Synthetic N fertilizer use, partitioned between cropland and permanent grassland.

    Manure N flows: excreted N, partitioned into grazing and in-house excretion; losses of N (mainly ammonia) from manure management in animal houses and manure storage facilities; quantities of manure N that finally reach cropland and permanent grassland from grazing animals or through field application.

    Atmospheric N deposition to cropland and permanent grassland.

    Figure 5 gives a high-level overview of the main data sources, transformation steps, and results.Fig. 5Overview of the main data sources and results of this study. Major input datasets colored blue. Major derived results colored orange.Full size imageTable 2 gives a list of symbols and abbreviations.Table 2 This list covers all non-standard abbreviations and symbols used in more than one subsection of the paper or in the data record.Full size tableThe input data and results described in this paper as well as source code for all the calculations have been archived as a public data record12.Crop areas and harvests: overviewWe combined a range of data sources to estimate crop areas and N harvests in 17 crop categories.For most arable and permanent crops, we used crop harvests and areas from the FAOSTAT database (see Table 1).The only major crops missing from the FAOSTAT database are fodder crops such as temporary grassland, forage legumes, cereal crops harvested green, and fodder roots and cabbages24. For these fodder crops, we instead used data from Eurostat’s Annual Crop Statistics (ACS) (see Table 1) and a range of other sources discussed in detail below. The few data on green fodder crops reported in FAOSTAT’s database (only green maize) were excluded to avoid any double-counting.These data were processed in several steps to identify and address data quality issues. The main steps of this process are illustrated in Fig. 6, and a detailed description is given in the following sections. The final result of the process is a dataset of areas and N harvests in 17 crop categories (Table 3), covering the entire time period defined for each of country. The resulting dataset is found in the data record12.Fig. 6Illustration of main data sources and transformation steps used to estimate areas and N harvests in 17 crop categories. Major input datasets colored blue. Major derived results colored orange. Intermediate transformation steps in gray.Full size imageTable 3 The crop categories resulting from the crop data processing.Full size tableRationale for crop categorizationThe crop categories listed in Table 3 were chosen to produce a dataset that gives comprehensible and agronomically relevant information about major trends in crop mix and productivity. Specifically, the categorization was made based on the following considerations:

    The categories should contain crops with similar N yields and similar N yield changes over time. For example, wheat and grain maize have had a substantially steeper yield increase than other cereals over time and are therefore reported separately.

    The categories should be well-known categories of crops to simplify interpretation and comparison with other datasets. For example, although sugar beets and potatoes are comparable to cereals in terms of N yields and could possibly be grouped according to the previous criterion, we separated them since they are typically separated in agricultural statistics and models.

    The categories should make visible the characteristic differences in crop mix between countries. For example, we separated olives and grapes from other permanent crops since these two crops cover 20–30% of the cropland in four Mediterranean countries, compared to about 9% of the total European cropland.

    The number of categories should not be too large because with very small categories the signal-to-noise ratio declines, making any statistical analysis more difficult.

    Crop areas and harvests except fodder cropsFrom the FAOSTAT crop database, we extracted data on 121 arable and permanent crops. We converted the reported harvests to N quantities assuming crop N contents from Lassaletta et al.17. The same crop N contents are assumed for all countries in the whole period 1961–2019, even if it is likely that crop N contents have varied both geographically and over time, for several reasons. For example, N contents will tend to increase with dry weather and low yields; N contents will tend to increase with fertilizer rates; and N contents may be may be changed both up and down because of crop breeding. However, a more detailed analysis of these effects has not been possible in the scope of this study.As described above, we then aggregated these crops to 12 categories listed in Table 3.While the FAOSTAT database has a good level of coverage for most major crops (apart from fodder), there are sometimes data gaps where harvested quantities are reported without corresponding areas. This occurs mainly before 1985 in permanent crops (including olives) and to a lesser extent in the categories “Oilseeds” and “Vegetables and other”. For these crops, summing the individual crop areas in a country-year would sometimes result in considerable underestimation of the total harvested area. Similarly, summing areas and harvests separately for available crops would result in incompatible estimates of category areas and harvests.In order to ensure consistent estimates of category areas and harvests, we therefore took the following approach.For each crop category, we calculated the sums of available crop areas (Asum) and crop N harvests (Hsum). In addition, based on the crops where both areas and harvests were available, we calculated category-level weighted average N yields (Yest). Since some area data were missing for individual crops, we additionally calculated an estimate of total category area Aest = Hsum/Yest. When crop harvests are available but some crop areas are missing, the estimate Aest is equivalent to assuming that the weighted average N yield of crops with missing areas is equal to the weighted average yield of the other crops in the same category.For each country and category, we then generated and inspected figures showing the time series of the category-level variables Asum, Aest, Hsum, and Yest along with available crop-level data on area, harvests, and yields. If all crop categories had been present in all countries, there would have been 28 × 12 = 336 figures. However, some country × category combinations are absent (e.g., olive trees in Sweden) and in total there were 304 figures to inspect (available in the data record12).Based on visual inspection of these figures, and in some cases cross-checking against other sources, we then chose on a case-by-case basis how to estimate the category area and harvest from available data. Note that we used available crop-level data to fill category-level data gaps. Filling all the data gaps on crop level have been very laborious and also unnecessary, since the aim was merely to make complete and consistent estimates on crop category level.By default we used Hsum as estimate of the category harvest and Aest as estimate of the category area. We used these default estimates when they looked fairly smooth over time and no obvious and serious data gaps were present in the crop-level data.However, sometimes there were data gaps in the crop-level data that motivated other adjustments on the category level. For example, one type of problem is when a category-level variable is completely missing. The most important example is that for olives in Greece, Portugal, and Spain, harvest data are available since 1961 but area data only since the 1980s. For olives in these three countries, we collected data based on national statistics12,25,26,27,28. In addition, we filled a small number of minor data gaps (e.g., harvests or areas completely missing during one or a few years) using constant extrapolation. A similar type of problem is when some crops within a category lack area and/or harvest data during a part of the period. This can create a false impression that the category’s N yield has suddenly changed abruptly. In a few such cases it appeared more plausible to extrapolate (or sometimes interpolate) areas or yields to obtain a dataset covering the whole period. Generally, these various adjustments made to the default estimates Hsum and Aest were rather small. The missing olive areas in Greece, Portugal, and Spain were by far the most important adjustments made in this process. In each country, there was less than one percent difference between the default estimate Hsum and the final adjusted estimate of total N harvests. Prior to 1985, these adjustments increased the total crop area by about 1.5% on average across all countries.Finally, we also explicitly assigned zero harvests and areas where data were completely missing, to obtain a full dataset of the 12 FAOSTAT crop categories in each country.Fodder crop areas and harvests: overviewAs mentioned above, the FAOSTAT crop production database excludes most fodder crops. Here, we instead assembled a fodder crop dataset using Eurostat crop statistics (see Table 1) and other sources.Eurostat reports areas and harvests of arable and permanent crops in a hierarchy of crop codes23,29. The most important category of fodder crops in this hierarchy is “Plants harvested green from arable land” (crop code G0000), which is further subdivided in a number of subcrops as shown in Table 4. In addition, we included the Eurostat category “Other root crops n.e.c.” (R9000) which mainly accounts for fodder roots, for example Beta vulgaris (known by many names, including fodder or forage beet, or mangold, mangelwurzel, etc.) and several Brassica species (rutabaga/swede, turnip, etc.)23. Crop code R9000 does not include roots for seed or human consumption. Other crops used completely or partly for animal feed, including grain legumes, cereals harvested for grain, sugar beets, potatoes, etc., are accounted for in the FAOSTAT database30. To our knowledge, the crop codes G0000 and R9000 together account for all the major European fodder crops not included in the FAOSTAT database.Table 4 Excerpt from Eurostat’s hierarchy of crop codes23.Full size tableIn the following sections we describe in detail how we combined data from Eurostat with other sources corresponding to the Eurostat crop codes G1000, G2100, G2900, G3000, G9000, and R9000. The results are available in the data record12.Fodder crop areasData extraction from the Eurostat ACS databaseHarvested areas are reported in Eurostat’s ACS database23. Data coverage in the Eurostat ACS data varies widely. For some countries, especially the early members of EEC and EU, the area data are complete back to the 1950s. For the more recent EU member states, the data coverage typically starts around the time of their accession application to the EU. For Croatia, Estonia, Latvia, Lithuania, and Slovenia, which in this study are covered starting in 1992, the data coverage is fairly complete. However, for the former communist states of Bulgaria, Czechoslovakia, East Germany, Hungary, Poland, and Romania, which in this study are covered starting in 1961, there are no data prior to 1987 in the Eurostat database. Eurostat also lacks data during some periods for several countries in western Europe. In a few cases we cross-checked suspected errors and gap-filled area data from the Eurostat Farm Structure Survey (FSS)29 (see Table 1). However, the FSS data generally have smaller coverage than the ACS and are not entirely comparable in scope and methods, so we used it very sparingly.Some special treatment was needed for crops G9000 and R9000. Crop code G9000 is not reported in the Eurostat ACS, but we summed it from available data of G9100 and G9900. The reason to merge these two crop codes is that their reported areas often fluctuate in such a way to suggest that the same crops have been reported variably as G9100 or G9900; thus, data gaps for the combined G9000 area are fewer and easier to fill than for the individual G9100 and G9900 areas. For R9000 (fodder roots), data were almost never explicitly stated, but could in many cases be calculated as R0000–R1000–R200023.Gap-filling and other adjustments to the Eurostat fodder crop areasIn most countries, the Eurostat data on fodder areas are fairly smooth, complete, and internally consistent since around year 2000. Before this period, several countries have data gaps and/or report large, sudden changes which we intepreted as potential errors. A reason to expect some errors and inconsistencies is that the nomenclature used in older national statistics likely is incompatible with the current Eurostat crop nomenclature, which may cause problems in the translation of old data to the Eurostat database. However, since considerable shifts in fodder crop areas actually have occurred since 1961 in most European countries, it is not always straightforward to determine whether abrupt changes in reported areas are reporting errors or accurate representations of historical developments. We therefore scrutinized and cross-checked the Eurostat data against other sources, filling data gaps and making other adjustments to reconcile major discrepancies. The collected dataset on fodder crop areas, as well as figures showing the stepwise gap-filling of fodder crop areas, are available in the data record12.Fodder roots were important crops in the first half of the 20th century in several European countries, but areas then declined as they were replaced by other, less labor-intensive fodder crops31. Therefore we paid special attention to filling data gaps in R9000 areas during the 1960s–1970s.Before listing the data sources and adjustments country by country, we specifically mention the common approach used to fill the long 1961–1986 data gaps in Bulgaria, Czechoslovakia, East Germany, Hungary, Poland, and Romania. We mainly used data from reports of the Economic Research Service of the US Department of Agriculture (USDA ERS), which during the 1960s–1980s collated information from the statistical yearbooks of the socialist states in a series of reports32,33,34,35,36. These reports cover the years 1960 and 1965–1987, and give areas for three categories of fodder crops: “feed roots”, “corn silage”, and “hay”. We assigned the former two crop codes R9000 and G3000, which in the overlapping year 1987 agreed perfectly with the Eurostat ACS. The last category, “hay”, is more complicated: it may refer to a mix of annual and perennial crops harvested green, predominantly forage legumes in pure stands or mixed with grass, cereals and cereal/legume mixtures, and possibly pure grass cultivation on arable land. The USDA ERS “hay” category clearly excludes harvests from permanent grassland. Since there is usually one year of data overlap between the USDA ERS statistics and the Eurostat ACS in 1987 for these countries, we could usually conclude that the “hay” area then corresponded to combination of crop codes G2100, G2900, and sometimes G9000. Country by country, we decided on a combination of Eurostat crop codes to match against the “hay” area, and then divided the 1960–1986 “hay” area between them in proportion to their their 1987 areas. Temporary grasslands (G1000) account only for a few percent of the fodder crops in most of Eastern Europe, and we therefore mostly extrapolated the earliest available G1000 areas back to 1961. Country-specific details are elaborated below.The remainder of this section lists the data sources and adjustments country by country. For brevity, we omit some descriptions of the following minor adjustments: interpolation of minor data gaps, sometimes using data from 1960 or 2020; extrapolation of minor fodder crops accounting for a small share of the total fodder area; removal of obvious outliers.Austria. The Eurostat ACS data are complete and consistent since the start in 1980. Data for 1960 and 1970–79 were filled using national statistics and data from the FAO 1960 World Census of Agriculture37,38,39. Remaining data gaps interpolated and extrapolated.Belgium and Luxembourg. The Eurostat ACS data are almost perfectly complete and consistent since the start in 1955. Minor data gaps in G2100 and G2900 areas interpolated.Bulgaria. The Eurostat ACS data are almost complete and consistent since the start in 1987. In 1960–1986, we used G3000 and R9000 areas from USDA ERS publications34,35,36 as explained above. In 1987, Eurostat’s combined area of G2100 and G2900 matches the USDA ERS “hay” area, so we divided the 1960–1986 hay area between these crop codes in proportion to their 1987 shares, and extrapolated G9000 and G1000 values constant to 1961.Croatia. The Eurostat ACS data are complete from the start in 2000. We extrapolated the areas back to 1992.Czechia. The Eurostat ACS data are complete from the start in 1987 apart from G1000 which is reported at around 10% of the fodder area since 2011. A lone G1000 value in 1999 is conspicously close to the temporary 1995–1999 decrease in G2900, suggesting a temporary classification change in temporary grassland and legume-dominated crops. Since G1000 data are largely missing, we chose to discard the 1995–1999 decrease in G2900 area and interpolate surrounding values while extrapolating the 2011 G1000 area constant back in time to 1987.Slovakia. The Eurostat ACS data are complete and consistent from the start in 1987 except a minor gap in G2100 and G2900 areas which we interpolated.Czechoslovakia. Areas for 1987–1992 were taken as the sum of the adjusted Eurostat ACS data for Czechia and Slovakia. In 1987, Eurostat’s combined area of G2100, G2900, and G9000 matches the USDA ERS “hay” area, so we divided the 1960–1986 hay area between these crop codes in proportion to their 1987 shares. G1000 was extrapolated constant to 1961.Denmark. The Eurostat ACS data are fairly complete and consistent back to 1955, except for the G9000 area which fluctuates considerably. A closer inspection shows that 1973–2009, G9900 has a large area share while G9100 is not reported; and from 2010 the G9900 area is zero while G9100 has a smaller share. Before 1973, the Eurostat database does not report G9100 or G9900 areas. National statistics show that this apparent discontinuity arises because the reported G9900 area for some years, in addition to cereals harvested green, also includes aftermath, i.e. late season harvests or grazing after other crops. The aftermath at its peak in year 2000 accounted for about half the reported fodder area on arable land but less than 10% of the harvested feed value40,41. Considering the incomplete data coverage and the minor importance of the aftermath in terms of harvested quantities, we replaced the Eurostat G9000 area by a complete record of cereals harvested green (i.e., corresponding to G9100) based on national statistics42,43. The national statistics prior to 1982 report the combined area corresponding to G3000 + G9100, so to avoid double counting we subtracted the G3000 area as reported by Eurostat.Estonia. The Eurostat ACS data are mostly complete from 1991. A data gap in the G1000 area was filled by difference since the G2000 data clearly includes the later G1000 area until 2003. We also filled minor data gaps in G2100 and G2900 using national statistics44.Finland. The Eurostat ACS data are complete since 1998. We used national statistics45 to fill the G1000 area, which in 1998 covered more than 95% of the fodder area. Fodder roots made a minor contribution in Finland even in the 1950s and 1960s46 when they were much more common in other countries. The main feed root seemingly was potatoes46,47, which is already accounted for in the FAOSTAT crop database30. Considering the lack of further data and the dominance of G1000 in the fodder production we extrapolated the 1998 area of other fodder crops back to 1961.France. The Eurostat ACS data are complete since 1961 apart from a few minor data gaps which we interpolated.Germany. The Eurostat ACS data cover all the fodder crops starting in 1955, but geographically covers only West Germany until 1989. To complete the period 1961–1989, we estimated East Germany’s fodder crop areas in 1989 data from Eurostat ACS and USDA ERS34,35,36 data as follows. We estimated East Germany’s fodder crop areas in 1989 as the Eurostat increment in fodder crop areas 1989–1990, an estimate which builds on the assumption that both West and East German fodder crop areas were approximately constant 1989–1990. East Germany’s 1960–1987 areas of G3000 and R9000 were then taken from the USDA ERS data. The estimated 1989 area of G1000, G2100, G2900, and G9000 matched the USDA ERS 1987 “hay” area, so we divided the 1960–1987 hay area between these crop codes in proportion to their 1989 shares. The remaining minor data gaps were interpolated.Greece. While the Eurostat ACS data appear complete and internally consistent since year 2000, they are difficult to reconcile with older Eurostat data and data from other sources. For the period 1969–1986, the Eurostat ACS data suggest that alfalfa (G2100) is the dominating arable fodder crop varying around 150–200 thousand hectares (kha). Similar alfalfa areas are reported from the 1950s until 2017 in multiple overviews of Greek fodder production as well as recent national statistical publications48,49,50,51,52,53. In the Eurostat database, however, there are almost no fodder crop data reported for 1988–1999, and from year 2000 the alfalfa area is reported around 10–15 kha. For other forage legumes, the Eurostat data consistently report zero area, while several sources report fairly constant areas in the range 35–70 kha48,50,51,52,53. The area data for fodder crops in the Eurostat FSS 1990–201629 are incomplete and offer few additional insights. Based on these considerations, we made the following adjustments to the Eurostat data: discarded the Eurostat area data for G2100 starting in year 2000 and extrapolated the fairly constant areas 1969–1986 to the whole period 1961–2019; assumed a constant G2900 area of 50 kha; interpolated the G9000 area data which appear broadly consistent with various data sources; and extrapolated the average G1000 area from year 2000 back to 1961.Hungary. The Eurostat ACS data on Hungary’s main fodder crops, forage legumes and green maize, are mostly complete from 1987. The USDA ERS data34,35,36 provide areas of G3000 and R9000 back to 1960. The G9000 area is missing 1987–1997. Since the USDA ERS “hay” area exceeds the combined G2100 and G2900 area in 1987, we estimated the G9000 area as the remainder of the hay area in 1987, and then divided the 1960–1986 hay area between G2100, G2900, and G9000 in proportion to their 1987 shares. The resulting G2100 area estimate in 1961–1987 agrees remarkably well with data from the Hungarian Central Statistical Office54. For G1000, the few available data indicate that the area is very small compared to the total fodder area, and we extrapolated the 2003 area constant back to 1961.Ireland. Temporary grassland (G1000) clearly dominates fodder production on cropland. However, a considerable area of temporary grassland has been reclassified to permanent grassland around 1996–2016 although sources disagree on the exact size and timing of this shift. The Eurostat ACS data reports a drop from about 0.7 Mha to 0.1 Mha between 1996 and 2000, and a corresponding increase in permanent grassland (crop code J0000) between 2001 and 2002. The Eurostat FSS reports a similar drop in G1000 area between 2013 and 2016. The FAOSTAT database reports a simultaneous change in both temporary and permanent grassland between 2006 and 2008. For consistency with the FAOSTAT permanent grassland area (details below) we used FAOSTAT’s temporary grassland area data during 2000–2006, the only period where it disagrees with the Eurostat ACS. For the remaining fodder crops, area data are almost complete from the start in 1955.Italy. The Eurostat ACS data are mostly complete for G2100, G2900, G3000, and R9000 since 1970. We filled a 1989–2013 data gap in G2900 by interpolation. For G1000 and G9000, the Eurostat annual crop statistics are incomplete and somewhat erratic. In the years where data for both G1000 and G9000 are available, their sum shows a smooth decline from around 1.2 Mha in the 1980s to about 0.8 Mha in 2014, suggesting that crop classifications may have changed more than actual areas. National statistics from 2006–201855 support this interpretation, although the national statistics cannot easily be matched to the Eurostat crop codes. A detailed expert overview from 197756 agrees well with the G1000 and G9000 data for the 1970s. Considering all this, we filled the G1000 and G9000 data gaps by interpolating their area sum between 1986 and 2014 and dividing the resulting total in proportion to their 1986 areas. Remaining data gaps were filled by extrapolation from the 1970s back to 1961.Latvia. The Eurostat ACS data are mostly complete for G3000, G9000, and R9000 since 1987. Areas of G1000 are reported since year 2000, completely dominating the fodder area. In 1987–1999, reported areas for G2000 appear to include what is later reported as G1000, suggesting a change in crop classifications. We extrapolated the small areas of G2100 and G2900 reported since 2015 constant back to 1992 and assigned the remainder of the 1992–1999 G2000 area to G1000.Lithuania. Like in Latvia, the Eurostat ACS data reports no G1000 area prior to 2000, but a G2000 area which likely includes the later G1000 area. We extrapolated the 2001 areas of G2100 and G2900 back to 1992 and assigned the remainder of the 1992–1999 G2000 area to G1000. This leads to an estimated 95% decrease in the G1000 area between 1999 and 2001, which agrees with an increase in Lithuania’s permanent grassland area registered in the Eurostat and FAOSTAT databases since 2001. Most likely this accurately reflects a heavy decrease in the resowing of grasslands in Lithuania which occured after the collapse of the Soviet Union57. We extrapolated the G1000 area from 1999 to 2000 to agree with the timing of this change in the reported permanent grassland areas.Netherlands. The Eurostat ACS data are complete since 1955 except for minor gaps which we interpolated. The data show an abrupt increase in temporary grassland (G1000) area from less than 40 kha in the mid-1990s to about 200 kha in the mid-2000s, a change which is also reflected in decreased permanent grassland areas in the Eurostat and FAOSTAT databases.Poland. The Eurostat ACS data are mostly complete since the start in 1987, but only partly consistent with other data sources. The USDA ERS data for 1960–198734,35,36 precisely match the G3000 and R9000 areas in 1987 and we therefore used these without modification. However, the USDA ERS “hay” area in 1987 is considerably smaller than the combined Eurostat G2100 and G2900 areas, which suggests a reporting error in at least one of the datasets. Several factors strongly suggest that the Eurostat G2100 and G2900 areas are both incorrect in 1987–2001. First, Eurostat’s reported G2900 area falls by more than 1.5 Mha in 1987–1998, a change rate which appears unlikely even given the rapid changes taking place in Poland starting in the late 1980s. Second, national statistics from Polish statistical yearbooks58,59,60,61 report that the area of fodder legumes has never been as high as 1.5 Mha in Poland, and specifically suggest a mistake in crop code assignments since Eurostat’s G2100 areas from 1987 to 2001 exactly equal the total areas of perennial legumes according to national statistics. Third, area data from the Eurostat FSS in 2003 also suggest that Eurostat ACS data for G2100 and G2900 are incorrect until year 2001. Based on this, we discarded the Eurostat ACS areas of G2100, G2900, and G9000 in 1987–2001.The USDA ERS “hay” area varies between 1.4 and 1.8 Mha in 1960–1987. The Polish statistical yearbooks do not give sufficient information to divide this between crop codes G1000, G2100, G2900, and G9000. However, data from the FAO 1960 World Census of Agriculture39 shows that clover was the main forage legume (around 0.6 Mha). Alfalfa covered about 0.13 Mha, about 13% of the pure forage legume area. Throughout the 1960s–1980s, the statistical yearbooks show that perennial legumes covered about half the hay area58,59. Similarly, an expert summary of national statistics in 1965 loosely described the fodder area except maize and fodder roots as consisting of 48% pure forage legumes and 52% of clover/grass and pure grass62. In line with this, Eurostat’s G1000 area for 1987 corresponds to 29% of the 1987 hay area. Based on these data, we divided the USDA ERS hay area in 1960–1987 using the following fixed proportions: 29% G1000, 6% G2100, 42% G2900, and the remaining 23% G9000. We interpolated the remaining data gaps.Portugal. The Eurostat ACS data are incomplete. The total G0000 area is reported roughly constant since 1978. An almost constant share of around 9% G1000 and 18% G3000 is reported since 1991. The major area appears to be G9000, but its area is only reported since 2011. Quantitative data from other sources appear to be scarce63. A paper64 from 1990 describes cereal/legume mixtures (i.e., G9000) as the main arable fodder, which agrees with recent Eurostat data. Considering the near-constant G0000 area since 1978 and the near-constant shares of different fodder crops, and that few additional data could be found, we extrapolated available data constant back to 1961.Romania. The Eurostat ACS data are complete and consistent since the start in 1987. Areas of G3000 and R9000 1960–1987 were filled from USDA ERS publications34,35,36 as explained above. In 1987, Eurostat’s combined area of G2100, G2900, and G9000 precisely matches the USDA ERS “hay” area, so we divided the 1960–1986 hay area between these crop codes in proportion to their 1987 shares, and extrapolated the very small G1000 area constant from 2005 back to 1961.Slovenia. The Eurostat ACS data are almost complete from the start in 1991. Some minor data gaps were filled using national statistics65.Spain. The Eurostat ACS data are almost complete since 1965. We extrapolated the 1965–66 areas back to 1961.Sweden. The Eurostat ACS data are fairly complete from the start in 1992. We filled the data gaps back to 1961 using national statistics66.United Kingdom. The Eurostat ACS data are almost complete since the start in 1955. The main exception is fodder roots (R9000), for which Eurostat data are available since 2000. We completed the record using data from national surveys67,68,69 covering all major fodder roots: turnips, swedes, and fodder beets (including mangolds, which are distinct from fodder beets in British terminology31). For G9000, a minor issue is that the 2010–2011 areas are reported identical to G3000 areas, likely by mistake. We discarded these data and filled by interpolation. We also extrapolated G9000 constant from 1970 back to 1961.Fodder crop harvestsWe used Eurostat ACS data to estimate fodder crop N harvests. The ACS production data have two important limitations: (1) they are incomplete, even more so than the area data, and (2) they have several inconsistencies related to the water content of the harvest (details below). For these reasons, it was not possible to establish time series of fodder crop yields for more than a few crop/country combinations. In the few cases where long-term time series are available, they however show that fodder crop yields on average have increased relatively slowly.Considering the lack of data, and that the period 2000–2019 has the best data coverage, we decided to estimate country/crop specific yields in 2010 from available data. Between 1961 and 2019, we assumed based on long-term statistics from Austria, France, Hungary, Italy, Poland, and Sweden28,37,56,58,59,70,71,72,73,74 a linear increase such that the 1961 yields are 75% of the 2010 yields. The following subsections describe the estimation of 2010 yields from available data.Estimating the dry matter yields of fodder cropsThe nominal water content of the harvest data is a central concern since harvests may be reported with different nominal water content. Some countries report in dry matter (0% water), while others use crop-specific nominal water content typically between 12% (hay) and 65–80% (e.g., green maize and other silage crops).The Eurostat ACS data since year 2000 accounts for these differences by reporting the water content (“humidity” in Eurostat’s nomenclature) along with the harvested quantities28. There are two different datasets, one in national humidity (0–88% water content) and one in EU standard humidity (always 65% water content for plants harvested green)28. The coverage of humidity values (since year 2000) is not complete, which means that sometimes the harvest is only given in national humidity basis. Prior to year 2000, all the data are given in national humidity basis, without corresponding humidity values.To construct a harmonized dataset of dry matter yields for each country/crop combination, we studied and compared four time series of available yield values: (1) in national humidity, (2) in EU standard humidity, (3) in dry matter based on the national humidity dataset, and (4) in dry matter based on standard EU humidity dataset. For crop code G9000, we calculated production-weighted average yields from data on crop codes G9100 and G9900. Based on a close inspection of these data and sometimes cross-checking against national data sources, we identified the following types of possible reporting errors:

    Some crop yields are reported equal in national and EU humidity, although the corresponding humidity values differ. This creates two different dry matter yields, at most one of which could be accurate. This possibly reflects a mistake in the conversion between national and EU humidity. In these cases we typically used the national data.

    Sometimes, the reported yield of a crop changes drastically from one year to the next such that the older yield in national or EU humidity roughly equals the new yield in dry matter, or vice versa. In these cases, one of the two yield levels could typically be ruled out as implausible.

    Some calculated dry matter yields seem implausibly low. In some of these cases it could be deduced that a yield reported with a nominal water content of 65% was actually in dry matter or hay basis (i.e., about 85% dry matter).

    Some yield values seem implausibly high or low compared to neighboring years or similar countries without any apparent reason.

    Based on these considerations, we selected a subset of yield values for each country/crop, which we then averaged to an estimate of the 2010 (reference year) yield. We aimed to use 2000–2019 data if possible, not only because they best represent the 2010 yields but also because the more recent data are more complete and consistent. We used data from the 1990s in a few cases where no 2000–2019 were available. We used dry matter yields with only a few exceptions: for the crops G1000, G2100, and G2900, we sometimes used yields in national humidity basis if these appeared to be reported as hay or dry matter, assuming a dry matter content of 85%; and for fodder roots (R9000), humidity values are not reported and we uniformly assumed a dry matter content of 16%75. Figures illustrating data selection are available in the data record12.Accounting for grazing on croplandWe considered the complication that cropland, especially temporary grasslands (G1000), to some extent is grazed in addition to the mechanical harvest. The harvest statistics for temporary grassland appear to account only for mowing which means that they underestimate the total crop production. Mixed mowing and grazing appears to occur in temporary grasslands throughout Europe76,77, but quantitatively it is probably most important in the Nordic countries where temporary grassland occupies a considerable share of the cropland and is grazed fairly commonly.In fact, grazing may also occur in several of the fodder crops as well as in other crops, between or after harvests. However, we consider temporary grasslands as the probable main source of grazing intake on cropland, and considering the lack of data on this topic we make an estimate of grazing intake on cropland accounting only for temporary grassland.Relevant data to accurately estimate the grazing component of temporary grassland production are very scarce, but a recent investigation of Swedish data shows that grazing contributes about 20% in addition to the mechanical harvest of temporary grassland78. At least in Finland and Sweden, similar proportions of temporary grassland are used exclusively for grazing45,79. Considering that no further information could be found, we inflated the G1000 yield estimates by 20% in all the countries.Filling of remaining data gaps in fodder crop yieldsA few remaining data gaps were filled by extrapolating crop yields from neighboring countries with similar climate and agricultural productivity. In some cases we averaged yield values from multiple neighboring countries. The following data were extrapolated from other countries:

    Belgium (G9000): average of Luxembourg and the Netherlands

    Czechoslovakia (all fodder crops): average of Czechia and Slovakia

    Germany (G9000): from France

    Greece (all fodder crops): from Bulgaria

    Ireland (G1000, G2900, G9000): averages of Belgium, France, Germany, and the Netherlands

    Italy (G2100, G9000): average of France and Austria

    Latvia: (G2100, G2900): from Lithuania

    Luxembourg (R9000): average of France and the Netherlands

    Portugal (G1000, G9000, R9000): from Spain

    Sweden (R9000): average of Latvia and Lithuania

    United Kingdom (G1000, G2100, G9000, R9000): averages of Belgium, France, Germany, and the Netherlands

    Estimation of fodder crop harvestsWe finally averaged the selected yield values for each country/crop, thus producing the estimates of 2010 dry matter yields. From these, we estimated N yields assuming N contents listed in Table 4.Finally, we calculated the N harvest for each of the fodder crops in each country by multiplying the estimated yield time series by the gap-filled area time series.Aggregation of fodder areas and harvests to categoriesThe six Eurostat fodder crops (Table 4) were ultimately aggregated to five categories following the considerations described above. Specifically, we aggregated alfalfa (G2100) and other forage legumes (G2900) into one category.Cropland and grassland areasSeveral of the results calculated in later sections depend on the total cropland area as well as on the areas of permanent and temporary grassland. This section explains how we estimated these areas from available data sources.As an estimate of cropland area in use, Lassaletta et al.17 used the sum of harvested crop areas, except when this sum exceeded the total cropland area reported in the FAOSTAT database, in which case the FAOSTAT cropland area was used instead. The sum of harvested crop areas can exceed the FAOSTAT cropland area due to multicropping or intercropping. We followed the same approach as Lassaletta et al., with the only difference that we used the sum of adjusted crop categories (henceforth called the crop area sum) rather than the sum of individual FAOSTAT crops. Some examples of the crop area sum compared to the FAOSTAT cropland are shown in Fig. 2.Permanent grassland areas, following Lassaletta et al.17, were taken from the FAOSTAT database (“Land under perm. meadows and pastures” in the land use dataset). These areas were used in estimates of the allocation of synthetic and manure N inputs to cropland (details below).The FAOSTAT land use data at the time of our study was only available until 2018. We considered two options for filling the data gap in 2019. Eurostat data on permanent grassland and cropland could be used, or the FAOSTAT data for 2018 could be extrapolated constant to 2019. In general, the two methods would produce very similar results. However, since the FAOSTAT and Eurostat data series on permanent grasslands in a few countries differ substantially, we chose to extrapolate the FAOSTAT land use data constant from 2018 to 2019 as it probably gives the most internally consistent results.Temporary grassland areas were taken equal to the gap-filled crop category Temporary grassland (see Table 3).Assessing the accuracy of the estimates of cropland in useOur approach may lead to an overestimate or an underestimate of actual cropland area in use. The crop area sum may overestimate the cropland in use if some areas are counted more than once due to multicropping or intercropping. This is the reason to set the FAOSTAT cropland area as an upper limit to cropland in use. The reported cropland area may also overestimate the cropland in use since it may include considerable areas of fallow land.In principle, a more direct estimate of cropland in use would be the reported cropland area minus fallow land, but this is not an option since available data on fallow land areas in the FAOSTAT and Eurostat databases are much too incomplete to cover the whole 1961–2019 period.Instead, to test the accuracy of our estimates, we compared the crop area sum to cropland minus fallow land using data from the Eurostat database which has the most complete records of fallow land. The most common result of this test is that the areas match within a few percent (less than ± 5% difference in 75% of 971 country-years with data on fallow area). See figures in the data record12 for details.We also compared the crop area sum to the FAOSTAT cropland area. The crop area sum exceeds the FAOSTAT cropland by at least 1% in 163 country-years, about 12% of all 1308 country-years. The only major exceedances, e.g., more than 10% exceedance for at least 3 years, are in the 1960s and 1970s in Italy, Portugal, and Romania. See figures in the data record12 for details. At least in Romania most of the difference is explained by inter- and multicropping of cereals with beans and squash or pumpkins which was a common but gradually decreasing practice during the 1960s and 1970s80,81. In general, however, we have not been able to systematically determine the extent of inter- and multicropping.In summary, these results show that the crop area sum is usually a good approximation of cropland in use in Europe, especially after the 1970s.Symbiotic N fixationWe estimated symbiotic N fixation in pulses and forage legumes using the same method as Lassaletta et al.17, i.e., assuming a linear relationship between the fixed N and the harvested N yield Y,$${rm{BNF}}=Ycdot {rm{Ndfa}}cdot frac{{rm{BGN}}}{{rm{NHI}}},$$
    (1)
    where Ndfa is the share of plant N derived from the atmosphere, BGN is the ratio between total and above-ground plant N, and NHI is the N harvest index, i.e., the ratio between harvested and total above-ground plant N. We used the same crop-specific parameter values as Lassaletta et al.17,82.For crop code G1000 (temporary grassland, including grass/clover mixtures) we assumed that 25% of the dry matter harvest was forage legumes;5,17 for G2900 (pure forage legumes, sometimes mixed with grass) 90% forage legumes; and for G9000 (a variety of crops, including cereal/legume mixtures) 25% forage legumes (see also Table 5, and Eurostat’s ACS handbook23).Table 5 N contents assumed for the fodder crops.Full size tableThe resulting parameter values are available in the data record12.Atmospheric N depositionLassaletta et al.17 developed country-specific time series of N deposition rates in 1961–2013. We used these, extrapolated constant in the period 2014–2019, to estimate atmospheric N deposition input to cropland by multiplying by the estimated cropland area in use.Synthetic N fertilizer consumptionTime series of agricultural synthetic N fertilizer consumption in European countries are available in several international databases:

    FAOSTAT has a dataset representing agricultural use of fertilizers, based primarily on data reported by countries. FAOSTAT fills data gaps using various imputation methods, e.g., based on trade and production data83.

    Eurostat has a dataset representing agricultural use of fertilizers, based on data reported by countries84.

    Eurostat additionally publishes a dataset representing total national fertilizer sales, using data from Fertilizers Europe85.

    The International Fertilizer Association (IFA) publishes a dataset representing national consumption. It is based on sales data from the fertilizer industry, and gap-filled using production and trade data and N budgeting86.

    In principle, national sales of fertilizers is not equal to agricultural use because (1) stock changes between years can shift agricultural use compared to sales, and (2) sales also cover non-agricultural uses such as parks and lawns. However, due to lack of data, the international datasets of agricultural use do not systematically account for stock changes and non-agricultural use83,84.There is no immediately apparent way to judge which of these four datasets best represents the actual history of agricultural use. We therefore inspected and compared the four datasets to assess (1) to what extent they cover the studied time period, (2) whether the data seem consistent and plausible, and (3) how well the datasets agree with each other. The four datasets sometimes disagree considerably. By far, the longest and most complete datasets are provided by FAOSTAT and IFA. While these two datasets are mostly smooth and consistent, in a few exceptional cases they exhibit implausible jumps in consumption (e.g., by 50% or more) from one year to another. The FAOSTAT database has somewhat more such episodes, and moreover disagrees with the other three datasets more often than the others. While agreement between several datasets is no guarantee for their correctness, it seems likely that the majority vote of these partly independent estimates is the best estimate. The IFA dataset generally agrees well with the Fertilizers Europe dataset. The Fertilizers Europe dataset arguably exhibits the fewest implausible jumps but covers a shorter time period.Based on these results, we chose to use the Fertilizers Europe dataset85 where possible, and in second place the IFA dataset86. These datasets together (1) have almost complete coverage of the country-years in this study, (2) mostly appear consistent and plausible, and (3) usually agree well with a majority vote of the four datasets. We used FAOSTAT data for Finland 1961–1984 and Slovenia 1992–2006 because the IFA data appeared implausible. In addition, a small number of data gaps were filled using FAOSTAT data (Croatia 1992–1993, Czechia and Slovakia 1993) and Eurostat country consumption data (Belgium and Luxembourg 2000–2019).Figures showing the comparison of the four datasets and the selected data are available in the data record12.Share of synthetic N fertilizer to croplandFollowing Lassaletta et al.17 we calculated the synthetic N fertilizer input to cropland by multiplying the total consumption by country-specific time series of the shares applied to cropland. Although data on these shares are scarce, it is well known that several European countries have long histories of large synthetic N inputs to permanent grassland. In this paper, we present a revised estimate of these shares, following same approach as Lassaletta et al., but adding a significant amount of additional data for several countries. To make sense of the incomplete and sometimes contradictory data we used several techniques to interpret and gap-fill the data on a country-by-country basis. Figure 7 gives an overview of the steps taken. In the following subsections we first describe the overall process in more detail and then provide country-specific details.Fig. 7Illustration of main data sources and transformation steps used to estimate the application of synthetic fertilizer on cropland and permanent grassland. Major input datasets colored blue. Major derived results colored orange. Intermediate transformation steps in gray.Full size imageIdeally, the inputs of synthetic N fertilizer would be further disaggregated, e.g., to the 17 crop categories for which areas and harvests are reported in this study (Table 3). However, there are no datasets that enable comprehensive and reliable estimates of synthetic N fertilizer inputs on the level of individual crops or crop categories for the time period and countries covered in this study. As will be detailed below, pan-European crop-and-country-specific data on fertilizer inputs are available only for a few years between the 1990s and today. Further research on this topic would be valuable, and could use, for example, a combination of mass balances, statistical surveys, expert estimates, and crop-specific fertilizer recommendations from throughout the decades in different countries to establish plausible estimates of how synthetic N fertilizers have been allocated between crops. Such an exercise, however, could prove rather labor-intensive and moreover would necessarily imply a level of uncertainty which we have decided is not acceptable in the scope of this study.Figures illustrating the final results are available in the data record12.Data collection and interpretationThree main categories of data were used. First, a few countries (France, Ireland, the Netherlands, and the United Kingdom) have carried out repeated national statistical surveys on grassland and cropland fertilization. Although each of these datasets have their idiosyncracies and required additional data processing (see below), we used them where possible since they are the most consistent and complete datasets available. Second, crop-specific fertilizer datasets from the fertilizer industry have been published in collaboration with FAO during the 1990s and early 2000s87,88,89. In addition, a similar unpublished dataset from EFMA90 (now Fertilizers Europe) gives crop-specific fertilizer rates for the crop year 2005/2006. These datasets are extremely valuable since they cover all the countries in this study during at least one year. The crop-specific datasets were later aggregated and further processed as explained in the following paragraphs. Third and last, we collected a range of expert estimates and other data from various literature sources. This last category is the least dependable type of data, and we used it only after exhausting other possibilities. However, given the general and long-running scarcity of data on grassland fertilization91, there are in many cases no alternatives to expert estimates prior to the 1990s. For all these data sources, as far as possible, we followed the references to the original data source and collected the data from there. In each case the citations in this paper point to the publication from which we collected the data.All the data were collected in a table with columns for fertilizer rate R, fertilizer quantity Q, and crop area A12. Sometimes only one or two of these numbers were available. Each row in the table contains data from one publication concerning one country-year-crop combination.We classified each row as belonging to one of the land categories cropland (C), permanent grassland (PG), temporary grassland (TG), non-grass cropland (C–TG), total grassland (PG + TG), or in some cases only the fertilized portion of grasslands (denoted PGf, TGf, or PGf + TGf). The main point of the land categories is that they enable the calculation of the share applied to cropland, QC/(QC + QPG). However, they also helped with the intermediate tasks of interpretation and consistency checking, aggregation, and gap-filling:

    Consistency checking and interpretation. A significant complication in using the FAO/EFMA datasets87,88,89,90 is that the crop categorization for some countries includes temporary grassland under the “grassland” item and for other countries under another item, usually “fodder (other)”. Since the final goal was to separate permanent grassland from cropland, it was necessary to work out what each “grassland” item refers to. This issue could usually be resolved since these datasets also report crop areas corresponding to the different fertilizer rates. To assign the appropriate land category to such items, we cross-checked the fertilized areas according to the FAO/EFMA datasets against the areas according to FAOSTAT/Eurostat, and thus could in most cases unambigously determine whether the “grassland” fertilizer rates referred to permanent (PG) or total grassland (PG + TG). We used the FAO/EFMA fertilizer data only when the summed fertilizer quantities and areas matched the previously estimated country-level fertilizer quantities and land category areas.

    Aggregation. Several data sources, including national statistical surveys and the FAO/EFMA datasets, report data for individual crops. In these cases, we assigned land category labels (C, PG, C−TG, etc.) to the individual crops that together constitute such a land category. We then calculated, e.g., the average application rate on cropland RC as an area-weighted average of the individual crop rates. Where possible, we assigned labels C and PG since this allows the direct calculation of the final result QC/(QC + QPG). Several datasets, however, only distinguish total grassland (PG + TG) from non-grass cropland (C − TG); in these cases we later estimated permanent grassland rates from total grassland rates (see details below).

    Gap-filling. Several data sources only report rates, not areas or quantities. In such cases, when the land category was known, we used the previously collected area data as necessary to calculate quantities, e.g., QC = RCAC. Furthermore, when rates and areas were reported for fertilized grassland (e.g., RPGf and APGf), we calculated the average rate on grassland using grassland areas from FAOSTAT/Eurostat (e.g., RPG = RPGfAPGf/APG).

    Estimating the share of synthetic N fertilizer applied to croplandIn principle, the share of synthetic N applied to cropland can be calculated as QC/Qtot using a QC value from one of the above data sources and Qtot from the IFA/FAOSTAT databases. Alternatively, if only QPG is known, the quantity to cropland can be estimated as QC = Qtot − QPG. However, this method is not ideal since it introduces noise and possibly bias stemming from the different methods used to estimate total fertilizer consumption and land category fertilizer use. We therefore designed the following method to estimate the share to permanent grassland, which uses the total IFA/FAOSTAT fertilizer quantities only as a last resort.As a primary option, we used the following estimate. Assuming that all synthetic fertilizer is used on agricultural land (see Section “Synthetic N fertilizer consumption” above), the share applied to cropland can be expressed as$$frac{{Q}_{{rm{C}}}}{{Q}_{{rm{C}}}+{Q}_{{rm{PG}}}}=frac{{R}_{{rm{C}}}{A}_{{rm{C}}}}{{R}_{{rm{C}}}{A}_{{rm{C}}}+{R}_{{rm{PG}}}{A}_{{rm{PG}}}}={left(1+frac{{R}_{{rm{PG}}}}{{R}_{{rm{C}}}}frac{{A}_{{rm{PG}}}}{{A}_{{rm{C}}}}right)}^{-1}.$$
    (2)
    Note that this expression does not depend on the total fertilizer consumption (e.g., from IFA/FAOSTAT) but only on the area ratio APG/AC and the rate ratio RPG/RC. Whenever possible, we used this expression with area rate data from our own study and rate ratios calculated without mixing data from different publications. In the very few cases where different publications produced different rate ratio estimates for the same country/year, we used the average of the estimates.A variant of this primary option was used in the cases where rates are available only for total grassland (PG + TG) and non-grass cropland (C − TG). To estimate the rate ratio RPG/RC in these cases, additional data or assumptions are necessary. In some cases, we simply assumed that the average rate is equal on permanent and temporary grassland, but in other cases there was clear evidence of higher rates on temporary grassland. If the fertilizer rate on temporary grassland is k times the rate on permanent grassland, RTG = kRPG, then$${R}_{{rm{PG}}+{rm{TG}}}={R}_{{rm{PG}}}frac{{A}_{{rm{PG}}}+k{A}_{{rm{TG}}}}{{A}_{{rm{PG}}+{rm{TG}}}},$$
    (3)
    and similarly,$${R}_{{rm{C}}}=frac{{R}_{{rm{C}}-{rm{TG}}}{A}_{{rm{C}}-{rm{TG}}}+k{R}_{{rm{PG}}}{A}_{{rm{TG}}}}{{A}_{{rm{C}}}}.$$
    (4)
    It follows from Eqs. (3) and (4), using the shorthand m = APG+TG/(APG + kATG) and noting that RPG = mRPG+TG, that$$frac{{R}_{{rm{PG}}}}{{R}_{{rm{C}}}}=frac{m{R}_{{rm{PG}}+{rm{TG}}}{A}_{{rm{C}}}}{{R}_{{rm{C}}-{rm{TG}}}{A}_{{rm{C}}-{rm{TG}}}+km{R}_{{rm{PG}}+{rm{TG}}}{A}_{{rm{TG}}}}={left(frac{{R}_{{rm{C}}-{rm{TG}}}}{m{R}_{{rm{PG}}+{rm{TG}}}}frac{{A}_{{rm{C}}-{rm{TG}}}}{{A}_{{rm{C}}}}+kfrac{{A}_{{rm{TG}}}}{{A}_{{rm{C}}}}right)}^{-1}.$$
    (5)
    Note that this expression, like Eq. (2), does not depend on the total quantity Qtot but only on various rate ratios and area ratios.As a secondary and last option, if only one of the rates was known, we instead estimated the rate ratio using the IFA/FAOSTAT total fertilizer quantity Qtot. For example,$$frac{{R}_{{rm{PG}}}}{{R}_{{rm{C}}}}=frac{{R}_{{rm{PG}}}{A}_{{rm{C}}}}{{Q}_{{rm{C}}}}=frac{{R}_{{rm{PG}}}({A}_{{rm{tot}}}-{A}_{{rm{PG}}})}{{Q}_{{rm{tot}}}-{R}_{{rm{PG}}}{A}_{{rm{PG}}}},$$where Atot is the total agricultural land area. More generally, the rate ratio ({R}_{x}/{R}_{widetilde{x}}) between a land category x and its complement (widetilde{x}) (i.e., the remainder of the agricultural land) can be written as$$frac{{R}_{x}}{{R}_{widetilde{x}}}=frac{{R}_{x}({A}_{{rm{tot}}}-{A}_{x})}{{Q}_{{rm{tot}}}-{R}_{x}{A}_{x}}.$$
    (6)
    We used this equation in several cases (details below) but only when the rate ratio could not be calculated directly as the ratio of two rates.As a final step, for years without data, we gap-filled the estimated rate ratios RPG/RC using linear interpolation between data points and constant extrapolation before and after the first and last data point. We then estimated the share applied on cropland using Eq. (2) with area data (AC and APG) from FAOSTAT/Eurostat. We finally calculated the fertilizer quantities to cropland and permanent grassland using these shares so that the quantities agree with total fertilizer consumption statistics but not necessarily with the collected data on crop-specific rates and quantities. The resulting estimates of synthetic N application on cropland and permanent grassland are provided in the data record12.Countries not using synthetic N fertilizer on permanent grasslandBased on the FAO/EFMA data and other publications we concluded that the following countries have zero or negligible fertilizer rates to permanent grassland: Bulgaria, Croatia, Estonia, Finland, Greece, Hungary, Latvia, Lithuania, Portugal, Romania, Spain, and Sweden. These countries can conceptually be divided into two groups. One group is the countries in north-east Europe (Sweden, Finland, Estonia, Latvia, Lithuania) which all have substantial areas of temporary grassland. In Sweden and Finland, intensive grassland cultivation is a central part of agriculture but occurs almost exclusively on arable land92,93. In Estonia, Latvia, and Lithuania, the grassland cultivation practices have varied substantially through the decades; since 1992 much grassland has been abandoned or very extensively managed57,77,94,95, and the FAO/EFMA datasets87,90 suggest that if any synthetic N inputs have been applied to grassland, it has been to temporary grassland. The other group (Bulgaria, Croatia, Greece, Hungary, Portugal, Romania, and Spain) has small or negligible areas of temporary grassland but considerable areas of permanent grassland, which however are mostly extensively managed due to a combination of economic and climatic conditions17,77,96,97,98,99,100,101,102,103. For both groups, the quantitative data87,88,89,90 suggests that synthetic N inputs to permanent grassland was negligible at least during the 1990s and 2000s, and since no data could be found before the 1990s we have set it to zero for the whole period 1961–2019.A caveat to this assumption applies especially to the former communist states, which before the 1990s typically had much higher fertilizer inputs and different patterns of agricultural land use. In this study this is not an issue for Croatia, Estonia, Latvia, and Lithuania, which are included only from 1992, but for the former socialist republics of Bulgaria, Hungary, and Romania, which are included from 1961, we emphasize that data are very scarce.Countries using synthetic N fertilizer on permanent grasslandAustria. Most of Austria’s permanent grasslands are unfertilized104. Available quantitative data87,88,89,90,105 show that the rate ratio RPG/RC has been fairly constant in the period 1993–2006. Intensification of grassland fertilization is reported to have occured in the 1970s–1980s104 which roughly coincides with the steepest increase in total synthetic fertilizer use in Austria. We therefore extrapolated the average of the available rate ratio (RPG/RC ≈ 0.11) to the whole period 1961–2019.Belgium and Luxembourg. Belgium has a long history of heavy grassland fertilization106 but quantitative data are surprisingly scarce. The FAO/EFMA data shows that the average synthetic N rate on grassland (PG + TG) was on average 20% higher than on non-grass (C − TG) in the 1990s and 2000s. Since the fertilizer rates in Belgium-Luxembourg88 are practically indistinguishable from those referring to Belgium alone87,89,90, we pooled all the data and assumed the same rate ratio RPG/RC for Luxembourg ( More

  • in

    Intraspecific variation in thermal tolerance differs between tropical and temperate fishes

    1.Frölicher, T. L., Fischer, E. M. & Gruber, N. Marine heatwaves under global warming. Nature 560, 360–364 (2018).ADS 
    Article 

    Google Scholar 
    2.Stillman, J. H. Heat waves, the new normal: Summertime temperature extremes will impact animals, ecosystems, and human communities. Physiology 34, 86–100 (2019).CAS 
    Article 

    Google Scholar 
    3.Williams, C. M. et al. Biological impacts of thermal extremes: Mechanisms and costs of functional responses matter. Integr. Comp. Biol. 56, 73–84 (2016).Article 

    Google Scholar 
    4.Comte, L. & Olden, J. D. Climatic vulnerability of the world’s freshwater and marine fishes. Nat. Clim. Chang. 7, 718–722 (2017).ADS 
    Article 

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

    Google Scholar 
    6.Sunday, J. M., Bates, A. E. & Dulvy, N. K. Global analysis of thermal tolerance and latitude in ectotherms. Proc. R. Soc. B Biol. Sci. 278, 1823–1830 (2011).Article 

    Google Scholar 
    7.Rummer, J. L. et al. Life on the edge: Thermal optima for aerobic scope of equatorial reef fishes are close to current day temperatures. Glob. Chang. Biol. 20, 1055–1066 (2014).ADS 
    Article 

    Google Scholar 
    8.Stuart-Smith, R. D., Edgar, G. J., Barrett, N. S., Kininmonth, S. J. & Bates, A. E. Thermal biases and vulnerability to warming in the world’s marine fauna. Nature 528, 88–92 (2015).ADS 
    CAS 
    Article 

    Google Scholar 
    9.Pinsky, M. L., Eikeset, A. M., McCauley, D. J., Payne, J. L. & Sunday, J. M. Greater vulnerability to warming of marine versus terrestrial ectotherms. Nature 569, 108–111 (2019).ADS 
    CAS 
    Article 

    Google Scholar 
    10.Bennett, S., Duarte, C. M., Marbà, N. & Wernberg, T. Integrating within-species variation in thermal physiology into climate change ecology. Philos. Trans. R. Soc. B Biol. Sci. 374, 20180550 (2019).Article 

    Google Scholar 
    11.Comte, L. & Olden, J. D. Evolutionary and environmental determinants of freshwater fish thermal tolerance and plasticity. Glob. Chang. Biol. 23, 728–736 (2017).ADS 
    Article 

    Google Scholar 
    12.Bolnick, D. I. et al. Why intraspecific trait variation matters in community ecology. Trends Ecol. Evol. 26, 183–192 (2011).Article 

    Google Scholar 
    13.Pacifici, M. et al. Assessing species vulnerability to climate change. Nat. Clim. Chang. 5, 215–225 (2015).ADS 
    Article 

    Google Scholar 
    14.Moran, E. V., Hartig, F. & Bell, D. M. Intraspecific trait variation across scales: Implications for understanding global change responses. Glob. Change Biol. 22, 137–150 (2016).ADS 
    Article 

    Google Scholar 
    15.McKenzie, D. J. et al. Intraspecific variation in tolerance of warming in fishes. J. Fish Biol. 98, 1–20 (2020).

    Google Scholar 
    16.Mimura, M. et al. Understanding and monitoring the consequences of human impacts on intraspecific variation. Evol. Appl. 10, 121–139 (2017).Article 

    Google Scholar 
    17.Doyle, C. M., Leberg, P. L. & Klerks, P. L. Heritability of heat tolerance in a small livebearing fish, Heterandria formosa. Ecotoxicology 20, 535–542 (2011).CAS 
    Article 

    Google Scholar 
    18.Meffe, G. K., Weeks, S. C., Mulvey, M. & Kandl, K. L. Genetic differences in thermal tolerance of eastern mosquitofish (Gambusia holbrooki; Poeciliidae) from ambient and thermal ponds. Can. J. Fish. Aquat. Sci. 52, 2704–2711 (1995).Article 

    Google Scholar 
    19.Gradil, K. J., Garner, S. R., Wilson, C. C., Farrell, A. P. & Neff, B. D. Relationship between cardiac performance and environment across populations of Atlantic salmon (Salmo salar): A common garden experiment implicates local adaptation. Evol. Ecol. 30, 877–886 (2016).Article 

    Google Scholar 
    20.Fangue, N. A., Hofmeister, M. & Schulte, P. M. Intraspecific variation in thermal tolerance and heat shock protein gene expression in common killifish, Fundulus heteroclitus. J. Exp. Biol. 209, 2859–2872 (2006).CAS 
    Article 

    Google Scholar 
    21.Chen, Z., Farrell, A. P., Matala, A. & Narum, S. R. Mechanisms of thermal adaptation and evolutionary potential of conspecific populations to changing environments. Mol. Ecol. 27, 659–674 (2018).Article 

    Google Scholar 
    22.Chen, Z., Farrell, A. P., Matala, A., Hoffman, N. & Narum, S. R. Physiological and genomic signatures of evolutionary thermal adaptation in redband trout from extreme climates. Evol. Appl. 11, 1686–1699 (2018).CAS 
    Article 

    Google Scholar 
    23.Oliver, T. H. et al. Biodiversity and resilience of ecosystem functions. Trends Ecol. Evol. 30, 673–684 (2015).Article 

    Google Scholar 
    24.Mori, A. S., Furukawa, T. & Sasaki, T. Response diversity determines the resilience of ecosystems to environmental change. Biol. Rev. 88, 349–364 (2013).Article 

    Google Scholar 
    25.Rodgers, G. G., Donelson, J. M., McCormick, M. I. & Munday, P. L. In hot water: Sustained ocean warming reduces survival of a low-latitude coral reef fish. Mar. Biol. 165, 1–10 (2018).Article 

    Google Scholar 
    26.Seebacher, F., White, C. R. & Franklin, C. E. Physiological plasticity increases resilience of ectothermic animals to climate change. Nat. Clim. Chang. 5, 61–66 (2015).ADS 
    Article 

    Google Scholar 
    27.Meffe, G. K., Weeks, S. C., Mulvey, M. & Kandl, K. L. Genetic differences in thermal tolerance of eastern mosquitofish (Gambusia holbrooki; Poeciliidae) from ambient and thermal ponds. Can. J. Fish. Aquat. Sci. 52(12), 2704–2711 (1995).Article 

    Google Scholar 
    28.Baer, C. F. & Travis, J. Direct and correlated responses to artificial selection on acute thermal stress tolerance in a livebearing fish. Evolution 54(1), 238–244 (2000).CAS 
    Article 

    Google Scholar 
    29.Klerks, P. L., Athrey, G. N. & Leberg, P. L. Response to selection for increased heat tolerance in a small fish species, with the response decreased by a population bottleneck. Front. Ecol. Evol. 7, 270 (2019).Article 

    Google Scholar 
    30.Morgan, R., Finnøen, M. H., Jensen, H., Pélabon, C. & Jutfelt, F. Low potential for evolutionary rescue from climate change in a tropical fish. Proc. R. Soc. B Biol. Sci. 117(52), 33365–33372 (2020).CAS 

    Google Scholar 
    31.Frölicher, T. L. Extreme climatic events in the ocean. In Predicting Future Oceans: Sustainability of Ocean and Human Systems Amidst Global Environmental Change (eds Cisneros-Montemayor A. M. et al.) 53–60 (2019).32.Wernberg, T. et al. An extreme climatic event alters marine ecosystem structure in a global biodiversity hotspot. Nat. Clim. Chang. 3(1), 78–82 (2013).ADS 
    Article 

    Google Scholar 
    33.Burrows, M. T. et al. Geographical limits to species-range shifts are suggested by climate velocity. Nature 507(7493), 492–495 (2014).ADS 
    CAS 
    Article 

    Google Scholar 
    34.Molinos, J. G. et al. Climate velocity and the future global redistribution of marine biodiversity. Nat. Clim. Chang. 6(1), 83–88 (2016).ADS 
    Article 

    Google Scholar 
    35.Cox, D. K. Effects of three heating rates on the critical thermal maximum of bluegill. In Thermal Ecology (eds Gibbons, J. W. & Sharitz, R. R.) (National Technical Information Service, 1974).
    Google Scholar 
    36.Currie, S. & Schulte, P. M. Thermal stress. In The Physiology of Fishes 4th edn (eds Evans, D. H. et al.) 257–279 (CRC Press, 2014).
    Google Scholar 
    37.Grafen, A. The phlyogenetic regression. Philos. Trans. R. Soc. Lond. 326, 119–157 (1989).ADS 
    CAS 

    Google Scholar 
    38.Garland, T. Jr. & Ives, A. R. Using the past to predict the present: Confidence intervals for regression equations in phylogenetic comparative methods. Am. Nat. 155, 346–364 (2000).Article 

    Google Scholar 
    39.Orme, D. The caper package: comparative analysis of phylogenetics and evolution in R http://cran.r-project.org/web/packages/caper/index.html (2013).40.Hinchliff, C. E. et al. Synthesis of phylogeny and taxonomy into a comprehensive tree of life. Proc. Natl. Acad. Sci. 112(41), 12764–12769 (2015).ADS 
    CAS 
    Article 

    Google Scholar 
    41.Michonneau, F., Brown, J. W. & Winter, D. J. rotl: An R package to interact with the Open Tree of Life data. Methods Ecol. Evol. 7, 1476–1481 (2016).Article 

    Google Scholar 
    42.Freckleton, R. P. The seven deadly sins of comparative analysis. J. Evol. Biol. 22, 1367–1375 (2009).CAS 
    Article 

    Google Scholar  More

  • in

    Modeling characterization of the vertical and temporal variability of environmental DNA in the mesopelagic ocean

    To capture the first-order characteristics of the distribution of eDNA shed by mesopelagic species, we used a mechanistic model with one spatial dimension (the vertical dimension). It treats migrating organisms as a continuous point source of eDNA and simulates temporal evolution of the eDNA concentration vertical profiles in a one-dimensional ocean with climatological conditions. Horizontal variation of the distribution of eDNA shed by mesopelagic species is not considered in this study. The point source in the model represents a number of individuals within a single species. The simulation includes relevant eDNA transport processes using a vertical Price-Weller-Pinkel (PWP) dynamical module42. We use this model to first explore the sensitivity of eDNA concentration profiles to a variety of biological and physical parameters and to assess vertical and temporal variability of the eDNA concentration. Table 1 outlines the parameters, separated by type, and ranges of values for each parameter used in the sensitivity analysis (see “Sensitivity analysis” section below) and in a case study of how the model can be used to determine the percent of individuals that migrate (see “Application to ecological questions” section below). We then discuss how the mechanistic model can be used in conjunction with field sampling to test ecological hypotheses, highlighting how analysis of eDNA concentration profiles can provide information regarding where and when organisms are located in the water column. Finally, we comment on implications of this work and how the model can be used to both design and interpret observations.Table 1 Parameters separated by category.Full size tableModel overviewGoverning equationsTo maintain simplicity and also provide a more realistic scenario where both settling and neutrally buoyant eDNA particles28 are considered, we assume that the eDNA material at any given time in the system consist of eDNA particles of two size classes. Here, large eDNA particles are on the order of 10 s of µm to 1 mm and represent materials such as fecal pellets, chunks of tissue, or gametes that would be subject to both settling and physical breakdown over time. Small eDNA particles are on the order of 0.1 µm to 10 s of µm and would be any eDNA shed from an organism that has a density such that they are near-neutrally buoyant including extracellular eDNA or forms such as sperm, urine, blood, or single cells. This size class of small particles is so small that we can neglect their settling. In practice, these two size fractions of eDNA would be captured using the common method of filtering water through a < 1 µm pore size filter43. Both small and large particles in the model are subject to the same eDNA decay rate constant and large particles break down over time into small particles.The equation describing the change of eDNA concentration in the form of large particles is:$$begin{aligned} frac{partial {C_{LP}}}{partial {t}} = - w_vfrac{partial {C_{LP}}}{partial {z}} + frac{partial {}}{partial {z}}left( kappa _zfrac{partial {C_{LP}}}{partial {z}}right) - kC_{LP} - delta C_{LP} + w_sfrac{partial {C_{LP}}}{partial {z}} + hat{S_{LP}} end{aligned}$$ (1) where (C_{LP}) is the concentration of large eDNA particles, (w_v) is the vertical velocity [L/T] with a positive value pointing upward, (kappa _z) is the vertical diffusivity [L(^{2})/T], k is the first order decay rate constant [1/T], (delta) is the breakdown rate of particles [1/T], (w_s) is the settling rate of eDNA [L/T], and (hat{S_{LP}}) is the shedding rate of new large eDNA particles.Similarly, the equation describing the temporal evolution of the concentration of small particles of eDNA, (C_{SP}), is:$$begin{aligned} frac{partial {C_{SP}}}{partial {t}} = - w_vfrac{partial {C_{SP}}}{partial {z}} + frac{partial {}}{partial {z}}left( kappa _zfrac{partial {C_{SP}}}{partial {z}}right) - kC_{SP} + delta C_{LP} + hat{S_{SP}} end{aligned}$$ (2) where the parameters are the same as Eq. (1). Note that there is no settling term in Eq. (2) as small eDNA particles are assumed to be neutrally buoyant and that the (delta) term in Eq. (2) has the opposite sign as that in Eq. (1), reflecting eDNA volume conservation during the breaking-down of the large particles into small particles.Physical modelThe 1-D vertical model is implemented in MATLAB. The physical module is based on a previously published 1-D model of physical processes in the mixed layer42 and is modified to include the relevant organism and eDNA parameters. The model spans from the surface down to 1500 m deep with a vertical resolution of 0.5 m and has a time step of 10 s. Each simulation runs for a total of 90 days, representing one of the four seasons: summer (July, August, September - JAS), fall (October, November, December - OND), winter (January, February, March - JFM), and spring (April, May, June - AMJ). The model has a zero-gradient open boundary condition at the bottom. The model is forced on the surface by climatological, 8-hourly meteorological conditions from the NCEP/NCAR Reanalysis at a location in the Northeast Atlantic Ocean (37.047(^circ) N, -71.25(^circ) W). Inputs include air temperature, short wave radiation, long wave radiation, precipitation, air pressure, air humidity, and winds. The wind stress and heat fluxes are calculated using the bulk formulae44.The initial vertical distribution of the temperature and salinity in the model are seasonal climatological profiles obtained from the NCEI World Ocean Atlas 2018 at a site in the Northwest Atlantic slope sea (39.125(^circ) N, −70.875(^circ) W; north of the Gulf Stream). The model also simulates the influences of vertical mixing and advection on the temperature, salinity and eDNA particles. The mixing influence is incorporated by prescribing synthetic vertical profiles of vertical diffusivity that combines observed seasonal mean mixed layer depth45 and simulated seasonal mean vertical diffusivity profile from an operational model46, both at a site in the Northwest Atlantic slope sea (see Supplemental Text S1 and Supplemental Fig. S1). Vertical advection profiles are set to increase linearly from 0 ms(^{-1}) at the surface to the maximum value of ((10^{-4}) ms(^{-1})) at 200 m and then decrease linearly back to 0 ms(^{-1}) at 400 m (Supplemental Fig. S2). This prescribed advection profile does not change seasonally and represents enhanced vertical motions associated with sub-mesoscale processes in the upper water column (e.g.,47). Note that the result of this study is not sensitive to the prescribed vertical profiles of diffusivity or vertical velocity (see below).The physical parameters that will be adjusted in this study include: the initial temperature and salinity profiles, the mixed layer depth, the vertical diffusivity profile, and the vertical velocity profile (see Ocean parameters in Table 1).Organism movementThis study focuses on eDNA shed by vertically migrating organisms living in the mesopelagic ocean. For simplicity, we assume the organisms reside at the constant depth of 500 m during the day and the constant depth of 50 m during the night. These depths fall within the range of where migrating organisms are found at these times41,48. As the objective of this study is to develop a qualitative understanding of the vertical distribution of the eDNA concentration, the exact depth a particular mesopelagic species resides at in the daytime is not crucial. Organisms begin migrating to the surface two hours before sunset, and the upward migration ends one hour after sunset. The time of sunset and sunrise are determined seasonally using NOAA’s ESL solar calculator for the year 2019 at 42.35(^circ) N, −71.05(^circ) W. Both upward and downward migrations are assumed to be in constant speed, and the migration times do not change within each season (Fig. 2).Figure 2Representative organism migration curves for each season. Migration times change with sunrise and sunset for each season. This illustration assumes that 50% of individuals migrate. The shading shows the percentage of individuals with black indicating 100%, i.e., all organisms, and grey indicating 50%. These numbers vary among the simulations.Full size imageThe migration parameters that can be adjusted in the model include: daytime residing depth, nighttime residing depth, the start and end times of migrations (and thus duration), the percent of individuals that migrate, and the layer thickness (i.e., “width” of school). (See Organism parameters in Table 1). In this study, we focus on the start and end times of migrations (and thus duration) and the percent of individuals that migrate.eDNA parametersWe assume organisms are continuously shedding eDNA, and the shedding rate is fixed in time in each simulation. The eDNA then remains in the water column, where it is subject to transport and decay as described by Eqs. (1) and (2). Although several studies have characterized eDNA shedding rates of different marine organisms29,49,50, shedding rates are highly variable across studies and have high error associated with them, likely due to high temporal variability29. Here, we use a constant shedding rate of 1 mass unit per time step (10 seconds). Note that this study focuses on the relative vertical distribution of the eDNA concentration and its temporal variability, rather than its absolute value. The particle size distribution of eDNA (and the breakdown rate of large to small particles) and the settling rate of eDNA are largely unknown and expected to be time-varying28. Here we use a breakdown rate of fecal pellets to apply to large eDNA particles breaking down into small eDNA particles51. The most well characterized parameter of eDNA is decay rate, which several studies have characterized as a function of water temperature29,49. The upper and lower limits of the decay rate have also been well established29,31. Finally, there are currently no estimates of eDNA settling rates. We use values of marine snow settling rates52 to apply to large eDNA particles that are subject to settling in the model.The eDNA-related parameters that are adjusted in this mechanistic model include: settling rate of large particles; ratio of large particles to small particles that are shed by an organism; eDNA breakdown rate (large particles to small particles); the eDNA decay rate. (See eDNA parameters in Table 1). Note that shedding rates in the model do not change over the course of each simulation.Sensitivity analysisA series of 90-day simulations were carried out with altered values of the aforementioned parameters to examine the sensitivity of the eDNA vertical distribution to the parameters. Table 1 provides a list of the parameters that were adjusted in the sensitivity analysis. Note that the diffusivity profile (both the shape of the profile and the mixed layer depth), the temperature profile (which impacts the eDNA decay rate), and the daytime length (which impacted migration times) all change with the season. Other sensitivity parameters include vertical advection, settling rate of large particles, the ratio of large to small eDNA particles shed by organisms, and the percent of individuals that migrate. A total of 972 simulations were conducted with the following sets of parameters: 4 mixing profiles (one for each season), 3 vertical advection profiles (upwelling, downwelling, none), 3 settling rates of large eDNA particles, 3 decay rate constant scenarios (two constant values representing high and low values in the literature31, one temperature-dependent rate29 for each season, Supplemental Fig. S3), 3 scenarios for the percent of individuals that migrate, and 3 ratios of large to small eDNA particles shed by the organisms.In order to assess the impact of the parameters on the eDNA profiles, several metrics were defined. First, the water column was divided into three sections: surface layer (0–100 m), mid-depth (100–450 m), and deep water (450–550 m), based on our representative daytime and nighttime depths (500 m and 50 m) used for the simulations. For each simulation, the mean and maximum eDNA concentration in each depth bin was recorded. Also, the cumulative eDNA concentration in each depth bin was normalized to calculate the proportion of eDNA in each depth bin at any given time during the simulation. For each simulation, the mean and standard deviation of the proportion of eDNA in each depth bin were calculated over the course of the 90 day simulation.A first-order comparison of the vertical length scales of eDNA transport by different processes were conducted. Here, we use a time scale of eDNA decay, T(_{90}), i.e., the time it takes for 90% of the released eDNA to decay, to determine the vertical length scale of each process. In particular, we would like to estimate, the vertical distances eDNA is transported by advection, mixing, and settling before 90% of the released eDNA has decayed, L(_{mix}), L(_{advect}) and L(_{settling}), respectively. These quantities are defined mathematically as,$$begin{aligned} T_{90} = frac{-ln(0.1)}{k_{avg}} end{aligned}$$ (3) $$begin{aligned} L_{mix} = sqrt{kappa _v T_{90}} end{aligned}$$ (4) $$begin{aligned} L_{advect} = w_{vm} T_{90} end{aligned}$$ (5) $$begin{aligned} L_{settle} = w_s T_{90} end{aligned}$$ (6) where (k_{avg}) is the average decay rate constant [T(^{-1})] for the whole water column over the 90 day simulation, (kappa _v) is the maximum vertical diffusivity coefficient [L(^{2})T(^{-1})] , (w_{vm}) is the maximum vertical velocity [LT(^{-1})] , and (w_s) is the settling rate [LT(^{-1})]. Note that because k is temperature dependent, the value of the decay rate constant will change both vertically (e.g., deeper water will be colder and have a lower decay rate constant) and temporally (e.g., at a given depth, the water temperature will be warmer in the middle of the day and thus will have a higher decay rate constant).Application to ecological questionsAfter obtaining a first-order understanding of the eDNA concentration profile and the impact of the parameters on the eDNA distribution, we ran another set of model simulations to examine if vertical and temporal eDNA concentration variability can shed light on an ecological question regarding what percentage of individuals within the same species migrate on a daily basis. To do this, for each season, all parameters were held constant except for the percent of individuals that migrate, P(_{m}). In each season, the value of P(_{m}) varies from 0 to 100% with an increment of 10%. A total of 44 such simulations (4 seasons x 11 percent migrate scenarios) were carried out. We then compared mean and maximum eDNA concentrations in the surface layer (0–100 m) to those in the deep water (450 to 550 m), and examined the relationship between the surface-to-deep ratios and P(_{m}). More

  • in

    Impact of intensifying nitrogen limitation on ocean net primary production is fingerprinted by nitrogen isotopes

    Modelling approachWe used the PISCES-v2 biogeochemical model, attached to the Nucleus for European Modelling of the Ocean version 4.0 (NEMO-v4) general ocean circulation model29. PISCES-v2 includes five nutrients pools (nitrate, ammonium, phosphate, silicic acid and dissolved iron), dissolved oxygen, the full carbon system and accounts for two phytoplankton (nanophytoplankton and diatoms) and two zooplankton types (microzooplankton and mesozooplankton). Bioavailable nitrogen in our simulations is considered to be the combination of nitrate and ammonium. Its nitrogen cycle includes nitrogen fixation, nitrification, burial, denitrification in both the water column and sediments, and coupled nitrification–denitrification. Nitrogen isotopes were integrated within PISCES-v2 for the purposes of this study, using nine new tracers (Supplementary Note 1). Horizontal model resolution varied between ~0.5° at the equator and poles, and 2° in the subtropics, whereas vertical resolution varied between 10 and 500 m thickness over 31 levels.We conducted simulations under both preindustrial control and climate change scenarios. The preindustrial control scenario from 1801 to 2100 maintained preindustrial greenhouse gas concentrations and only included internal modes of variability. The climate change simulation from 1851 to 2100 included natural variability, prescribed changes in land use, as well as historical changes in concentrations of greenhouse gases and aerosols until 2005, after which future concentrations associated with RCP8.5 were imposed30. The biogeochemical model (PISCES-v2) was run offline from the physical model (NEMO-v4) using monthly transports and other physical conditions generated by the low resolution version of the IPSL-CM5A ESM57.Experiments were initialized from biogeochemical fields created from an extensive spin-up of 5000 years under repeat physical forcing, followed by a 300-year simulation under the preindustrial control scenario. The preindustrial control simulation used in analysis was therefore the final 300 years of a 5600-year spin-up involving two repeat simulations of the preindustrial control scenario. We utilized a global compilation of δ15NNO320 supplemented with recent data to assess the isotopic routines in the model and conducted a thorough model-data skill assessment at replicating observed patterns in space (Supplementary Note 2 and Supplementary Figs. 1–3).Anthropogenic nitrogen depositionThe effect of increasing aeolian deposition of nitrogen was assessed in our simulations. Preindustrial nitrogen deposition was prescribed as the preindustrial estimate at 1850, whereas the historical to future deposition was created by linear interpolation between preindustrial (1850) and modern/future fields (2000, 2030, 2050 and 2100). These fields were provided by Hauglustaine et al.8. However, the rapid rise between 1950 and 2000 was maintained, such that 60% of the increase between the preindustrial and modern fields occurred after 1950 (Supplementary Fig. 4).The historical rise in anthropogenic nitrogen deposition was assessed by including it in additional simulations under both preindustrial control and climate change scenarios. Four initial experiments were therefore conducted: preindustrial control; preindustrial control plus anthropogenic nitrogen deposition; climate change; and climate change plus anthropogenic nitrogen deposition.Global model experimentsWe undertook four initial simulations to quantify the impacts of anthropogenic climate change and nitrogen deposition: a preindustrial control simulation from 1801 to 2100; a full anthropogenic scenario from 1851 to 2100; a climate change-only scenario without the increase in anthropogenic nitrogen deposition from 1851 to 2100; and a nitrogen deposition scenario without anthropogenic climate change from 1851 to 2100. Anthropogenic effects to nitrogen cycling were quantified by comparing mean conditions over the final 20 years of the twenty-first century (2081–2100) with mean conditions over the final 20 years of the preindustrial control simulation, whereas effects on nitrogen isotopes were quantified by comparing mean conditions over the final 20 years of the twenty-first century (2081–2100) with mean conditions over the historical period (1986–2005) from the same simulation.To understand the direct and indirect effects of climate change, we undertook two additional idealized simulations. First, we imposed temperature changes on biogeochemical rates, while maintaining ocean circulation associated with the preindustrial control scenario, to assess the direct effects of warming on biogeochemical processes. Second, we imposed the preindustrial control temperature field on biogeochemical processes, while altering the circulation in line with the climate change scenario, to assess the indirect effects of climate change (i.e., how changing circulation alters substrate supply to biogeochemical reactions). Each experiment was run from 1851 to 2100 and without the anthropogenic increase in atmospheric nitrogen deposition, parallel with the full climate change simulation.Agreement between the climate change simulation without anthropogenic nitrogen deposition was quantified using a pixel-by-pixel correlation analysis using Spearman’s rank correlation based on the non-parametric nature of the two-dimensional fields used for comparison. Fields were euphotic zone nitrate, twilight zone δ15NNO3, euphotic zone δ15NPOM, and vertically integrated NPP, zooplankton grazing, nitrogen fixation, water column denitrification and sedimentary denitrification.Depth zonesWe assessed changes in biogeochemical variables related to nitrogen cycling in two depth zones defined by light. The euphotic zone was defined by depths between the surface and 0.1% of incident irradiance as recommended by Buesseler et al.42. The twilight zone was also defined using light, as advocated by Kaartvedt et al.58. Depths between 0.1% and 0.0001% of incident irradiance defined the twilight zone. These definitions typically returned euphotic zone thicknesses of 137 ± 23 m (mean ± SD), and twilight zone thicknesses of 233 ± 37 m. The boundary between these depth zones were deepest in oligotrophic tropical and subtropical waters, and were shallowest in equatorial and temperate waters (Supplementary Fig. 7).Time of emergenceToE calculations determined when anthropogenic, anomalous trends emerged from the noise of background variability. ToE was calculated at each grid cell within both the euphotic and twilight zones (depth-averaged) and using annually averaged fields of ocean tracers. We therefore ignored temporal trends and variability at seasonal and sub-seasonal scales. Raw time series were first detrended and normalized using the linear slope and mean of the preindustrial control experiment, such that the preindustrial control time series varied about zero, while anomalous trends in experiments with climate change and/or nitrogen deposition deviated from zero. These detrended and normalized time series were smoothed using a boxcar (flat) moving average with a window of 11 years to filter decadal variability (Supplementary Fig. 12). Differences with the preindustrial control experiment were then computed.To determine whether the differences with the preindustrial control experiment were anomalous, we calculated a measure of noise from the raw, inter-annual time series of the preindustrial control experiment (1801–2100). A signal emerged from the noise if it exceeded 2 SDs, a threshold that represents with 95% confidence that a value was anomalous and is therefore a conservative envelope to distinguish normality from anomaly16.Furthermore, we required that anomalous values must consistently exceed the noise of the preindustrial control experiment until the end of the simulation (2100) to be registered as having emerged. Temporary emergences were therefore rejected, making our ToE estimates more conservative. A graphical representation of this process is shown in Supplementary Fig. 12.Isolating biogeochemical 15NO3 fluxesWe analysed the biogeochemical fluxes of 15NO3 and NO3 into and out of each model grid cell within the twilight zone, to determine whether the trends in δ15NNO3 were related to biogeochemical or physical changes. Fluxes of 15NO3 and NO3 included a net source from nitrification (NO3nitr) and net sinks due to new production (NO3new) and denitrification (NO3den). Although nitrification did not directly alter the 15N : 14N ratio in our simulations, the release of 15NO3 and NO3 by nitrification conveyed an isotopic signature determined by prior fractionation processes that produce ammonium (NH4). These processes include remineralization of particulate and dissolved organic matter, excretion by zooplankton and nitrogen fixation. The isotopic signatures of these processes were thus included implicitly in NO3nitr. For each grid cell, we calculated the biogeochemical tendency to alter δ15NNO3 based on the ratio of inputs minus outputs:$${Delta} {delta }^{15}{{{{{{rm{N}}}}}}}_{{{{{{rm{NO3}}}}}}}=left(frac{{,{!}^{15}{{{{{rm{N}}}}}}{{{{{rm{O}}}}}}}_{3}^{{{{{{rm{nitr}}}}}}}-{,{!}^{15}{{{{{rm{N}}}}}}{{{{{rm{O}}}}}}}_{3}^{{{{{{rm{new}}}}}}}-{,{!}^{15}{{{{{rm{N}}}}}}{{{{{rm{O}}}}}}}_{3}^{{{{{{rm{den}}}}}}}}{{,{!}^{14}{{{{{rm{N}}}}}}{{{{{rm{O}}}}}}}_{3}^{{{{{{rm{nitr}}}}}}}-{,{!}^{14}{{{{{rm{N}}}}}}{{{{{rm{O}}}}}}}_{3}^{{{{{{rm{new}}}}}}}-{,{!}^{14}{{{{{rm{N}}}}}}{{{{{rm{O}}}}}}}_{3}^{{{{{{rm{den}}}}}}}}-1right)cdot 1000$$
    (1)
    This calculation excluded any upstream biological changes and circulation changes that might have altered δ15NNO3.0D water parcel modelWe simulated the nitrogen isotope dynamics in a recently upwelled water parcel during transit to the subtropics by building a 0D model. The model simulates state variables of dissolved inorganic nitrogen (DIN), particulate organic nitrogen (PON) and exported particulate nitrogen (ExpN), as well as their heavy isotopes (DI15N, PO15N and Exp15N) in units of mmol N m−3 over 100 days given initial conditions and constants listed in Supplementary Table 1.$$frac{Delta {{{{{rm{DIN}}}}}}}{Delta t}=-{{{{{mathrm{N}}}}}}_{{{{{{rm{uptake}}}}}}}+{{{{{mathrm{N}}}}}}_{{{{{{rm{recycled}}}}}}}$$
    (2)
    $$frac{Delta {{{{{rm{PON}}}}}}}{Delta t}={{{{{mathrm{N}}}}}}_{{{{{{rm{uptake}}}}}}}-{{{{{mathrm{N}}}}}}_{{{{{{rm{recycled}}}}}}}-{{{{{mathrm{N}}}}}}_{{{{{{rm{exported}}}}}}}$$
    (3)
    $$frac{Delta {{{{{rm{ExpN}}}}}}}{Delta t}={{{{{mathrm{N}}}}}}_{{{{{{rm{exported}}}}}}}$$
    (4)
    $$frac{Delta {{{{{rm{DI1}}}}}}{}^{15}{{{{{rm{N}}}}}}}{Delta t}=-{}^{15}{{{{{rm{N}}}}}}_{{{{{{rm{uptake}}}}}}}+{}^{15}{{{{{rm{N}}}}}}_{{{{{{rm{recycled}}}}}}}$$
    (5)
    $$frac{Delta {{{{{rm{PO}}}}}}{}^{15}{{{{{rm{N}}}}}}}{Delta t}={}^{15}{{{{{rm{N}}}}}}_{{{{{{rm{uptake}}}}}}}-{}^{15}{{{{{rm{N}}}}}}_{{{{{{rm{recycled}}}}}}}-{}^{15}{{{{{rm{N}}}}}}_{{{{{{rm{exported}}}}}}}$$
    (6)
    $$frac{Delta {{{{mathrm{Exp}}}}}{}^{15}{{{{{rm{N}}}}}}}{Delta t}={}^{15}{{{{{rm{N}}}}}}_{{{{{{rm{exported}}}}}}}$$
    (7)
    First, the model calculates maximum potential growth rate of phytoplankton (μmax) in units of day−1 (Eq. 8) using temperature and then finds nitrogen uptake (Nuptake, Eq. 10) using PON and limitation terms for nitrogen (Nlim, Eq. 9), light (Llim, Supplementary Table 1) and iron (Felim, Supplementary Table 1).$${mu }_{{{max }}}=0.6,{{{{{rm{da}}}}}}{y}^{-1}cdot {e}^{Tcdot {T}_{{{{{{rm{growth}}}}}}}}$$
    (8)
    $${{{{{mathrm{N}}}}}}_{{{{{mathrm{lim}}}}}}=frac{{{{{{rm{DIN}}}}}}}{{{{{{rm{DIN}}}}}}+{{{{{mathrm{K}}}}}}_{{{{{{rm{DIN}}}}}}}}$$
    (9)
    $${{{{{mathrm{N}}}}}}_{{{{{mathrm{uptake}}}}}}={mu }_{max }cdot {{{{{mathrm{L}}}}}}_{{{{{mathrm{lim}}}}}}cdot ,min ({{{{{mathrm{Fe}}}}}}_{{{{{mathrm{lim}}}}}},{{{{{mathrm{N}}}}}}_{{{{{mathrm{lim}}}}}})cdot {{{{{mathrm{PON}}}}}}$$
    (10)
    At a constant temperature of 18 °C, μmax is equal to ~1.9 day−1. Limitation terms for light and iron are set as constant and are used to prevent unrealistically high nitrogen uptake when nitrogen is high, such as occurs immediately following upwelling in the high-nutrient low-chlorophyll regions of the tropics. Fractionation by phytoplankton is calculated assuming an open system21, in this case where nitrogen can be lost through export of organic matter. To calculate the fractionation associated with uptake (15Nuptake, Eq. 11), we multiply the total nitrogen uptake (Nuptake, Eq. 10) by the heavy to light isotope ratio (({r}_{{{{{{rm{DIN}}}}}}}^{15}), Eq. 12) and the fractionation factor (εphy, Supplementary Table 1), which is converted from units of per mil (‰) to a fraction relative to one. This fractionation factor (εphy) is constant at 5‰ but is decreased towards 0‰ by the nitrogen limitation term (Nlim, Eq. 9), such that when nitrogen is limiting to growth, the fractionation during uptake decreases (last term on the right-hand side approaches 1).$${}^{15}{{{{{rm{N}}}}}}_{{{{{{rm{uptake}}}}}}},=,{{{{{mathrm{N}}}}}}_{{{{{{rm{uptake}}}}}}}cdot {r}_{{{{{{rm{DIN}}}}}}}^{15}cdot left(1-frac{{{{{mathrm{N}}}}}_{{{{{mathrm{lim}}}}}}cdot {varepsilon }_{{{{{{rm{phy}}}}}}}}{1000}right)$$
    (11)
    $${r}_{{{{{{rm{DIN}}}}}}}^{15},=,frac{{{{mathrm{DI}}}}^{15}{{{{{{rm{N}}}}}}}}{{{{{{rm{DIN}}}}}}}$$
    (12)
    At each timestep, a fraction of the PON pool becomes detritus (Eq. 15) and this detritus is instantaneously recycled back to DIN or exported to ExpN and removed from the water parcel. The amount of detritus produced per timestep is calculated as the sum of linear respiration (Eq. 13) and quadratic mortality (Eq. 14) terms, where Presp (units of day−1), Kresp (units of mmol N m−3) and Pmort (units of (mmol N m−3)−1 day−1) are constants (Supplementary Table 1).$${{{{{rm{Respiration}}}}}},=,{{{{{mathrm{P}}}}}}_{{{{{{rm{resp}}}}}}}cdot {{{{{rm{PON}}}}}}cdot frac{{{{{{rm{PON}}}}}}}{{{{{{rm{PON}}}}}}+{{{{{mathrm{K}}}}}}_{{{{{{rm{resp}}}}}}}}$$
    (13)
    $${{{{{rm{Mortality}}}}}},=,{{{{{mathrm{P}}}}}}_{{{{{{rm{mort}}}}}}}cdot {{{{{rm{PON}}}}}}^{2}$$
    (14)
    $${{{{{rm{Detritus}}}}}},=,{{{{{rm{Respiration}}}}}},+,{{{{{rm{Mortality}}}}}}$$
    (15)
    Once we know the fraction of PON that becomes detritus at any given timestep, we must solve for the fraction of that detritus that becomes DIN through recycling (Eq. 17), and that which becomes ExpN through export (Eq. 18). The fraction of detritus that is recycled back into DIN is temperature dependent (Eq. 16), with higher temperatures increasing rates of recycling above a minimum fraction set by frecmin (Supplementary Table 1). The relationship with temperature is exponential, similar to phytoplankton maximum growth (μmax), but the degree of increase associated with warming is scaled down by a constant factor equal to Trec (Supplementary Table 1). The fraction that is exported to ExpN is the remainder (Eq. 18).$${f}_{{{{{{rm{recycled}}}}}}}={f}_{{{{{{rm{recmin}}}}}}}+{T}_{{{{{{rm{rec}}}}}}}cdot {e}^{Tcdot {T}_{{{{{{rm{growth}}}}}}}}$$
    (16)
    $${{{{{mathrm{N}}}}}}_{{{{{{rm{recycled}}}}}}}={{{{{rm{Detritus}}}}}}cdot {f}_{{{{{{rm{recycled}}}}}}}$$
    (17)
    $${{{{{mathrm{N}}}}}}_{{{{{{rm{exported}}}}}}}={{{{{rm{Detritus}}}}}}cdot (1-{f}_{{{{{{rm{recycled}}}}}}})$$
    (18)
    The major fluxes of Nuptake, Nrecycled and Nexported are now solved for. All that remains is to calculate the isotopic signatures of the recycling (Eq. 19) and export (Eq. 20) fluxes. These, similar to 15Nuptake (Eq. 11), are solved by multiplying against a standard ratio of heavy to light isotope (({r}_{{{{{{rm{PON}}}}}}}^{15}), Eq. 21).$${}^{15}{{{{{rm{N}}}}}}_{{{{{{rm{recycled}}}}}}}={{{{{mathrm{N}}}}}}_{{{{{{rm{recycled}}}}}}}cdot {r}_{{{{{{rm{PON}}}}}}}^{15}$$
    (19)
    $${}^{15}{{{{{rm{N}}}}}}_{{{{{{rm{exported}}}}}}}={{{{{mathrm{N}}}}}}_{{{{{{rm{exported}}}}}}}cdot {r}_{{{{{{rm{PON}}}}}}}^{15}$$
    (20)
    $${r}_{{{{{{rm{PON}}}}}}}^{15}=frac{{{{{{rm{PO}}}}}}{}^{15}{{{{{rm{N}}}}}}}{{{{{{rm{PON}}}}}}}$$
    (21)
    Finally, we calculate the δ15N values of the major pools in the model (DIN, PON and ExpN) as output (Eqs. 22–24). We assume in this model that the major pools of DIN, PON and ExpN represent the total amount of the light isotope (14N), whereas the DI15N, PO15N and Exp15N pools represent the relative enrichment in 15N compared to a standard ratio. For simplicity, we make the standard ratio equal to 1. Therefore, taking the ratio of the DI15N to DIN pools and subtracting one returns the isotopic signature. Multiplying this by 1000 converts this signature to per mil units (‰).$${delta }^{15}{{{{{{rm{N}}}}}}}_{{{{{{rm{DIN}}}}}}}=left(frac{{{{{{rm{DI}}}}}}{}^{15}{{{{{rm{N}}}}}}}{{{{{{rm{DIN}}}}}}}-1right)cdot 1000$$
    (22)
    $${}^{15}{{{{{rm{N}}}}}}_{{{{{{rm{PON}}}}}}}=left(frac{{{{{{{rm{PO}}}}}}}^{15}N}{{{{{{rm{PON}}}}}}}-1right)cdot 1000$$
    (23)
    $${}^{15}{{{{{rm{N}}}}}}_{{{{{{rm{ExpN}}}}}}}=left(frac{{{{{mathrm{Exp}}}}}{}^{15}{{{{{rm{N}}}}}}}{{{{{{rm{ExpN}}}}}}}-1right)cdot 1000$$
    (24) More

  • in

    Protected areas are not effective for the conservation of freshwater insects in Brazil

    1.Brooks, T. M. et al. Global biodiversity conservation priorities. Science (80-. ). 313, 58–61 (2006).2.Camacho-Sandoval, J. & Duque, H. Indicators for biodiversity assessment in Costa Rica. Agric. Ecosyst. Environ. 87, 141–150 (2001).Article 

    Google Scholar 
    3.Diniz-Filho, J. A. F. et al. Ensemble forecasting shifts in climatically suitable areas for Tropidacris cristata (Orthoptera: Acridoidea: Romaleidae). Insect Conserv. Divers. https://doi.org/10.1111/j.1752-4598.2010.00090.x (2010).Article 

    Google Scholar 
    4.Morse-Jones, S. et al. Stated preferences for tropical wildlife conservation amongst distant beneficiaries: Charisma, endemism, scope and substitution effects. Ecol. Econ. 78, (2012).5.Verissimo, D., MacMillan, D. C. & Smith, R. J. Toward a systematic approach for identifying conservation flagships. Conserv. Lett. vol. 4 (2011).6.Nóbrega, C. C. & De Marco, P. Unprotecting the rare species: a niche-based gap analysis for odonates in a core Cerrado area. Divers. Distrib. 17, 491–505 (2011).Article 

    Google Scholar 
    7.SNUC, (Sistema Nacional de Unidades de Conservação da Natureza). Lei no 9.985, de 18 de julho de 2000. Mma/Sbf (2000) doi:https://doi.org/10.1017/CBO9781107415324.004.8.Abell, R., Allan, J. D. & Lehner, B. Unlocking the potential of protected areas for freshwaters. Biol. Conserv. 134, 48–63 (2007).Article 

    Google Scholar 
    9.Monteiro, C. da S., Esposito, M. C. & Juen, L. Are the adult odonate species found in a protected area different from those present in the surrounding zone? A case study from eastern Amazonia. J. Insect Conserv. 20, 643–652 (2016).10.Margules, C. R. & Pressey, R. L. Systematic conservation planning. Nature 405, 243–253 (2000).CAS 
    Article 

    Google Scholar 
    11.Whittaker, R. J. et al. Conservation biogeography: assessment and prospect. Divers. Distrib. 11, 3–23 (2005).Article 

    Google Scholar 
    12.Bini, L. M., Diniz-Filho, J. A. F., Rangel, T. F. L. V. B., Bastos, R. P. & Pinto, M. P. Challenging Wallacean and Linnean shortfalls: knowledge gradients and conservation planning in a biodiversity hotspot. Divers. Distrib. https://doi.org/10.1111/j.1366-9516.2006.00286.x (2006).Article 

    Google Scholar 
    13.Rodrigues, A. S. L. & Gaston, K. J. Maximising phylogenetic diversity in the selection of networks of conservation areas. Biol. Conserv. https://doi.org/10.1016/S0006-3207(01)00208-7 (2002).Article 

    Google Scholar 
    14.Silva, D. C., Vieira, T. B., da Silva, J. M. & de Cassia Faria, K. Biogeography and priority areas for the conservation of bats in the Brazilian Cerrado. Biodivers. Conserv. 27, 815–828 (2018).15.Salkeld, D. J., Padgett, K. A. & Jones, J. H. A meta-analysis suggesting that the relationship between biodiversity and risk of zoonotic pathogen transmission is idiosyncratic. Ecol. Lett. 16, 679–686 (2013).Article 

    Google Scholar 
    16.Juen, L. & de Marco, P. Dragonfly endemism in the Brazilian Amazon: competing hypotheses for biogeographical patterns. Biodivers. Conserv. https://doi.org/10.1007/s10531-012-0377-0 (2012).Article 

    Google Scholar 
    17.Mendes, S. L. et al. Protected Areas for the Northern Muriqui, Brachyteles hypoxanthus (Primates, Atelidae). Neotrop. Primates 13, (2005).18.Serra, B. D. V., De Marco Júnior, P., Nóbrega, C. C. & Campos, L. A. D. O. Modeling potential geographical distribution of the wild nests of Melipona capixaba Moure & Camargo, 1994 (Hymenoptera, apidae): conserving isolated populations in mountain habitats. Nat. a Conserv. 10, 199–206 (2012).19.Mendes, P. & De Marco, P. Bat species vulnerability in Cerrado: integrating climatic suitability with sensitivity to land-use changes. Environ. Conserv. 45, 67–74 (2018).Article 

    Google Scholar 
    20.Brasil, L. S. et al. A niche‐based gap analysis for the conservation of odonate species in the Brazilian Amazon. Aquat. Conserv. Mar. Freshw. Ecosyst. aqc.3599 (2021) doi:https://doi.org/10.1002/aqc.3599.21.da Silva, J. G., Vieira, T. B. & Mews, H. A. Fine-scale effect of environmental variation and distance from watercourses on pteridophyte assemblage structure in the western Amazon. Folia Geobot. https://doi.org/10.1007/s12224-021-09390-y (2021).Article 

    Google Scholar 
    22.Doughty, C. R. Freshwater biomonitoring and benthic macroinvertebrates, edited by D. M. Rosenberg and V. H. Resh, Chapman and Hall, New York, 1993. ix + 488pp. ISBN 0412 02251 6. Aquat. Conserv. Mar. Freshw. Ecosyst. 4, 92–92 (1994).23.Harper, D. M., Rosenberg, D. A. & Resh, V. H. Freshwater biomonitoring and benthic macroinvertebrates. J. Appl. Ecol. 31, 790 (1994).Article 

    Google Scholar 
    24.Cunha, E. J. & Juen, L. Impacts of oil palm plantations on changes in environmental heterogeneity and Heteroptera (Gerromorpha and Nepomorpha) diversity. J. Insect Conserv. 21, 111–119 (2017).Article 

    Google Scholar 
    25.Schuh, R. T. & Slater, J. A. True bugs of the World (Hemiptera: Heteroptera). Classification and Natural History. (Cornell University Press, 1995).26.Giehl, N. F. da S., Dias-Silva, K., Juen, L., Batista, J. D. & Cabette, H. S. R. Taxonomic and Numerical Resolutions of Nepomorpha (Insecta: Heteroptera) in Cerrado Streams. PLoS One 9, e103623 (2014).27.Dias-Silva, K., Cabette, H. S. R., Juen, L. & Jr, P. D. M. The influence of habitat integrity and physical-chemical water variables on the structure of aquatic and semi-aquatic Heteroptera. Zool. 27, 918–930 (2010).28.Panizzi, A. R. & Grazia, J. True Bugs (Heteroptera) of the Neotropics. True Bugs (Heteroptera) of the Neotropics vol. 2 (Springer Netherlands, 2015).29.Polhemus, J. T. & Polhemus, D. A. Global diversity of true bugs (Heteroptera; Insecta) in freshwater. Hydrobiologia https://doi.org/10.1007/s10750-007-9033-1 (2008).Article 

    Google Scholar 
    30.Nieser, N. & Melo, A. L. Os Heterópteros Aquáticos de Minas Gerais. (UFMG, Belo Horizonte, 1997).31.Cunha, E. J., de Assis Montag, L. F. & Juen, L. Oil palm crops effects on environmental integrity of Amazonian streams and Heteropteran (Hemiptera) species diversity. Ecol. Indic. 52, 422–429 (2015).32.Cordeiro, I. & Moreira, F. New distributional data on aquatic and semiaquatic bugs (Hemiptera: Heteroptera: Gerromorpha & Nepomorpha) from South America. Biodivers. Data J. 3, e4913 (2015).33.Rodrigues, A. S. L. & Brooks, T. M. Shortcuts for biodiversity conservation planning: the effectiveness of surrogates. Annu. Rev. Ecol. Evol. Syst. 38, 713–737 (2007).Article 

    Google Scholar 
    34.Andelman, S. J. & Fagan, W. F. Umbrellas and flagships: Efficient conservation surrogates or expensive mistakes?. Proc. Natl. Acad. Sci. 97, 5954–5959 (2000).ADS 
    CAS 
    Article 

    Google Scholar 
    35.Fielding, A. H. & Bell, J. F. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ. Conserv. 24, 38–49 (1997).Article 

    Google Scholar 
    36.Abellan, P., Sanchez-Fernandez, D., Velasco, J. & Millan, A. Conservation of freshwater biodiversity: a comparison of different area selection methods. Biodivers. Conserv. 14, 3457–3474 (2005).Article 

    Google Scholar 
    37.Fearnside, P. M. Conservation policy in brazilian amazonia: understanding the dilemmas. World Dev. 31, 757–779 (2003).Article 

    Google Scholar 
    38.dos Santos, A. J., Vieira, T. B. & Faria, K. de C. Effects of vegetation structure on the diversity of bats in remnants of Brazilian Cerrado savanna. Basic Appl. Ecol. 17, 720–730 (2016).39.Groves, C. R. et al. Planning for biodiversity conservation: putting conservation science into practice. Bioscience https://doi.org/10.1641/0006-3568(2002)052[0499:pfbcpc]2.0.co;2 (2002).Article 

    Google Scholar 
    40.Fearnside, P. M. & Ferraz, J. A conservation gap analysis of Brazil’s Amazonian vegetation. Conserv. Biol. 9, 1134–1147 (1995).Article 

    Google Scholar 
    41.Fearnside, P. M. Introduction: strategies for social and environmental conservation in conservation units. In The Amazon Várzea 233–238 (Springer Netherlands, 2011). doi:https://doi.org/10.1007/978-94-007-0146-5_16.42.Cardoso, P., Erwin, T. L., Borges, P. A. V. & New, T. R. The seven impediments in invertebrate conservation and how to overcome them. Biol. Conserv. 144, 2647–2655 (2011).Article 

    Google Scholar 
    43.Marini, M. Â. & Garcia, F. I. Bird conservation in Brazil. Conserv. Biol. https://doi.org/10.1111/j.1523-1739.2005.00706.x (2005).Article 

    Google Scholar 
    44.Young, B. E. et al. Population declines and priorities for amphibian conservation in Latin America. Conserv. Biol. 15, 1213–1223 (2001).Article 

    Google Scholar 
    45.Dias-Silva, K., Moreira, F. F. F., Giehl, N. F. D. S., Nóbrega, C. C. & Cabette, H. S. R. Gerromorpha (Hemiptera: Heteroptera) of eastern Mato Grosso State, Brazil: checklist, new records, and species distribution modeling. Zootaxa https://doi.org/10.11646/zootaxa.3736.3.1 (2013).Article 
    PubMed 

    Google Scholar 
    46.Ferraz, K. M. P. M. de B., Ferraz, S. F. de B., Paula, R. C. de, Beisiegel, B. & Breitenmoser, C. Species Distribution Modeling for Conservation Purposes. Nat. Conserv. 10, 214–220 (2012).47.Marco-Júnior, P. & Siqueira, M. F. Como determinar a distribuição potencial de espécies sob uma abordagem conservacionista? Megadiversidade (2009).48.Hijmans, R. J. et al. DIVA-GIS, version 5.2. A geographic information system for the analysis of biodiversity data. Manual. . vol. 1 (International Potato Center, 2005).49.Borcard, D., Gillet, F. & Legendre, P. Numerical Ecology with R. Numerical Ecology with R (Springer New York, 2011). doi:https://doi.org/10.1007/978-1-4419-7976-6.50.Serra, B. D. V., De Marco, P. J., Nóbrega, C. C. & Campos, L. A. D. O. Modeling potential geographical distribution of the wild nests of Melipona capixaba Moure & Camargo, 1994 ( Hymenoptera, Apidae ): Conserving Isolated Populations in Mountain Habitats. Nat. e Conserv. 10, 199–206 (2012).Article 

    Google Scholar 
    51.Phillips, S. J., Anderson, R. P. & Schapire, R. E. Maximum entropy modeling of species geographic distributions. Ecol. Model. https://doi.org/10.1016/j.ecolmodel.2005.03.026 (2006).Article 

    Google Scholar 
    52.Swets, J. Measuring the accuracy of diagnostic systems. Science (80-. ). 240, 1285–1293 (1988).53.Girardello, M., Griggio, M., Whittingham, M. J. & Rushton, S. P. Identifying important areas for butterfly conservation in Italy. Anim. Conserv. https://doi.org/10.1111/j.1469-1795.2008.00216.x (2009).Article 

    Google Scholar 
    54.Allouche, O., Tsoar, A. & Kadmon, R. Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS). J. Appl. Ecol. 43, 1223–1232 (2006).Article 

    Google Scholar 
    55.Vieira, T. B., Mendes, P. & Oprea, M. Priority areas for bat conservation in the state of Espírito Santo, southeastern Brazil. Neotrop. Biol. Conserv. 7, 88–96 (2012).Article 

    Google Scholar 
    56.Delgado-Jaramillo, M., Aguiar, L. M. S., Machado, R. B. & Bernard, E. Assessing the distribution of a species-rich group in a continental-sized megadiverse country: Bats in Brazil. Divers. Distrib. 26, 632–643 (2020).Article 

    Google Scholar 
    57.Destro, G. F. G., de Fernandes, V., de Andrade, A. F. A., De Marco, P. & Terribile, L. C. Back home? Uncertainties for returning seized animals to the source-areas under climate change. Glob. Chang. Biol. 25, 3242–3253 (2019).ADS 
    Article 

    Google Scholar 
    58.Elith, J. et al. Novel methods improve prediction of species’ distributions from occurrence data. Ecography (Cop.). (2006) doi:https://doi.org/10.1111/j.2006.0906-7590.04596.x.59.Pearson, R. G., Raxworthy, C. J., Nakamura, M. & Townsend Peterson, A. Predicting species distributions from small numbers of occurrence records: a test case using cryptic geckos in Madagascar. J. Biogeogr. 34, 102–117 (2006).60.de Andrade, A. F. A., Velazco, S. J. E. & De Marco, P. Niche mismatches can impair our ability to predict potential invasions. Biol. Invasions 21, 3135–3150 (2019).Article 

    Google Scholar 
    61.Velazco, S. J. E., Villalobos, F., Galvão, F. & De Marco Júnior, P. A dark scenario for Cerrado plant species: Effects of future climate, land use and protected areas ineffectiveness. Divers. Distrib. 25, 660–673 (2019).62.Baselga, A. Partitioning the turnover and nestedness components of beta diversity. Glob. Ecol. Biogeogr. https://doi.org/10.1111/j.1466-8238.2009.00490.x (2010).Article 

    Google Scholar 
    63.Moilanen, A. et al. Prioritizing multiple-use landscapes for conservation : methods for large multi-species planning problems. Proc. R. Soc. 272, 1885–1891 (2005).
    Google Scholar 
    64.Moilanen, A. et al. Zonation spatial conservation planning framework and software v. 3.1, User manual. (2012).65.Moilanen, A. Landscape zonation, benefit functions and target-based planning: unifying reserve selection strategies. Biol. Conserv. 134, 571–579 (2007).Article 

    Google Scholar 
    66.Carvalho, A. R. de. Método de Monte Carlo e Aplicações. Repositório Inst. da Univ. Fed. Flum. 84 (2017).67.Feinleib, M. & Zar, J. H. Biostatistical analysis. J. Am. Stat. Assoc. https://doi.org/10.2307/2285423 (1975).Article 

    Google Scholar  More

  • in

    Rebound in China’s coastal wetlands following conservation and restoration

    1.Ma, Z. J. et al. Rethinking China’s new great wall. Science 346, 912–914 (2014).CAS 
    Article 

    Google Scholar 
    2.Murray, N. J. et al. The global distribution and trajectory of tidal flats. Nature 565, 222–225 (2019).Article 
    CAS 

    Google Scholar 
    3.Wang, X. et al. Tracking annual changes of coastal tidal flats in China during 1986–2016 through analyses of Landsat images with Google Earth Engine. Remote Sens. Environ. 238, 110987 (2020).Article 

    Google Scholar 
    4.Blum, M. D. & Roberts, H. H. Drowning of the Mississippi Delta due to insufficient sediment supply and global sea-level rise. Nat. Geosci. 2, 488–491 (2009).CAS 
    Article 

    Google Scholar 
    5.Murray, N. J., Clemens, R. S., Phinn, S. R., Possingham, H. P. & Fuller, R. A. Tracking the rapid loss of tidal wetlands in the Yellow Sea. Front. Ecol. Environ. 12, 267–272 (2014).Article 

    Google Scholar 
    6.Gedan, K. B., Silliman, B. R. & Bertness, M. D. Centuries of human-driven change in salt marsh ecosystems. Ann. Rev. Mar. Sci. 1, 117–141 (2009).Article 

    Google Scholar 
    7.Syvitski, J. P. M., Vörösmarty, C. J., Kettner, A. J. & Green, P. Impact of humans on the flux of terrestrial sediment to the global coastal ocean. Science 308, 376–380 (2005).CAS 
    Article 

    Google Scholar 
    8.Cui, B., He, Q., Gu, B., Bai, J. & Liu, X. China’s coastal wetlands: understanding environmental changes and human impacts for management and conservation. Wetlands 36, 1–9 (2016).Article 

    Google Scholar 
    9.Gong, P. et al. Stable classification with limited sample: transferring a 30-m resolution sample set collected in 2015 to mapping 10-m resolution global land cover in 2017. Sci. Bull. 64, 370–373 (2019).Article 

    Google Scholar 
    10.Han, Q., Niu, Z., Wu, M. & Wang, J. Remote-sensing monitoring and analysis of China intertidal zone changes based on tidal correction. Sci. Bull. 64, 456–473 (2019).
    Google Scholar 
    11.Mao, D. et al. National wetland mapping in China: a new product resulting from object-based and hierarchical classification of Landsat 8 OLI images. ISPRS J. Photogramm. Remote Sens. 164, 11–25 (2020).Article 

    Google Scholar 
    12.Wang, X. et al. Mapping coastal wetlands of China using time series Landsat images in 2018 and Google Earth Engine. ISPRS J. Photogramm. Remote Sens. 163, 312–326 (2020).Article 

    Google Scholar 
    13.Mcowen, C. J. et al. A global map of saltmarshes. Biodivers. Data J. 5, e11764 (2017).Article 

    Google Scholar 
    14.Giri, C. et al. Status and distribution of mangrove forests of the world using Earth observation satellite data. Glob. Ecol. Biogeogr. 20, 154–159 (2011).Article 

    Google Scholar 
    15.Chen, Y. et al. Effects of reclamation and natural changes on coastal wetlands bordering China’s Yellow Sea from 1984 to 2015. Land Degrad. Dev. 30, 1533–1544 (2019).Article 

    Google Scholar 
    16.Hu, Y. et al. Mapping coastal salt marshes in China using time series of Sentinel-1 SAR. ISPRS J. Photogramm. Remote Sens. 173, 122–134 (2021).Article 

    Google Scholar 
    17.Zhang, X. et al. Quantifying expansion and removal of Spartina alterniflora on Chongming Island, China, using time series Landsat images during 1995–2018. Remote Sens. Environ. 247, 111916 (2020).18.Chen, B. Q. et al. A mangrove forest map of China in 2015: analysis of time series Landsat 7/8 and Sentinel-1A imagery in Google Earth Engine cloud computing platform. ISPRS J. Photogramm. Remote Sens. 131, 104–120 (2017).Article 

    Google Scholar 
    19.Hu, L., Li, W. & Xu, B. Monitoring mangrove forest change in China from 1990 to 2015 using Landsat-derived spectral-temporal variability metrics. Int. J. Appl. Earth Obs. Geoinf. 73, 88–98 (2018).Article 

    Google Scholar 
    20.Jia, M., Wang, Z., Zhang, Y., Mao, D. & Wang, C. Monitoring loss and recovery of mangrove forests during 42 years: the achievements of mangrove conservation in China. Int. J. Appl. Earth Obs. Geoinf. 73, 535–545 (2018).Article 

    Google Scholar 
    21.Jia, M. et al. Rapid, robust, and automated mapping of tidal flats in China using time series Sentinel-2 images and Google Earth Engine. Remote Sens. Environ. 255, 112285 (2021).Article 

    Google Scholar 
    22.Ma, T., Li, X., Bai, J. & Cui, B. Tracking three decades of land use and land cover transformation trajectories in China’s large river deltas. Land Degrad. Dev. 30, 799–810 (2019).Article 

    Google Scholar 
    23.Wang, K. Evolution of Yellow River delta coastline based on remote sensing from 1976 to 2014, China. Chin. Geogr. Sci. 29, 181–191 (2019).Article 

    Google Scholar 
    24.Zhao, Y. F. et al. Assessing natural and anthropogenic influences on water discharge and sediment load in the Yangtze River, China. Sci. Total Environ. 607, 920–932 (2017).Article 
    CAS 

    Google Scholar 
    25.Yim, J. et al. Analysis of forty years long changes in coastal land use and land cover of the Yellow Sea: the gains or losses in ecosystem services. Environ. Pollut. 241, 74–84 (2018).CAS 
    Article 

    Google Scholar 
    26.Wang, S. et al. Reduced sediment transport in the Yellow River due to anthropogenic changes. Nat. Geosci. 9, 38–41 (2016).
    Google Scholar 
    27.Chen, Y. et al. Land claim and loss of tidal flats in the Yangtze Estuary. Sci. Rep. 6, 24018 (2016).CAS 
    Article 

    Google Scholar 
    28.Yang, M. et al. Spatio-temporal characterization of a reclamation settlement in the Shanghai coastal area with time series analyses of X-, C-, and L-band SAR datasets. Remote Sens. 10, 329 (2018).29.Han, X., Pan, J. & Devlin, A. T. Remote sensing study of wetlands in the Pearl River Delta during 1995–2015 with the support vector machine method. Front. Earth Sci. 12, 521–531 (2018).Article 

    Google Scholar 
    30.Liu, L., Xu, W., Yue, Q., Teng, X. & Hu, H. Problems and countermeasures of coastline protection and utilization in China. Ocean Coast. Manag. 153, 124–130 (2018).Article 

    Google Scholar 
    31.Yunxuan, Z. et al. Degradation of coastal wetland ecosystem in China: drivers, impacts, and strategies. Bull. Chin. Acad. Sci. 31, 1157–1166 (2016).
    Google Scholar 
    32.Jiang, T. T., Pan, J. F., Pu, X. M., Wang, B. & Pan, J. J. Current status of coastal wetlands in China: degradation, restoration, and future management. Estuar. Coast. Shelf Sci. 164, 265–275 (2015).Article 

    Google Scholar 
    33.Sun, Z. et al. China’s coastal wetlands: conservation history, implementation efforts, existing issues and strategies for future improvement. Environ. Int. 79, 25–41 (2015).Article 

    Google Scholar 
    34.Ren, C. et al. Rapid expansion of coastal aquaculture ponds in China from Landsat observations during 1984–2016. Int. J. Appl. Earth Obs. Geoinf. 82, 101902 (2019).35.Gu, J. et al. Losses of salt marsh in China: trends, threats and management. Estuar. Coast. Shelf Sci. 214, 98–109 (2018).Article 

    Google Scholar 
    36.Wang, W., Liu, H., Li, Y. & Su, J. Development and management of land reclamation in China. Ocean Coast. Manag. 102, 415–425 (2014).Article 

    Google Scholar 
    37.Barbier, E. B. et al. The value of estuarine and coastal ecosystem services. Ecol. Monogr. 81, 169–193 (2011).Article 

    Google Scholar 
    38.Barbier, E. B. A global strategy for protecting vulnerable coastal populations. Science 345, 1250–1251 (2014).CAS 
    Article 

    Google Scholar 
    39.He, Q. et al. Economic development and coastal ecosystem change in China. Sci. Rep. 4, 5995 (2014).40.Zhou, C. et al. Preliminary analysis of C sequestration potential of blue carbon ecosystems on Chinese coastal zone. Sci. China Life Sci. 46, 475–486 (2016).
    Google Scholar 
    41.Zhang, Q. et al. Propagule types and environmental stresses matter in saltmarsh plant restoration. Ecol. Eng. 143, 105693 (2020).Article 

    Google Scholar 
    42.Cui, B., Yang, Q., Yang, Z. & Zhang, K. Evaluating the ecological performance of wetland restoration in the Yellow River Delta, China. Ecol. Eng. 35, 1090–1103 (2009).Article 

    Google Scholar 
    43.Pan, X. Research on Xi Jinping’s thought of ecological civilization and environment sustainable development. IOP Conf. Ser. Earth Environ. Sci. 153, 062067 (2018).44.Hansen, M. H., Li, H. & Svarverud, R. Ecological civilization: interpreting the Chinese past, projecting the global future. Glob. Environ. Change. 53, 195–203 (2018).Article 

    Google Scholar 
    45.Moreno-Mateos, D., Power, M. E., Comín, F. A. & Yockteng, R. Structural and functional loss in restored wetland ecosystems. PLoS Biol. 10, e1001247 (2012).CAS 
    Article 

    Google Scholar 
    46.He, Q. Conservation: ‘No net loss’ of wetland quantity and quality. Curr. Biol. 29, R1070–R1072 (2019).CAS 
    Article 

    Google Scholar 
    47.Gong, P., Li, X. & Zhang, W. 40-year (1978-2017) human settlement changes in China reflected by impervious surfaces from satellite remote sensing. Sci. Bull. 64, 756–763 (2019).Article 

    Google Scholar 
    48.Wang, X. et al. Gainers and losers of surface and terrestrial water resources in China during 1989–2016. Nat. Commun. 11, 3471 (2020).49.Zou, Z. H. et al. Divergent trends of open-surface water body area in the contiguous United States from 1984 to 2016. Proc. Natl Acad. Sci. USA 115, 3810–3815 (2018).CAS 
    Article 

    Google Scholar  More

  • in

    Contribution of historical herbarium small RNAs to the reconstruction of a cassava mosaic geminivirus evolutionary history

    1.Stukenbrock, E. H. & McDonald, B. A. The origins of plant pathogens in agro-ecosystems. Annu. Rev. Phytopathol. https://doi.org/10.1146/annurev.phyto.010708.154114 (2008).Article 
    PubMed 

    Google Scholar 
    2.Savary, S., Ficke, A., Aubertot, J. N. & Hollier, C. Crop losses due to diseases and their implications for global food production losses and food security. Food Secur. https://doi.org/10.1007/s12571-012-0200-5 (2012).Article 

    Google Scholar 
    3.Strange, R. N. & Scott, P. R. Plant disease: a threat to global food security. Annu. Rev. Phytopathol. https://doi.org/10.1146/annurev.phyto.43.113004.133839 (2005).Article 
    PubMed 

    Google Scholar 
    4.Anderson, P. K. et al. Emerging infectious diseases of plants: pathogen pollution, climate change and agrotechnology drivers. Trends Ecol. Evol. https://doi.org/10.1016/j.tree.2004.07.021 (2004).Article 
    PubMed 

    Google Scholar 
    5.Scholthof, K. B. G. et al. Top 10 plant viruses in molecular plant pathology. Mol. Plant Pathol. https://doi.org/10.1111/j.1364-3703.2011.00752.x (2011).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    6.Stukenbrock, E. H. & Bataillon, T. A population genomics perspective on the emergence and adaptation of new plant pathogens in agro-ecosystems. PLoS Pathog. https://doi.org/10.1371/journal.ppat.1002893 (2012).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    7.Gilligan, C. A. Sustainable agriculture and plant diseases: an epidemiological perspective. Philos. Trans. R. Soc. B: Biol. Sci. https://doi.org/10.1098/rstb.2007.2181 (2008).Article 

    Google Scholar 
    8.Li, L. M., Grassly, N. C. & Fraser, C. Genomic analysis of emerging pathogens: methods, application and future trends. Genome Biol.ogy https://doi.org/10.1186/s13059-014-0541-9 (2014).Article 

    Google Scholar 
    9.Lemey, P., Rambaut, A., Drummond, A. J. & Suchard, M. A. Bayesian phylogeography finds its roots. PLoS Comput. Biol. https://doi.org/10.1371/journal.pcbi.1000520 (2009).MathSciNet 
    Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    10.Lefeuvre, P. et al. The spread of tomato yellow leaf curl virus from the middle east to the world. PLoS Pathog. https://doi.org/10.1371/journal.ppat.1001164 (2010).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    11.Monjane, A. L. et al. Reconstructing the history of maize streak virus strain A dispersal tor reveal diversification hot spots and its origin in southern Africa. J. Virol. https://doi.org/10.1128/jvi.00640-11 (2011).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    12.Trovao, N. S. et al. Host ecology determines the dispersal patterns of a plant virus. Virus Evol. https://doi.org/10.1093/ve/vev016 (2015).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    13.Rakotomalala, M. et al. Comparing patterns and scales of plant virus phylogeography: rice yellow mottle virus in Madagascar and in continental Africa. Virus Evol. https://doi.org/10.1093/ve/vez023 (2019).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    14.Gibbs, A. J., Fargette, D., García-Arenal, F. & Gibbs, M. J. Time – The emerging dimension of plant virus studies. J General Virol. https://doi.org/10.1099/vir.0.015925-0 (2010).Article 

    Google Scholar 
    15.Simmonds, P., Aiewsakun, P. & Katzourakis, A. Prisoners of war: host adaptation and its constraints on virus evolution. Nat. Rev. Microbiol. https://doi.org/10.1038/s41579-018-0120-2 (2019).Article 
    PubMed 

    Google Scholar 
    16.Jones, R. A. C., Boonham, N., Adams, I. P. & Fox, A. Historical virus isolate collections: an invaluable resource connecting plant virology’s pre-sequencing and post-sequencing eras. Plant Pathol. 70, 235–248 (2021).Article 

    Google Scholar 
    17.Smith, O. et al. A complete ancient RNA genome: Identification, reconstruction and evolutionary history of archaeological Barley Stripe Mosaic Virus. Sci. Rep. https://doi.org/10.1038/srep04003 (2014).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    18.Malmstrom, C. M., Shu, R., Linton, E. W., Newton, L. A. & Cook, M. A. Barley yellow dwarf viruses (BYDVs) preserved in herbarium specimens illuminate historical disease ecology of invasive and native grasses. J. Ecol. https://doi.org/10.1111/j.1365-2745.2007.01307.x (2007).Article 

    Google Scholar 
    19.Peyambari, M., Warner, S., Stoler, N., Rainer, D. & Roossinck, M. J. A 1000-Year-old RNA virus. J. Virol. 93, e01188-18 (2019).CAS 
    Article 

    Google Scholar 
    20.Adams, I. P. et al. Next-generation sequencing and metagenomic analysis: a universal diagnostic tool in plant virology. Mol. Plant Pathol. https://doi.org/10.1111/j.1364-3703.2009.00545.x (2009).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    21.Vayssier-Taussat, M. et al. Shifting the paradigm from pathogens to pathobiome new concepts in the light of meta-omics. Front. Cell. Infect. Microbiol. https://doi.org/10.3389/fcimb.2014.00029 (2014).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    22.Massart, S., Olmos, A., Jijakli, H. & Candresse, T. Current impact and future directions of high throughput sequencing in plant virus diagnostics. Virus Res. https://doi.org/10.1016/j.virusres.2014.03.029 (2014).Article 
    PubMed 

    Google Scholar 
    23.Roossinck, M. J., Martin, D. P. & Roumagnac, P. Plant virus metagenomics: advances in virus discovery. Phytopathology https://doi.org/10.1094/PHYTO-12-14-0356-RVW (2015).Article 
    PubMed 

    Google Scholar 
    24.Kreuze, J. F. et al. Complete viral genome sequence and discovery of novel viruses by deep sequencing of small RNAs: a generic method for diagnosis, discovery and sequencing of viruses. Virology https://doi.org/10.1016/j.virol.2009.03.024 (2009).Article 
    PubMed 

    Google Scholar 
    25.Pooggin, M. M. Small RNA-omics for plant virus identification, virome reconstruction, and antiviral defense characterization. Front. Microbiol. https://doi.org/10.3389/fmicb.2018.02779 (2018).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    26.Hartung, J. S. et al. History and diversity of Citrus Leprosis virus recorded in herbarium specimens. Phytopathology https://doi.org/10.1094/PHYTO-03-15-0064-R (2015).Article 
    PubMed 

    Google Scholar 
    27.Golyaev, V., Candresse, T., Rabenstein, F. & Pooggin, M. M. Plant virome reconstruction and antiviral RNAi characterization by deep sequencing of small RNAs from dried leaves. Sci. Rep. https://doi.org/10.1038/s41598-019-55547-3 (2019).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    28.Patil, B. L. & Fauquet, C. M. Cassava mosaic geminiviruses: actual knowledge and perspectives. Mol. Plant Pathol. https://doi.org/10.1111/j.1364-3703.2009.00559.x (2009).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    29.Legg, J. P., Owor, B., Sseruwagi, P. & Ndunguru, J. Cassava mosaic virus disease in east and central Africa: epidemiology and management of a regional pandemic. Adv. Virus Res. https://doi.org/10.1016/S0065-3527(06)67010-3 (2006).Article 
    PubMed 

    Google Scholar 
    30.Wang, H. L. et al. First report of Sri Lankan cassava mosaic virus infecting cassava in Cambodia. Plant Dis. https://doi.org/10.1094/PDIS-10-15-1228-PDN (2016).Article 
    PubMed 

    Google Scholar 
    31.Minato, N. et al. Surveillance for sri lankan cassava mosaic virus (SLCMV) in Cambodia and Vietnam one year after its initial detection in a single plantation in 2015. PLoS One https://doi.org/10.1371/journal.pone.0212780 (2019).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    32.Mugerwa, H., Wang, H. L., Sseruwagi, P., Seal, S. & Colvin, J. Whole-genome single nucleotide polymorphism and mating compatibility studies reveal the presence of distinct species in sub-Saharan Africa Bemisia tabaci whiteflies. Insect Sci. https://doi.org/10.1111/1744-7917.12881 (2020).Article 
    PubMed 

    Google Scholar 
    33.Ntawuruhunga, P. et al. Incidence and severity of cassava mosaic disease in the Republic of Congo. African Crop Sci. J. https://doi.org/10.4314/acsj.v15i1.54405 (2010).Article 

    Google Scholar 
    34.Zinga, I. et al. Epidemiological assessment of cassava mosaic disease in Central African Republic reveals the importance of mixed viral infection and poor health of plant cuttings. Crop Prot. https://doi.org/10.1016/j.cropro.2012.10.010 (2013).Article 

    Google Scholar 
    35.Jeske, H. Geminiviruses. Curr. Topics Microbiol. Immunol. https://doi.org/10.1007/978-3-540-70972-5_11 (2009).Article 

    Google Scholar 
    36.Vanitharani, R., Chellappan, P. & Fauquet, C. M. Geminiviruses and RNA silencing. Trends Plant Sci. https://doi.org/10.1016/j.tplants.2005.01.005 (2005).Article 
    PubMed 

    Google Scholar 
    37.Aregger, M. et al. Primary and secondary siRNAs in geminivirus-induced gene silencing. PLoS Pathog. https://doi.org/10.1371/journal.ppat.1002941 (2012).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    38.Olsen, K. M. & Schaal, B. A. Evidence on the origin of cassava: Phylogeography of Manihot esculenta. Proc. Natl. Acad. Sci. USA https://doi.org/10.1073/pnas.96.10.5586 (1999).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    39.Fauquet, C. African cassava mosaic virus: etiology, epidemiology, and control. Plant Dis. https://doi.org/10.1094/pd-74-0404 (1990).Article 

    Google Scholar 
    40.Legg, J. P. & Fauquet, C. M. Cassava mosaic geminiviruses in Africa. Plant Mol. Biol. https://doi.org/10.1007/s11103-004-1651-7 (2004).Article 
    PubMed 

    Google Scholar 
    41.De Bruyn, A. et al. Divergent evolutionary and epidemiological dynamics of cassava mosaic geminiviruses in Madagascar. BMC Evol. Biol. https://doi.org/10.1186/s12862-016-0749-2 (2016).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    42.Weiß, C. L. et al. Temporal patterns of damage and decay kinetics of dna retrieved from plant herbarium specimens. R. Soc. Open Sci. https://doi.org/10.1098/rsos.160239 (2016).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    43.Chellappan, P., Vanitharani, R., Ogbe, F. & Fauquet, C. M. Effect of temperature on geminivirus-induced RNA silencing in plants. Plant Physiol. https://doi.org/10.1104/pp.105.066563 (2005).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    44.Smith, O. & Gilbert, M. T. P. Ancient RNA. in (2018). doi:https://doi.org/10.1007/13836_2018_17.45.Filloux, D. et al. The genomes of many yam species contain transcriptionally active endogenous geminiviral sequences that may be functionally expressed. Virus Evol. https://doi.org/10.1093/ve/vev002 (2015).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    46.Sharma, V. et al. Large-scale survey reveals pervasiveness and potential function of endogenous geminiviral sequences in plants. Virus Evol. https://doi.org/10.1093/ve/veaa071 (2020).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    47.Bredeson, J. V. et al. Sequencing wild and cultivated cassava and related species reveals extensive interspecific hybridization and genetic diversity. Nat. Biotechnol. https://doi.org/10.1038/nbt.3535 (2016).Article 
    PubMed 

    Google Scholar 
    48.Serfraz, S. et al. Insertion of Badnaviral DNA in the Late Blight Resistance Gene (R1a) of Brinjal Eggplant (Solanum melongena). Front. Plant Sci. https://doi.org/10.3389/fpls.2021.683681 (2021).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    49.Lefeuvre, P. et al. Evolutionary time-scale of the begomoviruses: evidence from integrated sequences in the Nicotiana genome. PLoS One https://doi.org/10.1371/journal.pone.0019193 (2011).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    50.Martin, D. P., Murrell, B., Golden, M., Khoosal, A. & Muhire, B. RDP4: detection and analysis of recombination patterns in virus genomes. Virus Evol. https://doi.org/10.1093/ve/vev003 (2015).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    51.Murray, G. G. R. et al. The effect of genetic structure on molecular dating and tests for temporal signal. Methods Ecol. Evol. 7, 80–89 (2016).Article 

    Google Scholar 
    52.Drummond, A. J. & Rambaut, A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol. https://doi.org/10.1186/1471-2148-7-214 (2007).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    53.Yoshida, K. et al. Mining herbaria for plant pathogen genomes: back to the future. PLoS Pathog. https://doi.org/10.1371/journal.ppat.1004028 (2014).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    54.Dufrénoy, J. & Hédin, L. . La. Mosaïque des feuilles du Manioc au Cameroun. J. d’agriculture Tradit. Bot. appliquée 94, 361–365 (1929).
    Google Scholar 
    55.Duffy, S. & Holmes, E. C. Validation of high rates of nucleotide substitution in geminiviruses: phylogenetic evidence from East African cassava mosaic viruses. J. Gen. Virol. 90, 1539–1547 (2009).CAS 
    Article 

    Google Scholar 
    56.Worobey, M. et al. Direct evidence of extensive diversity of HIV-1 in Kinshasa by 1960. Nature https://doi.org/10.1038/nature07390 (2008).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    57.Mühlemann, B. et al. Ancient hepatitis B viruses from the Bronze Age to the Medieval period. Nature https://doi.org/10.1038/s41586-018-0097-z (2018).Article 
    PubMed 

    Google Scholar 
    58.Toppinen, M. et al. Bones hold the key to DNA virus history and epidemiology. Sci. Rep. https://doi.org/10.1038/srep17226 (2015).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    59.Gilbert, M. T. P., Bandelt, H. J., Hofreiter, M. & Barnes, I. Assessing ancient DNA studies. Trends Ecol. Evol. https://doi.org/10.1016/j.tree.2005.07.005 (2005).Article 
    PubMed 

    Google Scholar 
    60.Inoue-Nagata, A. K., Albuquerque, L. C., Rocha, W. B. & Nagata, T. A simple method for cloning the complete begomovirus genome using the bacteriophage φ29 DNA polymerase. J. Virol. Methods https://doi.org/10.1016/j.jviromet.2003.11.015 (2004).Article 
    PubMed 

    Google Scholar 
    61.Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics https://doi.org/10.1093/bioinformatics/btu170 (2014).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    62.Zheng, Y. et al. VirusDetect: An automated pipeline for efficient virus discovery using deep sequencing of small RNAs. Virology https://doi.org/10.1016/j.virol.2016.10.017 (2017).Article 
    PubMed 

    Google Scholar 
    63.Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics https://doi.org/10.1093/bioinformatics/btp324 (2009).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    64.Langmead, B., Trapnell, C., Pop, M. & Salzberg, S. L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. https://doi.org/10.1186/gb-2009-10-3-r25 (2009).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    65.Jónsson, H., Ginolhac, A., Schubert, M., Johnson, P. L. F. & Orlando, L. MapDamage2.0: Fast approximate Bayesian estimates of ancient DNA damage parameters. in Bioinformatics (2013). doi:https://doi.org/10.1093/bioinformatics/btt193.66.Broad Institute. Picard Tools – By Broad Institute. Github (2009).67.Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics https://doi.org/10.1093/bioinformatics/btq033 (2010).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    68.Krzywinski, M. et al. Circos: an information aesthetic for comparative genomics. Genome Res. https://doi.org/10.1101/gr.092759.109 (2009).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    69.Depristo, M. A. et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat. Genet. https://doi.org/10.1038/ng.806 (2011).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    70.Bankevich, A. et al. SPAdes: A new genome assembly algorithm and its applications to single-cell sequencing. J. Comput. Biol. https://doi.org/10.1089/cmb.2012.0021 (2012).MathSciNet 
    Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    71.Katoh, K. & Standley, D. M. MAFFT multiple sequence alignment software version 7: Improvements in performance and usability. Mol. Biol. Evol. https://doi.org/10.1093/molbev/mst010 (2013).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    72.Stamatakis, A. RAxML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313 (2014).CAS 
    Article 

    Google Scholar 
    73.Darriba, D., Taboada, G. L., Doallo, R. & Posada, D. JModelTest 2: More models, new heuristics and parallel computing. Nat. Methods https://doi.org/10.1038/nmeth.2109 (2012).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    74.Jombart, T. & Dray, S. Adephylo: Exploratory analyses for the phylogenetic comparative method. Bioinformatics (2010).75.Duchêne, S., Duchêne, D., Holmes, E. C. & Ho, S. Y. W. The performance of the date-randomization test in phylogenetic analyses of time-structured virus data. Mol. Biol. Evol. 32, 1895–1906 (2015).Article 

    Google Scholar 
    76.Rieux, A. & Khatchikian, C. E. Tipdatingbeast: an r package to assist the implementation of phylogenetic tip-dating tests using beast. Mol. Ecol. Resour. https://doi.org/10.1111/1755-0998.12603 (2017).Article 
    PubMed 

    Google Scholar 
    77.Raftery, A. E. Approximate Bayes factors and accounting for model uncertainty in generalised linear models. Biometrika https://doi.org/10.1093/biomet/83.2.251 (1996).MathSciNet 
    Article 
    MATH 

    Google Scholar 
    78.Ho, S. Y. W. & Shapiro, B. Skyline-plot methods for estimating demographic history from nucleotide sequences. Mol. Ecol. Resour. https://doi.org/10.1111/j.1755-0998.2011.02988.x (2011).Article 
    PubMed 

    Google Scholar 
    79.Rambaut, A., Drummond, A. J., Xie, D., Baele, G. & Suchard, M. A. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst. Biol. (2018) doi:https://doi.org/10.1093/sysbio/syy032. More