Ecology

1.
Chavez, F. P. & Messié, M. A comparison of eastern boundary upwelling ecosystems. Prog. Oceanogr. 83, 80–96 (2009).
ADS  Article  Google Scholar 
2.
Auger, P.-A., Gorgues, T., Machu, E., Aumont, O. & Brehmer, P. What drives the spatial variability of primary productivity and matter fluxes in the north-west African upwelling system? A modelling approach. Biogeosciences 13, 6419–6440 (2016).
ADS  Article  Google Scholar 

3.
Benazzouz, A. et al. An improved coastal upwelling index from sea surface temperature using satellite-based approach—The case of the Canary Current upwelling system. Cont. Shelf Res. 81, 38–54 (2014).
ADS  Article  Google Scholar 

4.
Citeau, J., Finaud, L., Cammas, J. & Demarcq, H. Questions relative to ITCZ migrations over the tropical Atlantic ocean, sea surface temperature and Senegal River runoff. Meteorol. Atmos. Phys. 41, 181–190 (1989).
ADS  Article  Google Scholar 

5.
Maloney, E. & Shaman, J. Intraseasonal variability of the West African Monsoon and Atlantic ITCZ. J. Clim. 21, 2898–2918 (2008).
ADS  Article  Google Scholar 

6.
Herbland, A. & Voituriez, B. L. production primaire dans l’upwelling mauritanien en mars 1973. Cahiers ORSTOM 12, 187–201 (1974).
CAS  Google Scholar 

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

8.
Hjermann, D. Ø., Ottersen, G. & Stenseth, N. C. Competition among fishermen and fish causes the collapse of Barents Sea capelin. PNAS 101, 11679–11684 (2004).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

9.
Fréon, P., Cury, P., Shannon, L. & Roy, C. Sustainable exploitation of small pelagic fish stocks challenged by environmental and ecosystem changes: A review. Bull. Mar. Sci. 76, 385–462 (2005).
Google Scholar 

10.
Schwartzlose, R. A. et al. Worldwide large-scale fluctuations of sardine and anchovy populations. Afr. J. Mar. Sci. 21, 289–347 (1999).
Article  Google Scholar 

11.
Hofstede, R. T., Dickey-Collas, M., Mantingh, I. T. & Wague, A. The link between migration, the reproductive cycle and condition of Sardinella aurita off Mauritania, north-west Africa. J. Fish Biol. 71, 1293–1302 (2007).
Article  Google Scholar 

12.
Zeeberg, J., Corten, A., Tjoe-Awie, P., Coca, J. & Hamady, B. Climate modulates the effects of Sardinella aurita fisheries off Northwest Africa. Fish. Res. 1, 65–75 (2008).
Article  Google Scholar 

13.
Gibson, R. N., Atkinson, R. J. A. & Gordon, J. D. M. Oceanography and Marine Biology, Vol. 47 (2009).

14.
Hays, G. C., Richardson, A. J. & Robinson, C. Climate change and marine plankton. Trends Ecol. Evol. 20, 337–344 (2005).
PubMed  Article  PubMed Central  Google Scholar 

15.
Ndoye, S. et al. Dynamics of a “low-enrichment high-retention” upwelling center over the southern Senegal shelf. Geophys. Res. Lett. 44, 5034–5043 (2017).
ADS  Article  Google Scholar 

16.
Behagle, N. et al. Acoustic distribution of discriminated micronektonic organisms from a bi-frequency processing: The case study of eastern Kerguelen oceanic waters. Prog. Oceanogr. 156, 276–289 (2017).
Article  Google Scholar 

17.
Benoit-Bird, K. & Au, W. Diel migration dynamics of an island-associated sound-scattering layer. Deep Sea Res. Part I 51, 707–719 (2004).
Article  Google Scholar 

18.
Sato, M. & Benoit-Bird, K. J. Spatial variability of deep scattering layers shapes the Bahamian mesopelagic ecosystem. Mar. Ecol. Prog. Ser. 580, 69–82 (2017).
ADS  Article  Google Scholar 

19.
Algueró-Muñiz, M. et al. Ocean acidification effects on mesozooplankton community development: Results from a long-term mesocosm experiment. PLoS ONE 12, e0175851 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

20.
Matlab R 2018a. The Math Works, Inc. (MATLAB & Simulink – MathWorks., 2018).

21.
Hegerl, G. C. et al. Causes of climate change over the historical record. Environ. Res. Lett. 14, 123006 (2019).
ADS  CAS  Article  Google Scholar 

22.
Belkin, I. M. Rapid warming of large marine ecosystems. Prog. Oceanogr. 81, 207–213 (2009).
ADS  Article  Google Scholar 

23.
Valdés, L. & Déniz-González, I. Oceanographic and biological features in the Canary Current Large Marine Ecosystem, Vol. 115 (2015).

24.
Gómez-Letona, M., Ramos, A. G., Coca, J. & Arístegui, J. Trends in primary production in the canary current upwelling system—A regional perspective comparing remote sensing models. Front. Mar. Sci. 4, 370 (2017).
Article  Google Scholar 

25.
Bakun, A. Global climate change and intensification of coastal ocean upwelling. Science 247, 198–201 (1990).
ADS  CAS  PubMed  Article  Google Scholar 

26.
Benazzouz, A., Demarcq, H. & González-Nuevo, G. Oceanographic and biological features in the Canary current large marine ecosystem. (2015).

27.
Behrenfeld, M. J. et al. Climate-driven trends in contemporary ocean productivity. Nature 444, 752–755 (2006).
ADS  CAS  PubMed  Article  Google Scholar 

28.
Boyce, D. G., Lewis, M. R. & Worm, B. Global phytoplankton decline over the past century. Nature 466, 591–596 (2010).
ADS  CAS  PubMed  Article  Google Scholar 

29.
Hofmann, M., Worm, B., Rahmstorf, S. & Schellnhuber, H. J. Declining ocean chlorophyll under unabated anthropogenic CO2 emissions. Environ. Res. Lett. 6, 034035 (2011).
ADS  Article  CAS  Google Scholar 

30.
Lewandowska, A. et al. Effects of sea surface warming on marine plankton. Ecol. Lett. 17, 614–623 (2014).
PubMed  Article  PubMed Central  Google Scholar 

31.
Arístegui, J. et al. Sub-regional ecosystem variability in the Canary Current upwelling. Prog. Oceanogr. 83, 33–48 (2009).
ADS  Article  Google Scholar 

32.
Demarcq, H. Trends in primary production, sea surface temperature and wind in upwelling systems (1998–2007). Prog. Oceanogr. 83, 376–385 (2009).
ADS  Article  Google Scholar 

33.
Barton, A. D., Irwin, A. J., Finkel, Z. V. & Stock, C. A. Anthropogenic climate change drives shift and shuffle in North Atlantic phytoplankton communities. PNAS 113, 2964–2969 (2016).
ADS  CAS  PubMed  Article  Google Scholar 

34.
Jacob, B. G. et al. Major changes in diatom abundance, productivity, and net community metabolism in a windier and dryer coastal climate in the southern Humboldt Current. Prog. Oceanogr. 168, 196–209 (2018).
ADS  Article  Google Scholar 

35.
Jacox, M. G., Hazen, E. L. & Bograd, S. J. optimal environmental conditions and anomalous ecosystem responses: Constraining bottom-up controls of phytoplankton biomass in the California current system. Sci. Rep. 6, 27612 (2016).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

36.
Botsford, L. W., Lawrence, C. A., Dever, E. P., Hastings, A. & Largier, J. Wind strength and biological productivity in upwelling systems: An idealized study. Fish. Oceanogr. 12, 245–259 (2003).
Article  Google Scholar 

37.
García-Reyes, M., Largier, J. L. & Sydeman, W. J. Synoptic-scale upwelling indices and predictions of phyto- and zooplankton populations. Prog. Oceanogr. 120, 177–188 (2014).
ADS  Article  Google Scholar 

38.
Libralato, S., Coll, M., Tudela, S., Palomera, I. & Pranovi, F. Novel index for quantification of ecosystem effects of fishing as removal of secondary production. Mar. Ecol. Prog. Ser. https://doi.org/10.3354/meps07224 (2008).
Article  Google Scholar 

39.
Gasol, J. M., del Giorgio, P. A. & Duarte, C. M. Biomass distribution in marine planktonic communities. Limnol. Oceanogr. 42, 1353–1363 (1997).
ADS  CAS  Article  Google Scholar 

40.
Harfoot, M. B. J. et al. Emergent global patterns of ecosystem structure and function from a mechanistic general ecosystem model. PLoS Biol. 12, e1001841 (2014).
PubMed  PubMed Central  Article  CAS  Google Scholar 

41.
Wang, H., Morrison, W., Singh, A. & Weiss, H. General Mechanisms for Inverted Biomass Pyramids in Ecosystems. arXiv:0811.3657. [q-bio] (2008).

42.
Benoit-Bird, K. J. & Lawson, G. L. Ecological insights from pelagic habitats acquired using active acoustic techniques. Ann. Rev. Mar. Sci. 8, 463–490 (2016).
PubMed  Article  PubMed Central  Google Scholar 

43.
Alcaraz, M., Felipe, J., Grote, U., Arashkevich, E. & Nikishina, A. Life in a warming ocean: Thermal thresholds and metabolic balance of arctic zooplankton. J. Plankton Res. 36, 3–10 (2014).
Article  Google Scholar 

44.
Brochier, T. et al. Complex small pelagic fish population patterns arising from individual behavioral responses to their environment. Prog. Oceanogr. 164, 12–27 (2018).
ADS  Article  Google Scholar 

45.
Richardson, A., Schoeman, D., Richardson, A. J. & Schoeman, D. S. Climate impact on plankton ecosystems in the Northeast Atlantic. Science 305, 1609–1612 (2004).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

46.
Braham, C.-B. & Corten, A. Pelagic fish stocks and their response to fisheries and environmental variation in the Canary Current large marine ecosystem. Oceanographic and biological features in the Canary Current Large Marine Ecosystem 197–213 (2015).

47.
Ba, K. et al. Resilience of key biological parameters of the Senegalese flat sardinella to overfishing and climate change. PLoS ONE 11, e0156143 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

48.
Thiaw, M. et al. Effect of environmental conditions on the seasonal and inter-annual variability of small pelagic fish abundance off North-West Africa: The case of both Senegalese sardinella. Fish. Oceanogr. https://doi.org/10.1111/fog.12218 (2017).
Article  Google Scholar 

49.
Sarré, A. et al. Climate-driven shift of Sardinella aurita stock in Northwest Africa ecosystem as evidenced by robust spatial indicators [résumé]. In International conference ICAWA 2016 : extended book of abstract : the AWA project : ecosystem approach to the management of fisheries and the marine environment in West African waters (eds. Brehmer, P. et al.) 67–68 (SRFC/CSRP, 2017).

50.
Richardson, A. J. In hot water: Zooplankton and climate change. ICES J. Mar. Sci. 65, 279–295 (2008).
Article  Google Scholar 

51.
Beaugrand, G., Reid, P., Ibañez, F., Lindley, J. & Edwards, M. Reorganization of North Atlantic marine copepod biodiversity and climate. Science 296, 1692–1694 (2002).
ADS  CAS  PubMed  Article  Google Scholar 

52.
Lindley, J. A. & Daykin, S. Variations in the distributions of Centropages chierchiae and Temora stylifera (Copepoda: Calanoida) in the north-eastern Atlantic Ocean and western European shelf waters. ICES J. Mar. Sci. 62, 869–877 (2005).
Article  Google Scholar 

53.
Chassot, E. et al. Global marine primary production constrains fisheries catches. Ecol. Lett. 13, 495–505 (2010).
PubMed  Article  Google Scholar 

54.
Diogoul, N. et al. Fine-scale vertical structure of sound-scattering layers over an east border upwelling system and its relationship to pelagic habitat characteristics. Ocean Sci. 16, 65–81 (2020).
ADS  CAS  Article  Google Scholar 

55.
Brehmer, P. et al. Schooling behaviour of small pelagic fish: Phenotypic expression of independent stimuli. Mar. Ecol. Prog. Ser. 334, 263–272 (2007).
ADS  Article  Google Scholar 

56.
Munday, P. L., Jones, G. P., Pratchett, M. S. & Williams, A. J. Climate change and the future for coral reef fishes. Fish Fish. 9, 261–285 (2008).
Article  Google Scholar 

57.
Costello, J. H., Sullivan, B. K. & Gifford, D. J. A physical–biological interaction underlying variable phenological responses to climate change by coastal zooplankton. J. Plankton Res. 28, 1099–1105 (2006).
Article  Google Scholar 

58.
Garzke, J., Ismar, S. M. H. & Sommer, U. Climate change affects low trophic level marine consumers: Warming decreases copepod size and abundance. Oecologia 177, 849–860 (2015).
ADS  PubMed  Article  PubMed Central  Google Scholar 

59.
Horne, C. R., Hirst, A. G., Atkinson, D., Neves, A. & Kiørboe, T. A global synthesis of seasonal temperature–size responses in copepods. Glob. Ecol. Biogeogr. 25, 988–999 (2016).
Article  Google Scholar 

60.
Clark, C. W. & Levy, D. A. Diel vertical migrations by juvenile sockeye salmon and the antipredation window. Am. Nat. 131, 271–290 (1988).
Article  Google Scholar 

61.
Lampert, W. The adaptive significance of diel vertical migration of zooplankton. Funct. Ecol. 3, 21–27 (1989).
Article  Google Scholar 

62.
Hansson, S. Variation in hydroacoustic abundance of pelagic fish. Fish. Res. 16, 203–222 (1993).
Article  Google Scholar 

63.
Tiedemann, M. & Brehmer, P. Larval fish assemblages across an upwelling front: Indication for active and passive retention. Estuar. Coast. Shelf Sci. 187, 118–133 (2017).
ADS  Article  Google Scholar 

64.
CCLME. Analyse diagnostique transfrontalière du Grand écosystème marin du courant des Canaries 144 (2016).

65.
Bernhardt, J. R. & Leslie, H. M. Resilience to climate change in coastal marine ecosystems. Annu. Rev. Mar. Sci. 5, 371–392 (2013).
Article  Google Scholar 

66.
Baldé, B. S. et al. Variability of key biological parameters of round sardinella Sardinella aurita and the effects of environmental changes. J. Fish Biol. 94, 391–401 (2019).
PubMed  Article  PubMed Central  Google Scholar 

67.
Binet, D. Rôle possible d’une intensification des alizés sur le changement de répartition des sardines et sardinelles le long de la côte Ouest africaine. Aquat. Living Resour. 1, 115–132 (1988).
Article  Google Scholar 

68.
Berraho, A., Somoue, L., Hernández‐León, S. & Valdés, L. Zooplankton in the canary current large marine ecosystem. In Oceanographic and biological features in the Canary Current Large Marine Ecosystem Vol. 115, 183‐195 (IOC Technical Series, 2015).

69.
Ndour, I., Berraho, A., Fall, M., Ettahiri, O. & Sambe, B. Composition, distribution and abundance of zooplankton and ichthyoplankton along the Senegal-Guinea maritime zone (West Africa). Egypt. J. Aquat. Res. 44, 109–124 (2018).
Article  Google Scholar 

70.
Sarré, A., Krakstad, J.-O., Brehmer, P. & Mbye, E. M. Spatial distribution of main clupeid species in relation to acoustic assessment surveys in the continental shelves of Senegal and The Gambia. Aquat. Living Resour. 31, 9 (2018).
Article  Google Scholar 

71.
Foote, K. G., Knudsen, H. P., Vestnes, G., MacLennan, D. N. & Simmonds, E. J. Technical Report: ‘“Calibration of acoustic instruments for fish density estimation: A practical guide”’. J. Acoust. Soc. Am. 83, 831–832 (1987).
Google Scholar 

72.
Perrot, Y. et al. Matecho: An open-source tool for processing fisheries acoustics data. Acoust. Aust. 46, 241–248 (2018).
Article  Google Scholar 

73.
MacLennan, D. N., Fernandes, P. G. & Dalen, J. A consistent approach to definitions and symbols in fisheries acoustics. ICES J. Mar. Sci. 59, 365–369 (2002).
Article  Google Scholar 

74.
Jech, J. M., Lawson, G. L. & Lavery, A. C. Wideband (15–260 kHz) acoustic volume backscattering spectra of Northern krill (Meganyctiphanes norvegica) and butterfish (Peprilus triacanthus). ICES J. Mar. Sci. 74, 2249–2261 (2017).
Article  Google Scholar 

75.
Madureira, L. S. P., Everson, I. & Murphy, E. J. Interpretation of acoustic data at two frequencies to discriminate between Antarctic krill (Euphausia superba Dana) and other scatterers. J. Plankton Res. 15, 787–802 (1993).
Article  Google Scholar 

76.
Brehmer, P., Georgakarakos, S., Josse, E., Trygonis, V. & Dalen, J. Adaptation of fisheries sonar for monitoring schools of large pelagic fish: Dependence of schooling behaviour on fish finding efficiency. Aquat. Living Resour. 20, 377–384 (2007).
Article  Google Scholar 

77.
D’Elia, L. et al. A longitudinal study of the teacch program in different settings: The potential benefits of low intensity intervention in preschool children with autism spectrum disorder. J. Autism Dev. Disord. 44, 615–626 (2014).
PubMed  Article  Google Scholar 

78.
Zwolinski, J., Morais, A., Marques, V., Stratoudakis, Y. & Fernandes, P. G. Diel variation in the vertical distribution and schooling behaviour of sardine (Sardina pilchardus) off Portugal. ICES J. Mar. Sci. 64, 963–972 (2007).
Article  Google Scholar 

79.
Ayoubi, S. E. et al. Estimation of target strength of Sardina pilchardus and Sardinella aurita by theoretical approach. Fish. Sci. 82, 417–423 (2016).
Article  CAS  Google Scholar 

80.
Saunders, R. A., Fielding, S., Thorpe, S. E. & Tarling, G. A. School characteristics of mesopelagic fish at South Georgia. Deep Sea Res. Part I 81, 62–77 (2013).
Article  Google Scholar 

81.
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2016).

82.
Suppiah, R. & Hennessy, K. J. Trends in total rainfall, heavy rain events and number of dry days in Australia, 1910–1990. Int. J. Climatol. 18, 1141–1164 (1998).
Article  Google Scholar 

83.
Cotter, J. A selection of nonparametric statistical methods for assessing trends in trawl survey indicators as part of an ecosystem approach to fisheries management (EAFM). Aquat. Living Resour. 22, 173–185 (2009).
Article  Google Scholar  More

We have used a multi-scale modelling approach combining multiple modelling and inference types and techniques: generalised linear mixed-effect models (GLMMs), a mechanistic BTV transmission simulation model, marginal maximum likelihood inference, and direct calculation of resultant BTV risk predictions. We have parameterised our models with data from multiple existing sources (spatio-temporal climate time series, livestock density estimates, midge life process rates, BTV incidence time series), as well as data from our recent field study on the activity of midges at different latitudes in Europe34. The workflow and data sources of this study are summarised in Fig. 1.
Biting midge collection and identification
The method of biting midge sampling and identification was described in our earlier study34. In short, three habitat types were defined in which Culicoides specimens were collected: “farm”, “peri-urban”, and “wetland”. Traps were placed within a 50 m radius of cattle, a house, or waterbody, for farm, peri-urban, and wetland habitats respectively. Habitat types generally matched the classification of the CORINE European Land cover database46. Collections were performed in Sweden (surroundings of Linköping N58.410808, E15.621532), The Netherlands (surroundings of Wageningen N51.964795, E5.662898), and Italy (surroundings of San Benedetto del Tronto N42.949483, E13.878503). Onderstepoort Veterinary Institute black light traps were placed at three independent locations for each selected habitat type. Traps were at least 100 m apart to prevent overlap of the active trapping radius47. More details and exact trap locations can be found in Vogels et al.48 and Möhlmann et al.34. Collections were performed for six consecutive days in each month in all three countries, during the period from July 2014 to June 2015 except the winter months of December, January and February (and March for Sweden). Traps were emptied and repositioned at different locations every 24 h. Collections were sorted and stored in 70% ethanol at − 20 °C. Samples were identified to species level using the Interactive Identification Key for Culicoides (IIKC) developed by Mathieu et al.49.
Predictor variables for midge activity
Tinytag® meteorological data loggers (Gemini Data Loggers, Chichester, UK) were used to record local temperature and relative humidity every 30 min during the collection period from 17th July 2014 until 3rd July 2015. For each habitat in each country, one data logger was used (3 countries × 3 habitats). Additional meteorological data was collected from a number of sources: hourly wind velocities and local temperatures were available from the weather station closest to the trap location in Sweden (Swedish Meteorological and Hydrological Institute (SMHI) Linköping weather station), The Netherlands (Koninklijk Nederlands Meteorologisch Instituut (KNMI) weather stations “De Bilt”, “Deelen”, “Cabauw”, and “Volkel”), and Italy (San Benedetto del Tronto weather station). Finally, additional minimum, maximum, and mean daily temperatures along with precipitation and air pressure were sourced from the E-OBS European climate database50. This data was available at a daily temporal resolution and a spatial resolution of 0.25 degrees latitude and longitude. For climatic variables such as wind only one source per area could be used, whereas for temperature we had the opportunity to explore the most predictive of a number of data sources per country. Habitat effect was included by using the habitat types (farm, peri-urban, wetland) as a categorical predictor in the regression models.
BTV-seroprevalence survey data from Dutch sentinel herds
After the initial BTV outbreaks in 2006, the Dutch government decided to establish a sentinel network of dairy cattle herds in the winter of 2006/2007 to monitor the re-emergence of BTV-8 in 2007, using repeated milk ELISA testing35,36. For study purposes, The Netherlands was divided into 20 compartments based on geographic boundaries as proposed in Commission Decision 2005/393/EC. In each compartment, at least 10 randomly selected herds had to be sampled (with at least sixteen cows per herd) to obtain the required sample size. Herds were not necessarily completely BTV-8 seronegative at initial investigation, but cows designated for the sentinel program had to be BTV-8 seronegative at the moment of selection in May 2007. Therefore, dairy herds to act as sentinels for BTV incidence were selected, that had at least sixteen seronegative cows and at least 50 cattle in total. Monthly milk samples were collected from the sentinel cows in each herd unless prevalence had already reached 100%. The first round of monthly testing of sentinel cows was done in June 2007 and continued until January 2008. The monthly milk samples were tested at GD Animal Health, Deventer The Netherlands for antibodies to BTV-8 using a commercially available ELISA test. For further details on the sampling protocol and commercial ELISA see Santman-Berends et al.35,36.
Regression model for biting midge catches
A preliminary inference investigation using a hurdle negative binomial model to explain trap catches, inferred using Metropolis–Hastings MCMC, suggested that the excess of zero catches could be explained without zero inflation using the hurdle mechanism. We therefore used GLMM regression for its convenience and flexibility. Regression for the fixed effect coefficients and variance parameters of the random effects was performed via maximum likelihood using the Laplace approximation method implemented by the fitglme MATLAB® function. We followed the recommendation guidelines of Bolker et al.51 for using generalized linear models in the context of applied ecology, starting from a model with a full set of predictors and performed systematic model reduction using AICc to score model improvement (backwards model selection). AICc is a well-established information criterion for model selection since it is easy to calculate and interpret for GLMMs. The AICc difference between two models (ΔAICc) estimates the relative likelihood of the two models (with (sim {text{exp(}} – Delta {text{AICc}}/2)) as relative likelihood for the model with greater AICc). Moreover, the model selected by lowest AICc is also the model that would be selected by leave-one-out cross-validation52.
Candidate predictor variables for removal from the model were chosen by assessing which fixed-effect coefficients had the greatest P value (for the null hypothesis that the coefficient is zero) and which random effect coefficients had the smallest predicted standard deviation. This was followed by trialling the removal of either of the variables and removing the variable that reduced AICc the most, until no further reduction could be made. Several trials were made where we started from a number of different highly over-parameterised models, which all ended with the same best model.
For any given combination of predictor variables, the catches were assumed to be conditionally Poisson distributed, with the conditional mean for each collection defined by the log-link relationship,

$$lnleft( {mu_{lct} } right) = beta cdot X_{lt} + b_{l} cdot Z_{lt} + rho_{ct} + varepsilon_{lt}$$
(1)

where (beta) is the fixed effect regression coefficient that applies to all catches, with ({X}_{lt}) denoting the fixed effect predictors at trap location l on day t. ({b}_{l}sim Normal(0,Sigma )) are the location-grouped random effect regression coefficients with covariance matrix (Sigma), with ({Z}_{lt}) denoting the random effect predictors at location l on day t. The number of midges per catch was highly variable. Therefore, we included an independent random effect for each catch location and day ({varepsilon }_{lt}sim Normal(0,{sigma }_{varepsilon })). This corresponds to assuming that the number of midges per catch is Poisson log-normally distributed (a standard distribution for over-dispersed count data53). Spatial autocorrelation has been found in other midge catch regression studies38, so we also included an autocorrelation random effect grouped by country of trapping and day; ({rho }_{ct}sim Normal(0,{sigma }_{varrho })). This effect accounts for events influencing trap catch at a wider scale than just the location definition of each trap. See supplementary information for full model variables, regression coefficients and random effect variance estimates as well as AICc improvements from other models.
Connecting the regression model for midge captures to a midge biting prediction for herds
We connected the regression model for daily midge trap catch to a prediction of daily biting on cattle. In the BTV modelling literature it is common to assume that the vector-to-host ratio, and therefore the biting rate per livestock host, is proportional to the expected capture rate of a midge trap catch regression model28,54. In this study, we assumed that the biting rate per livestock host is proportional to the regression model presented in Results, with the major difference that we infer the scaling parameter that best fitted the serological data from Dutch sentinel herds. The location-group random effects in the regression model modelled unexplained spatial variation in average midge capture counts between trapping locations. We assumed that this unexplained spatial variation in midge abundance (as measured by trapping) mirrored the spatial variation in midge biting between different herds. Combining the proportionality assumption with the spatial variation assumption gave the biting rate among herds as,

$${text{Biting }},{text{rate }},{text{of }},{text{susceptible}},{text{ midges }},{text{per}},{text{ cattle }},{text{at }},{text{herd}},h,{text{on}},{text{ day}},t = xi mu_{ht} .$$
(2)

where ({mu }_{ht}) is the expected trap catch from the regression model given the climate condition local to herd h on day t and (xi) is a scaling parameter between the mean catch prediction and the biting rate prediction. In order to infer the proportionality factor between the catch model and the daily biting rate ((xi)), the outcomes of a dynamic transmission model for each herd were linked to the observed BTV seroprevalence data (see—“Mechanistic transmission model for BTV transmission within herds” section).
Mechanistic transmission model for BTV transmission within herds
We used a dynamic and mechanistic model of BTV transmission within herds, which was then matched to the data from the Dutch sentinel study. The dynamic BTV transmission model was formulated using disease compartments and rate-based transitions (see Keeling and Rohani55 for further details on this class of transmission model). In addition, it is in most respects similar to the model presented by Gubbins et al.27 in treating infectiousness amongst cattle, and latency amongst midges, as multi-stage processes that evolve deterministically. In particular, we follow Gubbins et al. in modelling the life processes of the infectious and latent midges (mortality, extrinsic incubation of BTV, and biting) as temperature dependent, and therefore varying daily with local background atmospheric temperature. (see Fig. 6 for a schematic diagram of the transmission model).
Figure 6

Schematic representation of the cattle herd level BTV transmission model. The population of cattle and infected biting midges are divided amongst discrete disease compartments. BTV latent midges (EM) enter the model at a rate proportional to the daily prediction of the catch model. The extrinsic incubation period for latent midges is modelled as a multi-stage process before midges become infectious (IM). Susceptible cattle (SC) become infectious cattle (IC) after a bite from infectious midges (IM), to become resistant cattle (RC) in time (also modelled as a multi-stage process; red box). Transitions are shown as solid lines, coloured according to their dependence on environmental variables: constant per-capita (black), daily mean temperature dependent (red), all predictor variables of capture model and the catch-to-bite scale parameter ({varvec{xi}}) (blue). Dotted lines indicate where the number of infected individuals in one species increases the incidence rate in the other species. Outcomes of the model are linked to observed cattle milk serology time series by the first two infectious stages for cattle with virus being undetectable by ELISA (blue box), whereas subsequent infectious stages and the recovered stage are detectable by ELISA (green box). The likelihood function for (xi) was inferred by marginalisation over the latent stochastic variables affecting model outcomes (e.g. herd-specific random effects, daily fluctuations in midge activity). Used images were available under open licence Creative Commons Deed CC0.

Full size image

The dynamic and mechanistic BTV transmission model used in this study describes the evolution of the numbers of susceptible, infected, and recovered cattle as well as latent and infectious midges for each herd (Fig. 6). Temperature-dependent midge bionomic rates were used for biting frequency of individual infectious midges56, the incubation rate of BTV within the midge57, and for the midge mortality58. The midge bionomic rates at each herd on each day were determined by the local mean temperature day according to the E-OBS climate dataset (see above—Predictor variables for midge activity). We modelled the incubation period of BTV within the midge vector as a ten-stage process, which is within the range of best fit models found in meta-analysis and laboratory studies of BTV incubation57 (see supplementary information for complete model details and literature estimates for rates).
The period during which BTV-infected cattle are detectable using an ELISA test (typically from 8–9 days post-infection onwards59) does not match the period during which the cattle are infectious (rapidly post-infection and then for an average of 20.6 days27). BTV-infected cattle can be in four states that are relevant to transmission modelling and their milk serology: 1) uninfected and susceptible to BTV, 2) infectious but undetectable by milk ELISA, 3) infectious and detectable by milk ELISA or, 4) non-infectious recovered from BTV but still milk ELISA positive. The BTV infectious period for cattle is usually modelled as a 5-stage process27, therefore it was convenient to model cattle in the first two stages of their infectious period (an average duration of 8.2 days) as infectious but undetectable. Cattle in the final three stages of the BTV infectious period are infectious and detectable (see Fig. 6 for a schematic representation of the BTV transmission and serology model).
The rate at which new susceptible midges arrived to bite cattle was assumed to be proportional to the expected number of midges from the regression model (({mu }_{ht})), with predictions of midge catches on each day and in each regional compartment of The Netherlands whilst the sentinel herd study was ongoing using the E-OBS historic climate records (see above—“Connecting the regression model for midge captures to a midge biting prediction for herds” section). The random effects in the catch model imply that ({mu }_{ht}) is a daily varying random variable, and that our transmission model is in the class of piecewise-deterministic Markov processes60. We assumed that the location-grouped random effects observed in the catch model became herd-grouped random effects for the biting model. In other words, we assumed that the high variance in midge catching between trapping location reflects high variance in midge biting between different cattle herd locations. Although an assumption, this would explain the highly variable intensity of BTV transmission observed between different herds in The Netherlands sentinel survey despite each herd experiencing a similar climate36. Because the herd locations were known only as geographic compartment occupancy, the daily local climate variables from the gridded E-OBS data used for predicting ({mu }_{ht}) were averaged over all spatial grids overlapping the herd’s geographic compartment. Accessing detailed and spatio-temporally resolved wind data across Europe was challenging, therefore we used the long-term average wind velocity of The Netherlands weather stations (see above—“Predictor variables for biting midge activity” section) as a constant predictor.
When simulating an outbreak of BTV within a herd, we first determined all relevant climatic predictors for the herd’s regional compartment and the daily temperature dependent midge bionomic rates. Second, we generated the herd-grouped and daily varying random effect coefficients which determined how biting at the herd from susceptible midges varied from a median prediction. Third, we solved the resultant deterministic BTV transmission model for each farm using the ode45 MATLAB® function (see Fig. 6 for an overview).
Inferring the catch-to-biting scale parameter from serological data
We inferred a maximum likelihood estimator for the trap-to-bite scale parameter by repeated simulation of the percentage of cattle detectable by milk ELISA tests herds in the sentinel herd network. For this we used the climatic conditions of The Netherlands in 2007, and at each simulation repeated redrawing the unobserved random effects for each herd and day. The average likelihood over many repeated simulations corresponds to a Monte Carlo estimate of the true marginal likelihood of the parameter (xi). Estimating the likelihood over a range of values of (xi) allowed the construction of a log-likelihood profile.
The stochastic elements of the piecewise-deterministic BTV transmission model were (i) the herd-grouped random coefficients (this modelled how biting varied from herd-to-herd) and (ii) the daily varying random effects (this modelled how biting varied from day-to-day). It is convenient to denote ({W}_{h}= ({{b}_{l}}^{(h)},{{rho }_{0}}^{(h)},{{rho }_{1}}^{(h)},{{rho }_{2}}^{(h)},…,{{epsilon }_{0}}^{(h)},{{epsilon }_{1}}^{(h)},{{epsilon }_{2}}^{(h)},…)) as the collection of all stochastic elements for herd h. For each simulation of the transmission model in each herd, we first drew ({W}_{h}) from their inferred distribution (see supplementary information for distribution parameters of best fitting model).
The likelihood of ({W}_{h}) and (xi) for each herd h was the chance of selecting the numbers of ELISA seroconverted cattle observed at the herd each month by The Netherlands sentinel study from the underlying distribution of ELISA detectable cattle implied by simulating the transmission model conditional on ((xi ,{W}_{h})),

$$L_{h} left( {xi ,W_{h} } right) = P({text{Serology }},{text{data }},{text{collected }},{text{at }},{text{herd }},h| xi ,W_{h} ).$$
(3)

Since we were not interested in inferring the specific values of ({W}_{h}) for each herd, they were treated as “nuisance” parameters. We inferred a maximum likelihood estimate, with confidence intervals, for (xi) by first estimating the marginal likelihood for (xi) (that is the likelihood after integrating over all possible values of the nuisance parameters) at each herd h,

$${L}_{h}(xi ) = int {L}_{h}left(xi ,wright)fleft(wright)dw.$$
(4)

where (f) is the density function for the distribution of random effects derived from the trapping model. The marginal log-likelihood function for the trap-to-bite scaling parameter, (l(xi )), for serological data over a number of herds, was then just the sum of the individual herd marginal log-likelihoods,

$$l(xi ) = sum_{h}log {L}_{h}left(xi right).$$
(5)

The herds we chose to contribute to the log-likelihood were those where BTV was found to be already present at the beginning of the study (see supplementary information for more details), to avoid making further assumptions about the introduction mechanism into the herd.
In practice, the log-likelihood was estimated for a profile of values of (xi) by simulating multiple realisations of ({W}_{h}) for each herd, that is we estimated (4) by Monte Carlo integration for (3) over a range of values of (xi), and interpolating between points with polynomial regression. The maximum likelihood estimator, ({xi }^{*}), was the maximizer of the marginal log-likelihood function presented in the main text along with confidence intervals derived by a standard comparison to the ({chi }^{2}) distribution (see methods section of King et al.61 for a brief but comprehensive introduction to maximum likelihood estimation using log-likelihood profiles in the context of inference for dynamical systems).
Calculating and mapping the herd reproductive ratio for bluetongue
The reproductive ratio for BTV will differ from day to day and across space. This reflects seasonality and variation in both climatic trends, and the population density of midges and livestock hosts. We approached estimating the reproductive ratio for BTV in the spirit of the case reproductive ratio62 using a technique already developed for midges spreading BTV37. That is, we calculated the expected number of secondary cases amongst hosts due to a host initially infected on each day t in each grid cell x whilst taking into account how the conditions for BTV transmission at location x changed after time t, and using the maximum likelihood model for midge biting. The size of each grid cell was determined by the resolution of the relevant climate data. We used the 0.25 degrees longitude and latitude grid resolution of the E-OBS climate datasets to map reproductive ratio estimates for Europe across space and time. Cattle and sheep densities across Europe were drawn from the livestock Geo-Wiki dataset63. The livestock Geo-Wiki datasets were more finely resolved (0.0083 degrees grid resolution) than the E-OBS climate datasets. We calculated the reproductive ratio for each grid cell x on each day t at the resolution of the E-OBS datasets. The cattle and sheep densities for this coarser grained grid were the average densities over the Geo-Wiki cells contained within the coarser grained grid.
The average number of secondary BTV cases amongst all hosts (cattle and sheep) given a host initially infected on day t and at grid cell x, is denoted ({R}^{(C)}(x,t)) for an initial infected cow and ({R}^{(S)}(x,t)) for an initial infected sheep. For both host species, the average number of secondary cases can be calculated by considering; how many days the host’s viraemia will last, the rate at which the host is bitten each day, the percentage of the biting midges that will become infected, how many of these biting midges are expected to survive their EIP to become actively infectious, and how many livestock will be successfully infected by those actively infectious midges. The methodology for combining these estimates using information about midge bionomic rates, EIP distribution, and the daily temperatures on each day after the initial host was infected has been developed by Brand and Keeling37.
In this study, we adapted the Brand-method for calculating the reproductive ratio to two species, and used the catch-to-biting scalar derived from comparison between the mechanistic transmission model and the herd sentinel serological survey. The cross-transmission between host species depends on how midge bites are distributed between cattle and sheep. We estimated the proportion of midge bites on cattle at grid cell x, ({phi }^{(C)}(x)), given the availability of sheep using a common relative preference model, e.g. Szmaragd et al.64,

$${phi }^{left(Cright)}left(xright)=frac{{N}^{left(Cright)}left(xright)}{{N}^{left(Cright)}left(xright)+ pi {N}^{left(Sright)}left(xright)}.$$
(6)

where ({N}^{(C)}(x)) and ({N}^{(S)}(x)) are, the local density of cattle and sheep at grid cell x. Parameter (pi) is a measure of the vector preference for sheep compared to cattle; (pi 1) preference for sheep. A relative biting study for sheep and cattle has revealed a preference for biting cattle30, from which we derived an estimate (pi =0.115) for use in this study. We combined ({R}^{(C)}(x,t)) and ({R}^{(S)}(x,t)) into a single reproductive ratio by calculating the leading eigenvalue of the next-generation matrix65,

$$Rleft(x,tright)= sqrt{{phi }^{left(Cright)}left(xright){R}^{left(Cright)}left(x,tright)+left(1-{phi }^{left(Cright)}left(xright)right){R}^{left(Sright)}left(x,tright)}.$$
(7)

An attractive feature of using the reproductive ratio as a measure of transmission intensity is its uncomplicated relationship with the persistence of transmission; if (Rle 1) then the infectious pathogen cannot persist. However, we expect that the biting rate, and therefore the reproductive ratio, will vary from herd-to-herd. From Eqs. (1) and (5) we see that the rate of biting from the susceptible midge population at each herd on each day depends on the random coefficients, ({{b}_{l}}^{(h)}), and daily varying random effects,({{rho }_{ct}}^{(h)}) and ({{varepsilon }_{t}}^{(h)}),

$$text{Biting rate of susceptible midges per cattle at herd }htext{ on day }tpropto expleft({{b}_{l}}^{left(hright)}cdot {Z}_{lt} +{{rho }_{ct}}^{left(hright)} + {{varepsilon }_{t}}^{left(hright)} right).$$
(8)

The daily varying random effects ((rho) and (epsilon)) are averaged over our estimates for ({R}^{(C)}) and ({R}^{(S)}) (this can be achieved analytically; see supplementary information for further details), and therefore our estimate of the reproductive ratio does not depend on daily fluctuations in midge activity. However, variation in the herd-grouped random coefficients indicated systematic differences in midge activity between herds that will not ‘average out’ over time. The distribution of ({{b}_{l}}^{(h)}) therefore implied a distribution of biting rates for herds within each grid cell on each day, and therefore a distribution of values of (R) for herds in each grid cell and on each day.
We present the distribution of (R) for herds by considering the reproductive ratio that would be calculated if the random variable in Eq. (8), ({{b}_{l}}^{(h)}cdot {Z}_{lt}), took its pth percentile value every day, denoting this reproductive ratio, ({R}_{p}(x,t)). ({R}_{p}(x,t)) estimates the reproductive ratio that p% of herd reproductive ratios are less than in grid cell x on day t. Also, we numerically invert the threshold relationship to find the percentage value, ({varvec{P}}), such that ({R}_{p}(x,t)) satisfies the threshold quantity,

$${varvec{P}}left(x,tright)={1-p in [mathrm{0,1}] | {R}_{p}(x,t) = 1 }.$$
(9)

({varvec{P}}(x,t)) is therefore an estimate of the percentage of herds that could have a multiplying BTV outbreak in grid cell x if BTV was introduced on day t.
To enable spatially extending risk predictions for BTV across Europe certain assumptions about how the midge biting model could be interpolated between the trap locations were necessary. The best regression model for midge catches found significantly different seasonality at the Italian catch locations compared to Sweden and The Netherlands. We assumed that day length was a determinant of midge seasonality and overwintering at different latitudes; for example, the first day of the year with a day length shorter than 9 h has been associated with the onset of overwintering in UK midges66. We found no midges on days with less than 8.5 h of daylight when trapping, although catching was not attempted outside of March-November so the number of samples was small. Therefore, at latitudes where no day is ever shorter than 8.5 h (lower latitudes than 46oN, which is more or less the border of Switzerland and Italy) we used Italian seasonality to predict midge biting. At latitudes that have at least one day shorter than 8 h (higher latitudes than 49oN) we used the Swedish and Dutch midge seasonality. Between these two latitudes, a linear interpolation between the predictions of the two seasonal models was applied.
All map images were generated using MATLAB® contourf function. Both the livestock density and E-OBS datasets included grid cells that contain only water, these cells were coloured white.
The MATLAB® code for generating the reproductive ratio estimates is available at a github repository (https://github.com/SamuelBrand1/BTV-European-Reproductive-Ratios). More

1.
Wood, T. J. & Goulson, D. The environmental risks of neonicotinoid pesticides: a review of the evidence post 2013. Environ. Sci. Pollut. Res. Int. 24, 17285–17325. https://doi.org/10.1007/s11356-017-9240-x (2017).
CAS  Article  PubMed  PubMed Central  Google Scholar 
2.
Wagner, D. L. Insect declines in the anthropocene. Annu. Rev. Entomol. 65, 457–480. https://doi.org/10.1146/annurev-ento-011019-025151 (2020).
CAS  Article  PubMed  Google Scholar 

3.
Klein, A.-M. et al. Importance of pollinators in changing landscapes for world crops. Proc. R. Soc. B Biol. Sci. 274, 303–313. https://doi.org/10.1098/rspb.2006.3721 (2007).
Article  Google Scholar 

4.
Gallai, N., Salles, J.-M., Settele, J. & Vaissière, B. E. Economic valuation of the vulnerability of world agriculture confronted with pollinator decline. Ecol. Econ. 68, 810–821. https://doi.org/10.1016/j.ecolecon.2008.06.014 (2009).
Article  Google Scholar 

5.
Hallmann, C. A. et al. More than 75 percent decline over 27 years in total flying insect biomass in protected areas. PLoS ONE 12, e0185809. https://doi.org/10.1371/journal.pone.0185809 (2017).
CAS  Article  PubMed  PubMed Central  Google Scholar 

6.
Goulson, D. The insect apocalypse, and why it matters. Curr. Biol. 29, R967–R971. https://doi.org/10.1016/j.cub.2019.06.069 (2019).
CAS  Article  PubMed  Google Scholar 

7.
Casida, J. E. & Durkin, K. A. Neuroactive insecticides: targets, selectivity, resistance, and secondary effects. Annu. Rev. Entomol. 58, 99–117. https://doi.org/10.1146/annurev-ento-120811-153645 (2013).
CAS  Article  PubMed  Google Scholar 

8.
Popp, J., Pető, K. & Nagy, J. Pesticide productivity and food security. A review. . Agron. Sustain. Dev. 33, 243–255. https://doi.org/10.1007/s13593-012-0105-x (2013).
Article  Google Scholar 

9.
Casida, J. E. Neonicotinoids and other insect nicotinic receptor competitive modulators: progress and prospects. Annu. Rev. Entomol. 63, 125–144. https://doi.org/10.1146/annurev-ento-020117-043042 (2018).
CAS  Article  PubMed  Google Scholar 

10.
Matsuda, K., Ihara, M. & Sattelle, D. B. Neonicotinoid insecticides: molecular targets, resistance, and toxicity. Annu. Rev. Pharmacol. Toxicol. 60, 241–255. https://doi.org/10.1146/annurev-pharmtox-010818-021747 (2020).
CAS  Article  PubMed  Google Scholar 

11.
Goulson, D. REVIEW: an overview of the environmental risks posed by neonicotinoid insecticides. J. Appl. Ecol. 50, 977–987. https://doi.org/10.1111/1365-2664.12111 (2013).
Article  Google Scholar 

12.
Eng, M. L., Stutchbury, B. J. M. & Morrissey, C. A. A neonicotinoid insecticide reduces fueling and delays migration in songbirds. Science 365, 1177. https://doi.org/10.1126/science.aaw9419 (2019).
ADS  CAS  Article  PubMed  Google Scholar 

13.
Yamamuro, M. et al. Neonicotinoids disrupt aquatic food webs and decrease fishery yields. Science (New York, N.Y.) 366, 620. https://doi.org/10.1126/science.aax3442 (2019).
ADS  CAS  Article  Google Scholar 

14.
Han, W., Tian, Y. & Shen, X. Human exposure to neonicotinoid insecticides and the evaluation of their potential toxicity: an overview. Chemosphere 192, 59–65. https://doi.org/10.1016/j.chemosphere.2017.10.149 (2018).
ADS  CAS  Article  PubMed  Google Scholar 

15.
Nauen, R., Ebbinghaus-Kintscher, U., Salgado, V. L. & Kaussmann, M. Thiamethoxam is a neonicotinoid precursor converted to clothianidin in insects and plants. Pestic. Biochem. Physiol. 76, 55–69. https://doi.org/10.1016/S0048-3575(03)00065-8 (2003).
CAS  Article  Google Scholar 

16.
EFSA. Peer review of the pesticide risk assessment of the active substance thiacloprid. Eur. Food Saf. Auth. J. 17, e05595 https://doi.org/10.2903/j.efsa.2019.5595 (2019).
Article  Google Scholar 

17.
Nicholls, E. et al. Monitoring neonicotinoid exposure for bees in rural and peri-urban areas of the U.K. during the transition from pre- to post-moratorium. Environ. Sci. Technol. 52, 9391–9402. https://doi.org/10.1021/acs.est.7b06573 (2018).
ADS  CAS  Article  PubMed  Google Scholar 

18.
Wintermantel, D. et al. Neonicotinoid-induced mortality risk for bees foraging on oilseed rape nectar persists despite EU moratorium. Sci. Total Environ. 704, 135400. https://doi.org/10.1016/j.scitotenv.2019.135400 (2020).
ADS  CAS  Article  PubMed  Google Scholar 

19.
Cressey, D. Fears for bees as UK lifts insecticide ban. Nature News. https://www.nature.com/news/fears-for-bees-as-uk-lifts-insecticide-ban-1.18052 (2015).

20.
Lamsa, J., Kuusela, E., Tuomi, J., Juntunen, S. & Watts, P. C. Low dose of neonicotinoid insecticide reduces foraging motivation of bumblebees. Proc. R. Soc. B 285, 20180506 https://doi.org/10.1098/rspb.2018.0506 (2018).
Article  PubMed  PubMed Central  Google Scholar 

21.
Tasman, K., Rands, S. A. & Hodge, J. J. The neonicotinoid insecticide imidacloprid disrupts bumblebee foraging rhythms and sleep. iScience 23, 101827 https://doi.org/10.2139/ssrn.3586989 (2020).
Article  PubMed  PubMed Central  Google Scholar 

22.
Palmer, M. J. et al. Cholinergic pesticides cause mushroom body neuronal inactivation in honeybees. Nat. Commun. 4, 1634. https://doi.org/10.1038/ncomms2648 (2013).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

23.
Aso, Y. et al. Mushroom body output neurons encode valence and guide memory-based action selection in Drosophila. eLife 3, e04580. https://doi.org/10.7554/eLife.04580 (2014).
Article  PubMed  PubMed Central  Google Scholar 

24.
Barnstedt, O. et al. Memory-relevant mushroom body output synapses are cholinergic. Neuron 89, 1237–1247. https://doi.org/10.1016/j.neuron.2016.02.015 (2016).
CAS  Article  PubMed  PubMed Central  Google Scholar 

25.
Helfrich-Förster, C. Sleep in insects. Annu. Rev. Entomol. 63, 69–86. https://doi.org/10.1146/annurev-ento-020117-043201 (2018).
CAS  Article  PubMed  Google Scholar 

26.
Peng, Y. C. & Yang, E. C. Sublethal dosage of imidacloprid reduces the microglomerular density of honey bee mushroom bodies. Sci. Rep. 6, 19298. https://doi.org/10.1038/srep19298 (2016).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

27.
Smith, D. B. et al. Insecticide exposure during brood or early-adult development reduces brain growth and impairs adult learning in bumblebees. Proc. R. Soc. B 287, 20192442. https://doi.org/10.1098/rspb.2019.2442 (2020).
CAS  Article  PubMed  Google Scholar 

28.
Andrione, M., Vallortigara, G., Antolini, R. & Haase, A. Neonicotinoid-induced impairment of odour coding in the honeybee. Sci. Rep. 6, 38110. https://doi.org/10.1038/srep38110 (2016).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

29.
Chouhan, N. S., Wolf, R., Helfrich-Förster, C. & Heisenberg, M. Flies remember the time of day. Curr. Biol. 25, 1619–1624. https://doi.org/10.1016/j.cub.2015.04.032 (2015).
CAS  Article  PubMed  Google Scholar 

30.
Flyer-Adams, J. G. et al. Regulation of olfactory associative memory by the circadian clock output signal Pigment-dispersing factor (PDF). J. Neurosci. 40, 9066–9077https://doi.org/10.1523/JNEUROSCI.0782-20.2020 (2020).
Article  PubMed  PubMed Central  Google Scholar 

31.
Zwaka, H. et al. Context odor presentation during sleep enhances memory in honeybees. Curr. Biol. 25(21), 869–2874. https://doi.org/10.1016/j.cub.2015.09.069 (2015).
CAS  Article  Google Scholar 

32.
Seugnet, L., Suzuki, Y., Donlea, J. M., Gottschalk, L. & Shaw, P. J. Sleep deprivation during early-adult development results in long-lasting learning deficits in adult Drosophila. Sleep 34, 137–146. https://doi.org/10.1093/sleep/34.2.137 (2011).
Article  PubMed  PubMed Central  Google Scholar 

33.
Tackenberg, M. C. et al. Neonicotinoids disrupt circadian rhythms and sleep in honey bees. Sci. Rep. 10, 17929. https://doi.org/10.1038/s41598-020-72041-3 (2020).
CAS  Article  PubMed  PubMed Central  Google Scholar 

34.
Helfrich-Forster, C. et al. The extraretinal eyelet of Drosophila: development, ultrastructure, and putative circadian function. J. Neurosci. 22, 9255–9266. https://doi.org/10.1523/JNEUROSCI.22-21-09255.2002 (2002).
Article  PubMed  PubMed Central  Google Scholar 

35.
Muraro, N. I. & Ceriani, M. F. Acetylcholine from visual circuits modulates the activity of arousal neurons in Drosophila. J. Neurosci. 35, 16315. https://doi.org/10.1523/JNEUROSCI.1571-15.2015 (2015).
CAS  Article  PubMed  PubMed Central  Google Scholar 

36.
McCarthy, E. V. et al. Synchronized bilateral synaptic inputs to Drosophila melanogaster neuropeptidergic rest/arousal neurons. J. Neurosci. 31, 8181–8193. https://doi.org/10.1523/jneurosci.2017-10.2011 (2011).
CAS  Article  PubMed  PubMed Central  Google Scholar 

37.
Parisky, K. M. et al. PDF cells are a GABA-responsive wake-promoting component of the Drosophila sleep circuit. Neuron 60, 672–682. https://doi.org/10.1016/j.neuron.2008.10.042 (2008).
CAS  Article  PubMed  PubMed Central  Google Scholar 

38.
Ly, S., Pack, A. I. & Naidoo, N. The neurobiological basis of sleep: Insights from Drosophila. Neurosci. Biobehav. Rev. 87, 67–86. https://doi.org/10.1016/j.neubiorev.2018.01.015 (2018).
Article  PubMed  PubMed Central  Google Scholar 

39.
Wegener, C., Hamasaka, Y. & Nassel, D. R. Acetylcholine increases intracellular Ca2+ via nicotinic receptors in cultured PDF-containing clock neurons of Drosophila. J. Neurophysiol. 91, 912–923. https://doi.org/10.1152/jn.00678.2003 (2004).
CAS  Article  PubMed  Google Scholar 

40.
Renn, S. C., Park, J. H., Rosbash, M., Hall, J. C. & Taghert, P. H. A pdf neuropeptide gene mutation and ablation of PDF neurons each cause severe abnormalities of behavioral circadian rhythms in Drosophila. Cell 99, 791–802. https://doi.org/10.1016/s0092-8674(00)81676-1 (1999).
CAS  Article  PubMed  Google Scholar 

41.
Schlichting, M., Menegazzi, P., Rosbash, M. & Helfrich-Förster, C. A distinct visual pathway mediates high-intensity light adaptation of the circadian clock in Drosophila. J. Neurosci. 39, 1621. https://doi.org/10.1523/JNEUROSCI.1497-18.2018 (2019).
CAS  Article  PubMed  PubMed Central  Google Scholar 

42.
Lelito, K. & Shafer, O. Reciprocal cholinergic and GABAergic modulation of the small ventrolateral pacemaker neurons of Drosophila’s circadian clock neuron network. J. Neurophysiol. 107, 2096–2108. https://doi.org/10.1152/jn.00931.2011 (2012).
CAS  Article  PubMed  PubMed Central  Google Scholar 

43.
Nitabach, M. N. et al. Electrical hyperexcitation of lateral ventral pacemaker neurons desynchronizes downstream circadian oscillators in the fly circadian circuit and induces multiple behavioral periods. J. Neurosci. 26, 479–489. https://doi.org/10.1523/jneurosci.3915-05.2006 (2006).
CAS  Article  PubMed  PubMed Central  Google Scholar 

44.
Cao, G. & Nitabach, M. N. Circadian control of membrane excitability in Drosophila melanogaster lateral ventral clock neurons. J. Neurosci. 28, 6493–6501. https://doi.org/10.1523/jneurosci.1503-08.2008 (2008).
CAS  Article  PubMed  PubMed Central  Google Scholar 

45.
Fernández, M. P., Berni, J. & Ceriani, M. F. Circadian remodeling of neuronal circuits involved in rhythmic behavior. PLoS Biol. 6, e69. https://doi.org/10.1371/journal.pbio.0060069 (2008).
CAS  Article  PubMed  PubMed Central  Google Scholar 

46.
Park, J. H. et al. Differential regulation of circadian pacemaker output by separate clock genes in Drosophila. Proc. Natl. Acad. Sci. 97, 3608. https://doi.org/10.1073/pnas.97.7.3608 (2000).
ADS  CAS  Article  PubMed  Google Scholar 

47.
Martelli, F. et al. Low doses of the neonicotinoid insecticide imidacloprid induce ROS triggering neurological and metabolic impairments in Drosophila. Proc. Natl. Acad. Sci. 117(41), 25840–25850. https://doi.org/10.1073/pnas.2011828117 (2020).
CAS  Article  PubMed  Google Scholar 

48.
Numata, H., Miyazaki, Y. & Ikeno, T. Common features in diverse insect clocks. Zool. Lett. 1, 10. https://doi.org/10.1186/s40851-014-0003-y (2015).
Article  Google Scholar 

49.
Farris, S. & Sinakevitch, I. Development and evolution of the insect mushroom bodies: towards the understanding of conserved developmental mechanisms in a higher brain center. Arthropod Struct. Dev. 32, 79–101. https://doi.org/10.1016/S1467-8039(03)00009-4 (2003).
Article  PubMed  Google Scholar 

50.
Jones, A. K. & Sattelle, D. B. In Insect Nicotinic Acetylcholine Receptors (ed Thany, S. H.) 25–43 (Springer New York, 2010).

51.
Homem, R. A. et al. Evolutionary trade-offs of insecticide resistance—the fitness costs associated with target-site mutations in the nAChR of Drosophila melanogaster. Mol. Ecol. 29, 2661–2675. https://doi.org/10.1111/mec.15503 (2020).
CAS  Article  PubMed  PubMed Central  Google Scholar 

52.
Blacquiere, T., Smagghe, G., van Gestel, C. A. & Mommaerts, V. Neonicotinoids in bees: a review on concentrations, side-effects and risk assessment. Ecotoxicology 21, 973–992. https://doi.org/10.1007/s10646-012-0863-x (2012).
CAS  Article  PubMed  PubMed Central  Google Scholar 

53.
Stanley, D. A. & Raine, N. E. Bumblebee colony development following chronic exposure to field-realistic levels of the neonicotinoid pesticide thiamethoxam under laboratory conditions. Sci. Rep. 7, 8005. https://doi.org/10.1038/s41598-017-08752-x (2017).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

54.
Whitehorn, P. R., O’Connor, S., Wackers, F. L. & Goulson, D. Neonicotinoid pesticide reduces bumble bee colony growth and queen production. Science 336, 351. https://doi.org/10.1126/science.1215025 (2012).
ADS  CAS  Article  PubMed  Google Scholar 

55.
Williamson, S. M., Willis, S. J. & Wright, G. A. Exposure to neonicotinoids influences the motor function of adult worker honeybees. Ecotoxicology 23, 1409–1418. https://doi.org/10.1007/s10646-014-1283-x (2014).
CAS  Article  PubMed  PubMed Central  Google Scholar 

56.
Wright, G. A., Softley, S. & Earnshaw, H. Low doses of neonicotinoid pesticides in food rewards impair short-term olfactory memory in foraging-age honeybees. Sci. Rep. 5, 15322. https://doi.org/10.1038/srep15322 (2015).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

57.
Malik, B. R. & Hodge, J. J. Drosophila adult olfactory shock learning. J. Vis. Exp. 90, 50107. https://doi.org/10.3791/50107 (2014).
CAS  Article  Google Scholar 

58.
Hodge, J. J. & Stanewsky, R. Function of the Shaw potassium channel within the Drosophila circadian clock. PLoS ONE 3, e2274–e2274. https://doi.org/10.1371/journal.pone.0002274 (2008).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

59.
Moffat, C. et al. Neonicotinoids target distinct nicotinic acetylcholine receptors and neurons, leading to differential risks to bumblebees. Sci. Rep. 6, 24764. https://doi.org/10.1038/srep24764 (2016).
ADS  Article  PubMed  PubMed Central  Google Scholar 

60.
Busto, G. U., Cervantes-Sandoval, I. & Davis, R. L. Olfactory learning in Drosophila. Physiology 25, 338–346. https://doi.org/10.1152/physiol.00026.2010 (2010).
CAS  Article  PubMed  Google Scholar 

61.
Lyons, L. C. & Roman, G. Circadian modulation of short-term memory in Drosophila. Learn. Mem. 16, 19–27. https://doi.org/10.1101/lm.1146009 (2009).
Article  PubMed  PubMed Central  Google Scholar 

62.
Depetris-Chauvin, A. et al. Adult-specific electrical silencing of pacemaker neurons uncouples the molecular oscillator from circadian outputs. Curr. Biol. 21, 1783–1793. https://doi.org/10.1016/j.cub.2011.09.027 (2011).
CAS  Article  PubMed  PubMed Central  Google Scholar 

63.
Baz, E.-S., Wei, H., Grosshans, J. & Stengl, M. Calcium responses of circadian pacemaker neurons of the cockroach Rhyparobia maderae to acetylcholine and histamine. J. Comp. Physiol. A. 199, 365–374. https://doi.org/10.1007/s00359-013-0800-3 (2013).
CAS  Article  Google Scholar 

64.
Sheeba, V. et al. Large ventral lateral neurons modulate arousal and sleep in Drosophila. Curr. Biol. 18, 1537–1545. https://doi.org/10.1016/j.cub.2008.08.033 (2008).
CAS  Article  PubMed  PubMed Central  Google Scholar 

65.
Thany, S. H. Insect Nicotinic Acetylcholine Receptors (Springer, New York, 2011).
Google Scholar 

66.
Gill, R. J. & Raine, N. E. Chronic impairment of bumblebee natural foraging behaviour induced by sublethal pesticide exposure. Funct. Ecol. 28, 1459–1471. https://doi.org/10.1111/1365-2435.12292 (2014).
Article  Google Scholar 

67.
Bloch, G., Bar-Shai, N., Cytter, Y. & Green, R. Time is honey: circadian clocks of bees and flowers and how their interactions may influence ecological communities. Phil. Trans. R. Soc. B 372, 20160256. https://doi.org/10.1098/rstb.2016.0256 (2017).
CAS  Article  Google Scholar 

68.
van Alphen, B., Yap, M. H. W., Kirszenblat, L., Kottler, B. & van Swinderen, B. A dynamic deep sleep stage in Drosophila. J. Neurosci. 33, 6917. https://doi.org/10.1523/JNEUROSCI.0061-13.2013 (2013).
CAS  Article  PubMed  PubMed Central  Google Scholar 

69.
Buhl, E., Higham, J. P. & Hodge, J. J. L. Alzheimer’s disease-associated tau alters Drosophila circadian activity, sleep and clock neuron electrophysiology. Neurobiol. Dis. 130, 104507. https://doi.org/10.1016/j.nbd.2019.104507 (2019).
CAS  Article  PubMed  Google Scholar 

70.
Levine, J. D., Funes, P., Dowse, H. B. & Hall, J. C. Signal analysis of behavioral and molecular cycles. BMC Neurosci. 3, 1. https://doi.org/10.1186/1471-2202-3-1 (2002).
Article  PubMed  PubMed Central  Google Scholar 

71.
Faville, R., Kottler, B., Goodhill, G. J., Shaw, P. J. & van Swinderen, B. How deeply does your mutant sleep? Probing arousal to better understand sleep defects in Drosophila. Sci. Rep. 5, 8454. https://doi.org/10.1038/srep08454 (2015).
CAS  Article  PubMed  PubMed Central  Google Scholar 

72.
Donelson, N. C. et al. High-resolution positional tracking for long-term analysis of Drosophila sleep and locomotion using the “tracker” program. PLoS ONE 7, e37250. https://doi.org/10.1371/journal.pone.0037250 (2012).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

73.
Schindelin, J. et al. Fiji: an open-source platform for biological-image analysis. Nat. Methods 9, 676. https://doi.org/10.1038/nmeth.2019 (2012).
CAS  Article  Google Scholar  More

Determination of the mtDNA D-loop haplotypes of indigenous chicken breeds and red junglefowl in Thailand
We determined the nucleotide sequences of the 780 bp fragments of the mtDNA D-loop region, including the hypervariable segment I, in 125 individuals from 10 indigenous chicken breeds (a list of breeds is shown in Table 1), and 279 red junglefowls from two subspecies (G. g. gallus and G. g. spadiceus) within 12 populations in Thailand. A total of 44 haplotypes with 62 variable sites, consisting of 26 singletons and 36 parsimony informative sites were identified (Supplementary Tables S1, S2; accession No. LC542982 to LC543385). Table 2 summarizes the details of the haplotypes found in the 10 indigenous chicken breeds and 12 red junglefowl populations, and Fig. 1 shows the composition of each haplogroup of indigenous chickens and two subspecies of red junglefowl. Forty-four haplotypes were temporally classified into eight common haplogroups; A, B, C, D, E, F, H, and J (Supplementary Figs. S1, S2), according to Liu et al.4 and Miao et al.5. In the present study, we treated haplogroups C and D as one unit (CD) as they were not clearly separated. In addition, haplogroup J was closely related to haplogroup CD. Haplogroups A, B, and E were predominant in the Thai indigenous chickens (Table 2), and their frequencies were almost the same (Fig. 1). Haplogroup CD was predominant in the G. g. spadiceus population, but rare in indigenous chickens. In the G. g. gallus population, the haplogroups B, CD, and E were detected at almost the same frequency; however, haplogroup A was not detected. The frequency of haplogroup J, which was mainly found in the Si Sa Ket population, was much higher in G. g. gallus compared with indigenous chicken and G. g. spadiceus populations (Fig. 1). BT, NK-W, NK-B, LHK, CH, PHD, Decoy, fighting chicken, and seven red junglefowl populations (Huai Sai [Ggg], Huai Sai [Ggs], Sa Kaeo, Chanthaburi, Khao Kho, Chaiyaphum, and Khok Mai Rua) each exhibited breed- or population-specific haplotypes (Supplementary Table S2): Hap_04 (BT); Hap_07 and 08 (NK-W); Hap_09 (NK-B); Hap_10 and Hap_11 (LHK); Hap_14 (CH); Hap_17 (PHD); Hap_18 to 20 (Decoy); Hap_29 to 33 (fighting chicken); Hap_21 (Huai Sai [Ggg]); Hap_22 (Huai Sai [Ggs]); Hap_24 and 25 (Sa Kaeo); Hap_27 and 28 (Chanthaburi); Hap_36 (Khao Kho); Hap_38 and 39 (Chaiyaphum); and Hap_43 and 44 (Khok Mai Rua). Twenty-nine out of 44 D-loop haplotypes which had not been previously deposited in GenBank, were newly identified in the present study.
Table 1 List of indigenous chicken breeds and red junglefowl populations examined in the present study.
Full size table

Table 2 Summary of haplogroups of mtDNA D-loop sequences and haplotypes that were found in 10 indigenous chicken breeds and 12 populations of two Gallus gallus subspecies in Thailand and their distribution.
Full size table

Figure 1

Composition ratio of haplogroups of the mtDNA D-loop sequences in chickens indigenous to Thailand, G. g. spadiceus, and G. g. gallus. The haplogroup names were conformed to those described by Miao et al.5 The numbers in parentheses indicate the number of individuals examined.

Full size image

The topologies of the Bayesian tree and the maximum-likelihood (ML) tree based on the HKY + G + I model of evolution, which were selected as the best-fit substitution model, were fundamentally similar. Although the Bayesian posterior probability of the internal nodes and the ML bootstrap values were relatively low due to the short internal branches (multifurcations) of the phylogenetic trees, the haplogroups A, B, F–I, K, Y, and Z were supported by a Bayesian posterior probability of greater than 0.97 (Supplementary Figs. S1, S2). Both The trees revealed that the D-loop sequences obtained in this study could be classified into six haplogroups: A, B, CD, E, F, and J, and a complex group of rare haplogroups (H, I, K, W, and X), except for an unclassified haplotype, Hap_38.
Six haplotypes from seven indigenous and two red junglefowl populations (LHK, CH, PHD, KP, BT, Decoy, fighting chicken, Huai Sai [Ggs], and Petchaburi) belonged to haplogroup A (Fig. 2a). Seven haplotypes from seven indigenous chicken breeds and five red junglefowl populations (LHK, CH, PHD, KP, Decoy, fighting chicken, DT, Sa Kaeo, Huai Sai [Ggg], Huai Sai [Ggs], Khao Kho, and Petchaburi) were classified into haplogroup B (Fig. 2b). Haplogroup CD contained 12 haplotypes, which were identified in two Thai indigenous chicken breeds (PHD and fighting chicken) and seven red junglefowl populations (Chanthaburi, Khok Mai Rua, Chaing Rai, Huai Sai [Ggs], Khao Kho, Chaiyaphum, and Huai Yang Pan) (Fig. 2c). Eight haplotypes in four indigenous chicken breeds (CH, BT, NK-W, and NK-B) and four red junglefowl populations (Roi Et, Khok Mai Rua, Chaiyaphum, and Huai Yang Pan) belonged to haplogroup E (Fig. 2e). Haplogroup F contained one haplotype, which was only found in two indigenous chicken breeds (LHK and PHD) (Fig. 2f). Eight haplotypes of haplogroup J (Hap_10, Hap_11, Hap_24 to 26, Hap_34, Hap_40, and Hap_42) were found in two indigenous chicken breeds (LHK and fighting chicken) and seven red junglefowl populations (Sa Kaeo, Chabthaburi, Si Saket, Roi Et, Khok Mai Rua, Chaing Rai, and Huai Yang Pan) (Fig. 2c). Only one haplotype of haplogroup H (Hap_05) was detected in two indigenous chicken breeds (BT and fighting chicken) (Fig. 2d). Hap_38, which was found in three individuals of the Chaiyaphum population, did not belong to any known haplogroups; however, the haplotype was more closely related to haplogroup CD than the other haplogroups (Fig. 2c; Supplementary Figs. S1, S2).
Figure 2

Locations of mtDNA D-loop haplotypes of Thai red junglefowl and indigenous chicken populations in the global chicken population network. (a) Haplogroup A. (b) Haplogroup B. (c) Haplogroups CD, Y, Z, J, and an unclassified haplotype, Hap_38. (d) Haplogroups H, I, K, X, and W. (e) Haplogroup E. (f) Haplogroup F. Haplotypes that were found in the present study and representative haplotypes reported by Miao et al.5 are shown by magenta and yellow circles, respectively. Black nodes are the inferred intermediate haplotypes. The number of bars on the lines, which link haplotypes, represent the number of nucleotide substitutions that occurred between the haplotypes for comparison.

Full size image

Divergence times for each haplotype were determined using BEAST analysis (Supplementary Table S2) and were 0.24–0.45 kilo years ago (KYA) for haplogroup A, 0.15–0.39 KYA for haplogroup B, and 0.14–0.37 KYA for haplogroup CD; the haplotypes of haplogroup E exhibited a wide range of divergent times, ranging from 0.12 to 0.70 KYA (0.41 KYA on average). One haplotype in haplogroups F and H had possibly diverged at 0.33 and 0.34 KYA, respectively. The divergence times of haplotypes in haplogroup J ranged from 0.10 to 0.60 KYA. Hap_38 exhibited a markedly earlier divergence time, which was estimated to be approximately 12,000 years ago (Supplementary Table S2; Supplementary Fig. S1).
Genetic diversity of mtDNA D-loop sequences
The number of D-loop haplotypes in each population (H) ranged from 1 (G. g. gallus population at Huai Sai and Si Sa Ket) to 10 (fighting chicken) (Table 3). Among the Thai indigenous chicken breeds, LHK, CH, PHD, KP, BT, Decoy and fighting chicken exhibited relatively higher genetic diversity (pi, 0.005 for KP and Decoy to 0.009 for LHK, PHD, and fighting chicken; Theta-w, 2.86 for KP to 8.30 for fighting chickens) than in NK-W, NK-B, and DT (pi, 0.001 for NK-W, NK-B, and DT; Theta-w, 0.35 for NK-B to 0.71 for NK-W) (Table 3). With regard to the red junglefowl populations, all populations excluding Si Sa Ket and Petchaburi exhibited similar levels of genetic diversity. The Petchaburi population had two haplotypes, and the genetic diversity was relatively low. The low genetic diversity in the Si Sa Ket population was attributed to the fact that all 30 examined individuals shared only one haplotype. Seven out of the 10 indigenous chicken breeds (LHK, CH, PHD, Decoy, fighting chicken, NK-W, and DT) exhibited negative Tajima’s D values, suggesting that the chickens were bred under purifying selection within each population; however, even though the Tajima’s D values of all populations were not statistically significant (p > 0.05).
Table 3 Genetic diversity of indigenous chicken breeds and red junglefowl populations estimated using of mtDNA D-loop sequences and 28 microsatellite markers.
Full size table

Phylogenetic relationships among mtDNA D-loop haplotypes in red junglefowl in Asia
Haplogroups A, B, CD, E, and J were frequently identified in red junglefowl in Thailand (Fig. 1; Supplementary Table S2). Haplotypes from Thailand, Vietnam, Laos, and Myanmar were located in internal nodes of the haplogroup A, B, CD, and E, and haplotypes from China were derived from the haplotypes in Southeast Asia (Fig. 3). Haplogroup J exclusively consisted of haplotypes from Thailand, Vietnam, and Cambodia. In haplogroup F, haplotypes from Cambodia exhibited the ancestral haplotypes of Chinese red junglefowl. Three haplotypes of red junglefowl from Indonesia were observed in haplogroup K. A small number of haplotypes from the other rare haplogroups G and W to Z were only observed in Chinese red junglefowl.
Figure 3

Median-joining haplotype network of mtDNA D-loop sequences of red junglefowl. The haplotypes are approximately subdivided into 12 haplogroups, A, B, CD, E, F, G, J, K, W, X, Y, and Z, and unclassified haplotypes (U) in this network, according to the haplotype classification by Miao et al.5 and the present study. The sizes of circles indicate relative frequencies of haplotypes, and the number of bars on the lines, which link haplotypes, represent the number of nucleotide substitutions that occurred between the haplotypes for comparison. Black nodes are the inferred intermediate haplotypes. The geographic origins of haplotypes or subspecies names are shown using circles with different colours. The numbers in parentheses after the location or subspecies names indicate the numbers of sequences used for analyses. Detailed information on the sequences obtained from the database are listed in Supplementary Table S7.

Full size image

Genetic characteristics of indigenous chicken breeds and red junglefowl estimated by 28 microsatellite DNA markers
Two-hundred and ninety-eight red junglefowls in two subspecies (G. g. gallus and G. g. spadiceus) from 12 populations and 138 chickens from 10 indigenous chicken breeds, were used for the genetic diversity analyses using 28 microsatellite markers (Table 1; Supplementary Table S3). The allelic richness (AR) values ranged from 1.40 for MCW103 to 1.93 for MCW0014 (1.77 on average) (Supplementary Table S3). Na ranged from 2.14 for MCW0103 to 7.59 for LEI0192 (2.86 on average). FIS varied from – 0.08 for LEI0166 to 0.54 for MCW0014 (0.05 on average). The FST and FIT values fell within the 0.07 (MCW0098) to 0.31 (MCW0247) range and 0.09 (MCW0123) to 0.63 (MCW0014) range, respectively (FST = 0.17, FIT = 0.22 on average). Two markers, MCW0222 and MCW0014, showed the null allele frequency across all populations (NAF) higher than 0.2 (Supplementary Table S3). Looking at each population, null allele frequencies higher than 0.2 were detected in seven, three, and three populations for MCW0014, MCW0222, and LEI0192, respectively, and in less than one or two populations for the other 11 markers (Supplementary Table S4).
Significant departures from Hardy–Weinberg equilibrium were observed for LEI0234, MCW0014, and MCW0123 in more than 10 populations (p  More

Ice environment
Physicochemical analyses of individual ice blocks were conducted to observe eventual differences that could be attributed to spatially related gradual freezing–melting and fresh ice deposition, and to characterize the habitat that enables long-term survival of ice microbiota. All ice samples contained low concentrations of salts, indicating that they originated from recent clean snow. Concentrations of anions in the upper layers, Ice-1 and Ice-2, were similar. However, the bottom layer Ice-3 had distinctly higher electrical conductivity (EC), hardness and alkalinity, less nitrate, and more sulphate. This could indicate that this ice stratum includes a higher proportion of percolation water, which contains more ions than rain and snow as shown by the differences between the percolation water from the cave Planinska jama (that was used for preparing growth media) and the ice, as shown in Table 1. Total organic carbon (TOC) concentrations in the ice were in a range typical of karst streams22, and above the minimum values reported for surface streams, i.e. 0.1–36.6 mg/l23, indicating a significant input of organic matter for the underground ecosystem. TOC indicates an available in situ source of carbon for the ice microbiome. Nitrogen expressed as nitrate did not exhibit high values in ice samples (Table 1). In this respect, a parallel can be drawn with karst sediments, where microbes are commonly limited more by carbon and phosphorus than by nitrogen24.
Table 1 Characteristics of ice samples from Paradana.
Full size table

Besides EC and temperature, pH and dissolved oxygen are additionaly two influential parametres that can affect the abundance and taxonomic structure of microbial communities. pH was found to drive the shift in the community structure not only in habitats such as freshwater, marine sediments or soils but also in cold habitats as Antarctic soils25. In the current samples, the pH effect on the microbial community structure is less evident because all the values are rather similar (Table 1). Cave ice habitats with incoming waterflow are probably not oxygen depleted; on the contrary, for example in Antarctic lakes, glacial meltwater inflow is responsible for oxygen supersaturation26.
Isotopically, the Ice-3 stratum was significantly lighter than the stratum represented by Ice-1 and Ice-2 (Table 1). Correlation of δ2H and deuterium excess did not indicate any effect of kinetic fractionation during water freezing. Thus, intersection of the freezing-line determined by stable isotopes in samples Ice-1 to Ice-3 (δ2H = 6.48δ18O + 2.88) with the local meteoric-water line (LMWL) constructed for the precipitation station at Postojna (Supplementary Fig. S1) (δ2H = 7.95 δ18O + 12.13), provided the δ18O value − 6.3‰ for the original water before freezing. It represents relatively enriched water, but such a value is not uncommon in daily precipitation in Slovenia27. The ice lake in Paradana is presumably formed by the refreezing of water from melting snow accumulated during the winter months20, with some contribution of water dripping from the cave ceiling. November and December 2015 had only a few days with precipitation in Postojna (5 and 4, respectively). However, January and February 2016 had 12 and 20 days with precipitation and monthly totals were high, 152 mm and 312 mm, respectively. The air temperature data adjusted for the elevation difference between Postojna and the Trnovski gozd karst plateau (about 600 m) indicate that about one third of the precipitation in January and one half in February probably fell as snow. The rest was probably a mixture of solid and liquid precipitation, but heavy rains could have occurred as well (e.g. about 55.5 mm of precipitation was measured in Postojna on February 8–9, with mean daily air temperatures between 8 °C and 9 °C). Isotopic composition of precipitation varied significantly between and also during individual events. It is known that snow cover can preserve the isotopic composition of the original snowfalls for long periods28. However, individual snowfalls can mix at the entrance of the cave and the isotopic composition of snow accumulated in the cave can also be influenced by thaws caused by temporary increases of air temperature or rainfall. The isotopic composition of snowmelt water that eventually refreezes in the cave is therefore the result of many processes. Further research with better temporal and spatial resolution of samples and sampling of snowmelt water would be needed to improve knowledge on the dynamics and sources of ice formation. LMWLs known from the literature for other precipitation stations in Slovenia, i.e. Kozina, Portorož and Ljubljana that are given in Supplementary Fig. S1 provided δ18O values for the original water, which we consider too high (− 3.0‰ for LMWL from Portorož, − 3.8‰ for LMWL from Ljubljana and − 5,1‰ for LMWL from Kozina). Postojna is the closest precipitation station to the Paradana and the data on isotopic composition of precipitation cover the period of ice sampling (Supplementary Fig. S2). Therefore, the LMWL at Postojna could be the best representation of the isotopic composition of precipitation supplying water to the Paradana Ice Cave (after considering the elevation difference between the two sites, which is about 600 m).
When analysed in more detail, results obtained using the approach described above (to calculate the isotopic composition of the water that formed the sampled ice) also revealed the sensitivity of the constructed LMWL, the length of data series and extreme values. This is illustrated by records of isotopically very light precipitation in November and December 2015 (δ18O − 17.6‰ and − 14.2 δ18O, respectively). Although such isotopically light precipitation occurred in just two of the 27 months of the observation period, the two values changed the LMWL intercept significantly. However, because they did occur, they cannot be disregarded in the LMWL construction. Daily precipitation data indicate that in both cases monthly values were influenced dominantly by precipitation that fell during just one day (precipitation on those days represented almost the entire monthly precipitation). The LMWL intercept at Postojna without those two months would be 8.3, i.e. closely similar to values in Ljubljana and Kozina. Long-term data from Ljubljana show that the δ18O value of monthly precipitation was lower than − 16.0‰ (values around − 14.0‰ were quite abundant until 1986 and after 2004) in only 5 months in the years 1981–2010. Thus, precipitation with notable isotopically light values, as observed in Postojna between 21 and 23 November 2015 (92% of the precipitation fell on 22 November) appears to be rare in the study area. Nevertheless, it was observed, and it influenced the intercept of LMWL significantly.
It is worth noting that the δ18O values of Ice-1 and Ice-2 are higher than those reported for the Paradana Ice Cave by Carey et al.20. Deuterium excess is also significantly higher than the mean value reported for samples from different depths of ice by Carey et al.20. The difference in δ18O values could be related to different sampling sites. Carey et al.20 sampled the wall ice, whereas the samples collected during this study represent the frozen lake. Investigation of the difference in deuterium levels would be especially interesting. It could point at the input (either by overland flow from the cave entrance or by percolation from the vadose zone) of water from the autumn/winter months, with precipitation from the Eastern Mediterranean air masses having particularly high d-excess (up to 22‰). The Western Mediterranean air masses have d-excess of about 14‰, whereas air masses from the Atlantic have values of only about 10‰29. Late autumn to early winter precipitation in Slovenia (October to December) regularly exhibits high d-excess27. Unfortunately, the available data are insufficient to support analysis of the reason for high deuterium excess of the ice in detail. Study samples also display far lower concentrations of chloride, sulphate and nitrate than samples collected by Carey et al.20.
Concentration of microbes in cave ice
The upper ice stratum represented by Ice-1 and Ice-2 had comparable microbial load expressed in total ATP concentration and total cell counts, whereas the Ice-3 block exhibited significantly higher values (Table 1). Interestingly, the total cell counts of microorganisms in the ice samples was similar (4.67 × 104–15.15 × 104) to that recorded in the Pivka River (SW Slovenia) at the ponor connecting to the karst underground, i.e. 4.29 × 104–12.38 × 104, 30. A large proportion (51.0–85.4%) of entrapped microbes in the ice were viable, showing that they were able to survive ice formation and melting, or even several freezing–melting cycles. A relatively high cell viability can be linked to the availability of compatible solutes, indicated by correspondingly high TOC (Table 1). Not only do sugars and polyols increase microbial resistance to freezing, they can also be used inside the cell as carbon and nitrogen sources31. Higher concentration of salts in Ice-3 block was accompanied by the highest total cell counts and percentage of viable cells (Table 1). In ice from Scărişoara Cave total cell counts varied from 0.84 × 103 to 3.14 × 104 cells/ml with corresponding viability from 28.2 to 84.9%, but no correlation was observed between the ice age (0–13,000 years BP) or depth (0–25 m) and the total number of cells or viability14.
The media types used in this study differed in their ability to stimulate the growth of colonies. In general, nutrient-poor media and low temperatures resulted in higher colony counts in all samples. This phenomenon has been reported previously in cave microbiology, but was not correlated with phylogenetic diversity of microbes obtained on the growth media32. After 28 days of incubation, samples grown on the oligotrophic medium with percolation water (PWA) and cultivated at 10 °C produced the highest colony counts (Table 2). In context this indicates that cave percolation water contains soluble compounds that are not present in tap water and which support the growth of cave-ice microorganisms. With respect to individual samples, the highest colony counts were found in the Ice-3 sample, i.e., 167.37‰ of all cell biomass, determined by flow cytometry (Table 2), and this sample also contained the highest concentration of nutrients (Table 1). Cultivable anaerobic bacteria and fungi were detected in all the ice samples (Table 2).
Table 2 Colony counts (colony-forming units—CFU/ml) and their proportion to total cell counts determined by flow cytometry (‰) at different cultivation conditions and media.
Full size table

Communities in the ice blocks differed in the representation of r-strategists, with their predominance in the Ice-1, and a big difference between Ice-1 and Ice-2, the two ice samples from the same stratum. Interestingly, a more-uniform community structure in terms of r-strategists was displayed in ice block Ice-2–Ice-3 (Table 1). R-strategists commonly dominate in uncrowded and unstable habitats where resources are temporarily abundant and available; with development of a community, r-strategists are gradually replaced by the slow-growing equilibrium K-strategists33.
Cultivation on different media showed that the ice contained metabolically diverse microorganisms, aerobic and anaerobic bacteria and fungi. Two species of yellow-green algae were also recovered in cultures from samples Ice-2 and Ice-3. The two cultivated species, Chloridella glacialis and Ellipsoidion perminimum (for identification see Supplementary Fig. S3), were also found in green ice from Antarctica34. It is known from results of previous studies that algae in ice can survive and even grow under such adverse conditions34,35,36. They can also be well adapted to low light and low water temperature; for example they can thrive under ice- and snow-cover where the available photosynthetic photon flux density is only around the photosynthetic compensation point37. In these terms, and particularly in ice caves with available light, algae and cyanobacteria should not be overlooked as an important part of the ice microbial community. Interestingly, in Himalayan-type glaciers, the algae-rich layers in ice cores were suggested as providing accurate boundary markers of annual layers38. It remains unclear whether algae can be applied similarly as boundary markers in cave ice. Their existence is already known from some caves, for example in Hungary in a small ice cave colonizing surfaces of the ice39, Romania in Scarişoara Ice Cave at the ice/water interface40 and in New Mexico, USA, in Zuni Ice Cave giving the distinctive greenish patina of the layered ice35.
Bacterial community structure
Previous study of ice from the Paradana Ice Cave showed that it probably originates from local rainfall that reaches the cave as drip water after dissolving bedrock while percolating from the surface, and from snow that includes dust particles20. Thus, the largely impacted cave ice in Paradana has different sources, each bringing along a diverse and adaptable microbiota. 16S metagenomic analysis was conducted to describe the taxonomic composition of bacteria found in different ice blocks. Quality filtration of sequence readings gave a total number of 120,381 sequences in the three studied samples (Table 3). The number of operational taxonomic units (OTUs) varied from 185 in Ice-2 to 304 in Ice-1. This pattern was in alignment with values of alpha diversity parameters: extrapolated richness (Chao1), abundance-based coverage estimator (ACE) and Shannon index (Table 3). The rarefaction curves indicated that the diversity had been sampled sufficiently (Supplementary Fig. S4).
Table 3 Number of reads, OTUs, taxon richness and diversity indexes for cave ice samples.
Full size table

A Venn diagram of the distribution of 441 distinct OTUs found in the three studied samples is presented in Fig. 1. Observations showed that 119 OTUs (28.3%) occurred in all three samples and can be interpreted as “a core microbiome”. Three of these OTUs dominated microbial communities in individual samples (relative abundance range 14.5–56.5%) and corresponded to the members of the genera Pseudomonas, Lysobacter, and Sphingomonas, as discussed below. These were followed in abundance by Polaromonas, Flavobacterium, Rhodoferax, Nocardioides, and Pseudonocardia (relative abundance range 3.3–6.9%). Another 35 OTUs had relative abundance above 0.5% and the remaining 76 OTUs had relative abundance below 0.5%. The unique OTUs probably contribute to the variability due to internal variations within the ice block caused by incoming snow or the freezing of percolation water. For example, samples Ice-2 and Ice-3 were cut from the same ice block in a vertical ice profile, but differed in their content of dark, particulate, organic inclusions.
Figure 1

Prokaryotic OTU distribution in cave ice. The Venn diagram indicates the number of distinct and shared OTUs in ice samples Ice-1, Ice-2 and Ice-3.

Full size image

Members of 29 bacterial phyla were detected in the cave ice microbiome (Fig. 2, Supplementary Fig. S5). All samples were dominated by Proteobacteria, with relative abundances of 79.1% in Ice-2, 65.5% in Ice-3 and 55.9% in Ice-1.
Figure 2

Relative abundance of phyla in the cave-ice samples. Phyla with relative abundance  1% of phylotypes in at least one sample and corresponded to Firmicutes, Cyanobacteria and Gemmatimonadetes. Phototrophic bacterial phylotypes belonging to Cyanobacteria were recovered from all three samples. They represented 1.3% of phylotypes in sample Ice-1, but only 0.6% and 0.3% in samples Ice-2 and Ice-3 respectively, from where algae, C. glacialis and E. perminimum, were obtained via cultivation.
Phyla whose relative abundance was less than 1% were grouped together and classified as “Rare phyla”. These phyla comprised 2.2%, 1.5% and 1.2% of Ice-1, Ice-2, and Ice-3, respectively. Their relative abundance is presented in Supplementary Fig. S5.
Among the 31 classes detected in this study, members of Gammaproteobacteria were most abundant and represented 20.1% (Ice-1), 45.3% (Ice-2) and 42.5% (Ice-3) of total detected phylotypes (Fig. 3A). This proteobacterial group was also most abundant in the ice from Scărişoara Cave14. Actinobacteria represented the second most abundant group of phylotypes, with its relative abundances declining from 30.8% in Ice-1 to 26.2% in Ice-3 and 11.7% in Ice-2. Other notably abundant classes were Alpha- and Betaproteobacteria, whose abundances ranged from 9.6 to 26.3% and from 6.9 to 12.3%, respectively.
Figure 3

Heat-map analysis of the relative abundance of members of cave-ice prokaryotic communities at class (A) and genus (B) levels in Ice-1, Ice-2 and Ice-3. Phylotypes whose relative abundances at class level were  More

Construction of RUGs from rumen sequencing data
We produced 979G of Illumina sequencing data from 4 cows, 2 sheep, 4 red deer and 2 reindeer samples, then performed a metagenomic assembly of single samples and a co-assembly of all samples. This created a set of 391 dereplicated genomes (99% ANI (average nucleotide identity)) with estimated completeness ≥ 80% and estimated contamination ≤ 10% (Fig. 1). 284 of these genomes were produced from the single-sample assemblies and 107 were produced from the co-assemblies. 172 genomes were > 90% complete with contamination  90% complete with More

Tested animals
Altogether 73 individuals from seven species of subterranean rodents differing in body mass, phylogenetic relatedness, and sociality were studied (Table 1). All animals were adult non-breeders, or their breeding history was unknown in solitary species, but none of them showed signs of recent breeding, which may theoretically influence measured parameters. For the purpose of this study, we used the following taxa. African mole-rats (Bathyergidae): the social Ansell’s mole-rat Fukomys anselli (Burda, Zima, Scharff, Macholán & Kawalika 1999) occupies the miombo in a small area near Zambia’s capital Lusaka; another social species of the genus Fukomys is named here as Fukomys “Nsanje” because founders of the breeding colony were captured near town Nsanje in south Malawi. Although we used name Fukomys darlingi (Thomas 1895) for mole-rats from this population in previous studies (e.g.38,49), its taxonomic status is still not resolved; the social common mole-rat Cryptomys hottentotus hottentotus (Lesson, 1826) occurs in mesic and semi-arid regions of southern Africa; the solitary Cape dune mole-rat Bathyergus suillus (Schreber, 1782) inhabits sandy soils along the south-western coast of South Africa; and the solitary Cape mole-rat Georychus capensis (Pallas, 1778) occupies mesic areas of the South Africa50. In addition, we studied the social coruro Spalacopus cyanus (Molina, 1782) (Octodontidae) occupying various habitats in Chile51; and the solitary Upper Galilee Mountains blind mole rat Nannospalax galili (Nevo, Ivanitskaya & Beiles 2001) (Spalacidae) from Israel52. Further information about the species including number of individuals used in the study, their physiology and ecology is shown in Table 1.
All experiments were done on captive animals. Georychus capensis, C. hottentotus, and B. suillus, were captured about four months before the experiment, and kept in the animal facility at the University of Pretoria, South Africa (temperature: 23 °C; humidity: 40–60%, photoperiod: 12L:12D). The animals were housed in plastic boxes with wood shavings used as a bedding. Cryptomys hottentotus and G. capensis were fed with sweet potatoes; B. suillus with sweet potatoes, carrots, and fresh grass. Fukomys anselli, F. “Nsanje”, N. galili, and S. cyanus were kept for at least three years in captivity (or born in captivity) before the experiment in the animal facility at the University of South Bohemia in České Budějovice, Czech Republic (temperature: African mole-rats 25 °C, N. galili and S. cyanus 23 °C; humidity: 40–50%, photoperiod: 12L:12D). The animals were kept in terraria with peat as a substrate and fed with carrots, potatoes, sweet potatoes, beetroot, apple, and rodent dry food mix ad libitum.
Experimental design
We measured Tb and Ts in all species at six Tas (10, 15, 20, 25, 30 and 35 °C). Each individual of all species was measured only once in each Ta. Measurements were conducted in temperature controlled experimental rooms in České Budějovice and Pretoria. Each animal was tested on two experimental days.
The animals were placed in the experimental room individually in plastic buckets with wood shavings as bedding. On the first day, the experimental procedure started at Ta 25 °C. They spent 60 min of initial habituation in the first Ta after which Tb and Ts were measured as described in the following paragraphs. The Ta was then increased to 30 °C and 35 °C, respectively. After the experimental room reached the focal Ta, the animals were left minimally 30 min in each Ta to acclimate, and the measurements were repeated. Considering their relatively small body size, tested animals were very likely in thermal equilibrium after this period because mammals of a comparable body mass are thermally equilibrated after similar period of acclimation53,54,55,56. On the second day, the procedure was repeated with the initial Ta 20 °C and decreasing to 15 °C and 10 °C, respectively. The time span between the measurements of the same individual in different Ta was at least 150 min. Between experimental days, the animals were kept at 25 °C in the experimental room (individuals of social species were housed together with their family members).
Body temperature measurements
We used two sets of equipment to measure animal Tb and Ts. In B. suillus, G. capensis, and C. hottentotus, Tb was measured by intraperitoneally injected PIT tags ( More

1.
IPCC Climate Change 2014: Synthesis Report (eds Core Writing Team, Pachauri, R. K. & Meyer L. A.) (IPCC, 2014).
2.
IPCC Special Report on Climate Change, Desertification, Land Degradation, Sustainable Land Management, Food Security, and Greenhouse Gas Fluxes in Terrestrial Ecosystems (IPCC, 2019).

3.
Adoption of the Paris Agreement FCCC/CP/2015/10/Add.1 (UNFCCC, 2015).

4.
Klein Goldewijk, K., Beusen, A., Doelman, J. & Stehfest, E. New anthropogenic land use estimates for the Holocene: HYDE 3.2. Earth Syst. Sci. Data 9, 927–953 (2017).
Article  Google Scholar 

5.
Griscom, B. W. et al. National mitigation potential from natural climate solutions in the tropics. Philos. Trans. R. Soc. B 375, 20190126 (2020).
CAS  Article  Google Scholar 

6.
Friedlingstein, P. et al. Global carbon budget 2019. Earth Syst. Sci. Data 11, 1783–1838 (2019).
Article  Google Scholar 

7.
Houghton, R. A. & Nassikas, A. A. Global and regional fluxes of carbon from land use and land cover change 1850–2015. Glob. Biogeochem. Cycles 31, 456–472 (2017).
CAS  Article  Google Scholar 

8.
Hansis, E., Davis, S. J. & Pongratz, J. Relevance of methodological choices for accounting of land use change carbon fluxes. Glob. Biogeochem. Cycles 29, 1230–1246 (2015).
CAS  Article  Google Scholar 

9.
Pan, Y. et al. A large and persistent carbon sink in the world’s forests. Science 333, 988–993 (2011).
CAS  Article  Google Scholar 

10.
IPCC. 2006 IPCC Guidelines for National Greenhouse Gas Inventories Vol. 4 (eds Eggleston, S. et al.) (IGES, 2006).

11.
IPCC. 2019 Refinement to the 2006 IPCC Guidelines for National Greenhouse Gas Inventories Vol. 4 (eds Buendia, E. C. et al.) (IPCC, 2019).

12.
Grassi, G. et al. Reconciling global-model estimates and country reporting of anthropogenic forest CO2 sinks. Nat. Clim. Change 8, 914–920 (2018).
CAS  Article  Google Scholar 

13.
Lee, D., Llopis, P., Waterworth, R., Roberts, G. & Pearson, T. Approaches to REDD+ Nesting: Lessons Learned from Country Experiences (World Bank, 2018).

14.
Streck, C. et al. Options for Enhancing REDD+ Collaboration in the Context of Article 6 of the Paris Agreement (Meridian Institute, 2017).

15.
Hansen, M. C. et al. High-resolution global maps of 21st-century forest cover change. Science 342, 850–853 (2013).
CAS  Article  Google Scholar 

16.
World Database on Protected Areas User Manual (UNEP, 2016); https://www.protectedplanet.net/en/resources/wdpa-manual

17.
Curtis, P. G., Slay, C. M., Harris, N. L., Tyukavina, A. & Hansen, M. C. Classifying drivers of global forest loss. Science 361, 1108–1111 (2018).
CAS  Article  Google Scholar 

18.
Saatchi, S. S. et al. Benchmark map of forest carbon stocks in tropical regions across three continents. Proc. Natl Acad. Sci. USA 108, 9899–9904 (2011).
CAS  Article  Google Scholar 

19.
PRODES Deforestation (INPE, 2019); http://www.obt.inpe.br/OBT/assuntos/programas/amazonia/prodes

20.
Ogle, S. M.et al. Delineating managed land for reporting national greenhouse gas emissions and removals to the United Nations framework convention on climate change. Carbon Balance Manag. 13, 9 (2018).

21.
Pearson, T. R., Brown, S., Murray, L. & Sidman, G. Greenhouse gas emissions from tropical forest degradation: an underestimated source. Carbon Balance Manag. 12, 3 (2017).
Article  Google Scholar 

22.
Ahlström, A. et al. The dominant role of semi-arid ecosystems in the trend and variability of the land CO2 sink. Science 348, 895–899 (2015).
Article  Google Scholar 

23.
Van Der Werf, G. R. et al. Global fire emissions estimates during 1997–2016. Earth Syst. Sci. Data 9, 697–720 (2017).
Article  Google Scholar 

24.
Kirschbaum, M. U., Zeng, G., Ximenes, F., Giltrap, D. L. & Zeldis, J. R. Towards a more complete quantification of the global carbon cycle. Biogeosciences 16, 831–846 (2019).

25.
Global Forest Observations Initiative. Integration of Remote-sensing and Ground-based Observations for Estimation of Emissions and Removals of Greenhouse Gases in Forests 2nd edn (FAO, 2016).

26.
Potapov, P. et al. Annual continuous fields of woody vegetation structure in the Lower Mekong region from 2000–2017 Landsat time-series. Remote Sens. Environ. 232, 111278 (2019).
Article  Google Scholar 

27.
Federici, S., Lee, D. & Herold, M. Forest Mitigation: A Permanent Contribution to the Paris Agreement? (Climate and Land Use Alliance, 2017).

28.
Romijn, E. et al. Assessing change in national forest monitoring capacities of 99 tropical countries. Ecol. Manag. 352, 109–123 (2015).
Article  Google Scholar 

29.
Cook-Patton, S. Mapping potential carbon capture from global natural forest regrowth. Nature 585, 545–550 (2020).
CAS  Article  Google Scholar 

30.
The Global Stocktake (UNFCCC, 2015); https://unfccc.int/topics/science/workstreams/global-stocktake-referred-to-in-article-14-of-the-paris-agreement

31.
Austin, K. et al. Shifting patterns of oil palm driven deforestation in Indonesia and implications for zero-deforestation commitments. Land Use Policy 69, 41–48 (2017).
Article  Google Scholar 

32.
Gaveau, D. L. et al. Four decades of forest persistence, clearance and logging on Borneo. PLoS ONE 9, e101654 (2014).
Article  Google Scholar 

33.
Miettinen, J., Shi, C. & Liew, S. C. Land cover distribution in the peatlands of Peninsular Malaysia, Sumatra and Borneo in 2015 with changes since 1990. Glob. Ecol. Conserv. 6, 67–78 (2016).
Article  Google Scholar 

34.
Gunarso, P., Hartoyo, M., Agus, F. & Killeen, T. in Reports from the Technical Panels of the 2nd Greenhouse Gas Working Group of the Roundtable on Sustainable Palm Oil (eds Killeen, T. J. & Goon, J.) 29–64 (RSPO, 2013).

35.
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 

36.
Simard, M. et al. Mangrove canopy height globally related to precipitation, temperature and cyclone frequency. Nat. Geosci. 12, 40–45 (2019).
CAS  Article  Google Scholar 

37.
Baccini, A. et al. Estimated carbon dioxide emissions from tropical deforestation improved by carbon-density maps. Nat. Clim. Change 2, 182–185 (2012).
CAS  Article  Google Scholar 

38.
Zarin, D. J. et al. Can carbon emissions from tropical deforestation drop by 50% in 5 years? Glob. Change Biol. 22, 1336–1347 (2016).
Article  Google Scholar 

39.
Mokany, K., Raison, R. J. & Prokushkin, A. S. Critical analysis of root: shoot ratios in terrestrial biomes. Glob. Change Biol. 12, 84–96 (2006).
Article  Google Scholar 

40.
IPCC Supplement to the 2006 IPCC Guidelines for National Greenhouse Gas Inventories: Wetlands (eds Hiraishi, T. et al.) (IPCC, 2014).

41.
Methodological Tool: Estimation of Carbon Stocks and Change in Carbon Stocks in Dead Wood and Litter in A/R CDM Project Activities (UNFCCC, 2013); https://cdm.unfccc.int/methodologies/ARmethodologies/tools/ar-am-tool-12-v3.0.pdf

42.
Hengl, T. et al. SoilGrids250m: global gridded soil information based on machine learning. PLoS ONE 12, e0169748 (2017).
Article  Google Scholar 

43.
Sanderman, J. et al. A global map of mangrove forest soil carbon at 30 m spatial resolution. Environ. Res. Lett. 13, 055002 (2018).
Article  Google Scholar 

44.
Giglio, L., Boschetti, L., Roy, D. P., Humber, M. L. & Justice, C. O. The Collection 6 MODIS burned area mapping algorithm and product. Remote Sens. Environ. 217, 72–85 (2018).
Article  Google Scholar 

45.
Global Ecological Zones for FAO Forest Reporting: 2010 Update (FAO, 2012).

46.
Brus, D. et al. Statistical mapping of tree species over Europe. Eur. J. Res. 131, 145–157 (2012).
Article  Google Scholar 

47.
Del Lungo, A., Ball, J. & Carle, J. Global Planted Forests Thematic Study: Results and Analysis (FAO, 2006); http://www.fao.org/forestry/12139-03441d093f070ea7d7c4e3ec3f306507.pdf

48.
Portugal National Greenhouse Gas Inventory submitted to the UNFCCC, 1990–2018 (UNFCCC, 2020).

49.
Harris, N. L., Goldman, E. D. & Gibbes, S. Spatial Database on Planted Trees Version 1.0 https://www.wri.org/publication/spatialdatabase-planted-trees (WRI, 2019).

50.
Smith, J. E., Heath, L. S., Skog, K. E. & Birdsey, R. A. Methods for Calculating Forest Ecosystem and Harvested Carbon with Standard Estimates for Forest Types of the United States General Technical Report (USDA, Forest Service, 2006); https://doi.org/10.2737/NE-GTR-343

51.
Ruefenacht, B. et al. Conterminous US and Alaska forest type mapping using forest inventory and analysis data. Photogramm. Eng. Remote Sensing 74, 1379–1388 (2008).
Article  Google Scholar 

52.
Pan, Y. et al. Age structure and disturbance legacy of North American forests. Biogeosciences 8, 715–732 (2011) .
Article  Google Scholar 

53.
Turubanova, S., Potapov, P. V., Tyukavina, A. & Hansen, M. C. Ongoing primary forest loss in Brazil, Democratic Republic of the Congo, and Indonesia. Environ. Res. Lett. 13, 074028 (2018).
Article  Google Scholar 

54.
Potapov, P. et al. The last frontiers of wilderness: tracking loss of intact forest landscapes from 2000 to 2013. Sci. Adv. 3, e1600821 (2017).
Article  Google Scholar 

55.
Roman-Cuesta, R. M. et al. Hotspots of gross emissions from the land use sector: patterns, uncertainties, and leading emission sources for the period 2000–2005 in the tropics. Biogeosciences 13, 4253–4269 (2016).
CAS  Article  Google Scholar 

56.
Carter, S. et al. Agriculture-driven deforestation in the tropics from 1990–2015: emissions, trends and uncertainties. Environ. Res. Lett. 13, 014002 (2017).
Article  Google Scholar  More